sphereTriangleTest2.py

You can view and download this file on Github: sphereTriangleTest2.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test for SphereTrigContact, combining sphere-sphere and sphere-triangle contact
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2025-06-14
  8#
  9# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
 10#
 11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 12
 13import exudyn as exu
 14from exudyn.utilities import *
 15import exudyn.graphics as graphics
 16import numpy as np
 17
 18useGraphics = True #without test
 19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 21try: #only if called from test suite
 22    from modelUnitTests import exudynTestGlobals #for globally storing test results
 23    useGraphics = exudynTestGlobals.useGraphics
 24except:
 25    class ExudynTestGlobals:
 26        pass
 27    exudynTestGlobals = ExudynTestGlobals()
 28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 29#useGraphics = False
 30
 31testSolution = 0
 32
 33SC = exu.SystemContainer()
 34mbs = SC.AddSystem()
 35
 36solverList = [
 37              exu.DynamicSolverType.GeneralizedAlpha,
 38              exu.DynamicSolverType.VelocityVerlet,
 39              #exu.DynamicSolverType.TrapezoidalIndex2, #not used in test
 40              #exu.DynamicSolverType.ExplicitEuler,     #not used in test
 41              ]
 42
 43listSolutions = []
 44
 45for solverNum, solver in enumerate(solverList):
 46    mbs.Reset()
 47
 48    isExplicitSolver = (solver not in [exu.DynamicSolverType.GeneralizedAlpha,
 49                                       exu.DynamicSolverType.TrapezoidalIndex2])
 50
 51    exu.Print('\n\n***********************************')
 52    exu.Print('*** Test solver: '+str(solver)+' ***')
 53    exu.Print('*** is explicit='+str(isExplicitSolver)+' ***')
 54    exu.Print('***********************************\n')
 55    radius=0.1
 56    mass = 0.2                              #mass in kg
 57    contactStiffness = 2e4                  #stiffness of spring-damper in N/m
 58    contactDamping = 0*5e-4*contactStiffness  #damping constant in N/(m/s)
 59    dynamicFriction = 0.2
 60    restitutionCoefficient = 0.75
 61    impactModel = 2
 62
 63    tEnd = 0.25     #end time of simulation
 64    stepSize = 2e-4 #*10
 65    if isExplicitSolver:
 66        stepSize *= 0.1
 67
 68    if solver == exu.DynamicSolverType.RK44:
 69        stepSize *= 5e-6 #requires smaller step size
 70
 71    g = 9.81
 72
 73    size = 1
 74
 75
 76    oGround = mbs.CreateGround(graphicsDataList=[graphics.CheckerBoard(point=[0,0,0],size=2*size,
 77                                                                       color=graphics.color.lightgrey[0:3]+[graphics.material.indexChrome],
 78                                                                       alternatingColor=graphics.color.lightgrey2[0:3]+[graphics.material.indexChrome],
 79                                                                       ), ])
 80    mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
 81    mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,-radius]))
 82
 83    listMasses=[]
 84    sPosList=[]
 85
 86    ny = 4
 87    cnt = -1
 88    for jy in range(ny):
 89
 90        for ix in range(max(1,jy)):
 91            cnt+=1
 92            x = (ix-(jy-1)*0.5)*2*radius
 93            y = -4*radius + jy*radius*np.sqrt(3)
 94
 95            vy = 0
 96            vx = 0
 97            massFact = 1
 98            angX = 0
 99            if cnt == 0:
100                vy = 2
101                vx = 0.1
102                angX = -vy/radius
103                y -= 0.1
104                x -= radius
105                massFact = 2
106
107            #for explicit solver, we need Lie group node:
108            nodeType=exu.NodeType.RotationRotationVector if isExplicitSolver else exu.NodeType.RotationEulerParameters
109
110            oMass = mbs.CreateRigidBody(referencePosition=[x,y,radius],
111                                        initialVelocity=[vx,vy,0],
112                                        initialAngularVelocity=[angX,0,0],
113                                        nodeType=nodeType,
114                                        inertia=InertiaSphere(mass=massFact*mass, radius=radius),
115                                        gravity = [0,0,-g],
116                                        graphicsDataList=[graphics.Sphere(radius=radius,
117                                                                          color=graphics.colorList[cnt][0:3]+[graphics.material.indexDefault],
118                                                                          nTiles=48)],
119                                        )
120            listMasses.append(oMass)
121            mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
122
123            for oMass2 in listMasses[:-1]:
124                mMass2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass2))
125                nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
126                                                    numberOfDataCoordinates=4))
127                oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mMass, mMass2],
128                                                                nodeNumber=nData1,
129                                                                spheresRadii=[radius, radius],
130                                                                contactStiffness = contactStiffness,
131                                                                dynamicFriction=dynamicFriction,
132                                                                impactModel = impactModel,
133                                                                restitutionCoefficient = restitutionCoefficient,
134                                                                visualization=VObjectContactSphereSphere(show=True),
135                                                                ))
136            trianglePoints0 = exu.Vector3DList([[-size,-size,0],[size,-size,0],[-size,size,0]])
137            trianglePoints1 = exu.Vector3DList([[size,-size,0],[size,size,0],[-size,size,0]])
138            includeEdgesList = [5,3]
139
140            trigList = [trianglePoints0,trianglePoints1]
141            for k, trianglePoints in enumerate(trigList):
142                nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
143                                                    numberOfDataCoordinates=4))
144                oSSC = mbs.AddObject(ObjectContactSphereTriangle(markerNumbers=[mMass, mGround],
145                                                                nodeNumber=nData1,
146                                                                trianglePoints=trianglePoints,
147                                                                includeEdges=includeEdgesList[k],
148                                                                radiusSphere=radius,
149                                                                contactStiffness = contactStiffness,
150                                                                dynamicFriction=dynamicFriction,
151                                                                impactModel = impactModel,
152                                                                restitutionCoefficient = restitutionCoefficient,
153                                                                visualization=VObjectContactSphereSphere(show=True),
154                                                                ))
155
156            sPos=mbs.AddSensor(SensorBody(bodyNumber=oMass, storeInternal=True,
157                                          outputVariableType=exu.OutputVariableType.Position))
158            sPosList.append(sPos)
159
160    #exu.Print(mbs)
161    mbs.Assemble()
162
163    simulationSettings = exu.SimulationSettings()
164    simulationSettings.solutionSettings.writeSolutionToFile = True
165    simulationSettings.solutionSettings.solutionWritePeriod = 0.005
166    simulationSettings.solutionSettings.sensorsWritePeriod = 0.001  #output interval
167
168    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
169    simulationSettings.timeIntegration.endTime = tEnd
170    # simulationSettings.timeIntegration.numberOfSteps = 1
171    # simulationSettings.timeIntegration.endTime = stepSize
172    simulationSettings.timeIntegration.verboseMode = 1
173
174    #simulationSettings.timeIntegration.simulateInRealtime = True
175    simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
176    simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
177    #simulationSettings.timeIntegration.generalizedAlpha.computeInitialAccelerations = False
178    simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #speedup
179    simulationSettings.timeIntegration.explicitIntegration.computeMassMatrixInversePerBody = True #speedup
180    simulationSettings.timeIntegration.stepInformation = 3 #remove flag 64 which shows step reduction warnings
181
182    if isExplicitSolver:
183        simulationSettings.timeIntegration.discontinuous.useRecommendedStepSize = False #anyway do fine steps with explicit integrator
184
185    if solver == exu.DynamicSolverType.DOPRI5: #not recommended
186        simulationSettings.timeIntegration.absoluteTolerance = 0.25e-4 #default=1e-8 -> very accurate & small step size
187        simulationSettings.timeIntegration.relativeTolerance = 0.25e-4 #default=1e-8 -> very accurate & small step size
188        simulationSettings.timeIntegration.discontinuous.maxIterations = 1 #not used, as we anyway do step refinement
189        simulationSettings.timeIntegration.discontinuous.iterationTolerance = 1
190
191    simulationSettings.timeIntegration.newton.useModifiedNewton = True
192    simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
193
194    simulationSettings.displayStatistics = True
195    SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
196    SC.visualizationSettings.general.showSolverInformation = False
197
198    #++++++++++++++++++++++++++++++++++++++++++++++++++
199    #special visualization options
200    SC.visualizationSettings.openGL.multiSampling = 2
201    SC.visualizationSettings.openGL.light0.shadow = 0.2
202    SC.visualizationSettings.openGL.advanced.depthSorting = True
203    SC.visualizationSettings.openGL.light0.position = [3, -5, 10.0, 0.0]
204    SC.visualizationSettings.openGL.light1.enable = False
205    SC.visualizationSettings.view0.camera.perspective = 1
206    SC.visualizationSettings.raytracer.numberOfThreads = 16
207    SC.visualizationSettings.raytracer.keepWindowActive = True
208    SC.visualizationSettings.raytracer.imageSizeFactor = 7
209    SC.visualizationSettings.raytracer.maxTransparencyDepth = 2
210    SC.visualizationSettings.raytracer.maxReflectionDepth = 2
211    SC.visualizationSettings.raytracer.advanced.searchTreeFactor = 8
212    SC.visualizationSettings.raytracer.verbose = True
213
214    mat0 = SC.renderer.materials.Get(0)
215    mat0.alpha = 0.3
216    mat0.ior = 1.25
217    mat0.reflectivity = 0.3
218    mat0.shininess = 80
219    mat0.specular = [0.8]*3
220    SC.renderer.materials.Set(0,mat0)
221
222    mat1 = SC.renderer.materials.Get(4)
223    mat1.shininess = 40
224    mat1.specular = [0.8]*3
225    mat1.reflectivity = 0.2
226    SC.renderer.materials.Set(4,mat1)
227
228    SC.visualizationSettings.view0.window.renderWindowSize=[1280,1024]
229    SC.visualizationSettings.nodes.showBasis = True
230    SC.visualizationSettings.nodes.basisSize = radius*1.3
231    SC.visualizationSettings.exportImages.saveImageTimeOut = 500000
232    SC.visualizationSettings.connectors.show = False
233    #++++++++++++++++++++++++++++++++++++++++++++++++++
234
235    if useGraphics:
236        SC.renderer.Start()              #start graphics visualization
237        if solverNum == 0:
238            SC.renderer.DoIdleTasks()    #wait for pressing SPACE bar to continue
239
240    mbs.SolveDynamic(simulationSettings,
241                     solverType=solver)
242
243    if useGraphics:
244        #SC.renderer.DoIdleTasks()
245        SC.renderer.Stop()               #safely close rendering window!
246
247        if False:
248            mbs.PlotSensor(sPosList, components=[2]*len(sPosList))
249
250    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251    ode2 = mbs.systemData.GetODE2Coordinates()
252    listSolutions.append(ode2)
253
254    testSolution += np.linalg.norm(ode2)
255
256#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257exu.Print('solution of sphereTriangleTest2=',testSolution)
258exudynTestGlobals.testResult = testSolution
259#dense:  4.356119232234876 (since V1.10.78)
260#sparse: 4.356119232231812 (since V1.10.78)
261#OLD
262#sparse: 4.356128117693937 (until V1.10.77)
263#dense:  4.356119232234876 (until V1.10.77)
264#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
265
266for i, sol in enumerate(listSolutions):
267    exu.Print('solver=',str(solverList[i]),'\nsol=',sol[0:6])
268
269
270if useGraphics and False:
271    mbs.SolutionViewer()
272
273#convergence analysis:
274#NOTE: y-component is very sensitive to impact; would be better to check velocities
275
276#stepSize = 2e-4 / 2e-5 (implicit/explicit)
277# solver= DynamicSolverType.GeneralizedAlpha
278# sol= [ 2.28555584e-02 -8.46781176e-03 -1.90795166e-04 -4.21511032e-01 -8.08801267e-01  5.80026017e-02]
279# solver= DynamicSolverType.TrapezoidalIndex2
280# sol= [ 2.28534184e-02 -8.41688176e-03 -1.92925390e-04 -4.21674541e-01 -8.08919601e-01  5.79880986e-02]
281# solver= DynamicSolverType.ExplicitEuler
282# sol= [ 2.29470185e-02  3.61558459e-03 -1.95206592e-04 -1.85876084e+00  1.40580855e-01  1.99925961e-01]
283# solver= DynamicSolverType.VelocityVerlet
284# sol= [ 2.29144361e-02  3.49521933e-03 -1.94918274e-04 -1.85823086e+00  1.39871248e-01  2.00701750e-01]
285
286#stepSize = 1e-4 / 1e-5
287# solver= DynamicSolverType.GeneralizedAlpha
288# sol= [ 2.28801345e-02 -6.81190990e-03 -2.49749879e-04 -3.95284914e-01 -7.89223915e-01  5.91060718e-02]
289# solver= DynamicSolverType.TrapezoidalIndex2
290# sol= [ 2.28825196e-02 -6.80461763e-03 -2.51272130e-04 -3.95282395e-01 -7.89223494e-01  5.91315017e-02]
291# solver= DynamicSolverType.ExplicitEuler
292# sol= [ 2.30210743e-02 -5.30821458e-04 -1.96205957e-04 -1.84244869e+00  1.41086785e-01  1.99824875e-01]
293# solver= DynamicSolverType.VelocityVerlet
294# sol= [ 2.29556825e-02  3.24907957e-03 -1.96201003e-04 -1.85622032e+00  1.40718208e-01  1.99712086e-01]
295
296#stepSize = 0.5e-4 / 0.5e-5
297# solver= DynamicSolverType.GeneralizedAlpha
298# sol= [ 2.29061099e-02 -4.88224977e-03 -1.97367604e-04 -3.96792723e-01 -7.90436956e-01  5.93869433e-02]
299# solver= DynamicSolverType.TrapezoidalIndex2
300# sol= [ 2.29065522e-02 -4.89744183e-03 -1.97441120e-04 -3.96759579e-01 -7.90411593e-01  5.93912705e-02]
301# solver= DynamicSolverType.ExplicitEuler
302# sol= [ 2.33552784e-02  5.23855373e-02 -1.96200000e-04 -2.11720759e+00  1.61935043e-01  1.69430071e-01]
303# solver= DynamicSolverType.VelocityVerlet
304# sol= [ 2.29954307e-02  1.38827453e-03 -1.97018837e-04 -1.84873719e+00  1.41131642e-01  1.99464230e-01]
305
306#stepSize = 0.2e-4 / 0.2e-5
307# solver= DynamicSolverType.GeneralizedAlpha
308# sol= [ 2.29758255e-02 -3.22647701e-03 -1.93419162e-04 -3.98165065e-01 -7.91552245e-01  6.01495798e-02]
309# solver= DynamicSolverType.TrapezoidalIndex2
310# sol= [ 2.29758298e-02 -3.22677010e-03 -1.93419056e-04 -3.98164462e-01 -7.91551790e-01  6.01496166e-02]
311# solver= DynamicSolverType.ExplicitEuler
312# sol= [ 2.22822848e-02  4.66325857e-02 -1.96200000e-04 -2.01966751e+00  1.40862670e-01  1.94681099e-01]
313# solver= DynamicSolverType.VelocityVerlet
314# sol= [ 2.28667496e-02  8.44902056e-05 -1.96212507e-04 -1.84355659e+00  1.38051928e-01  2.03132629e-01]
315
316#stepSize = 0.1e-4 / 0.1e-5
317# solver= DynamicSolverType.GeneralizedAlpha
318# sol= [ 2.30115801e-02 -2.49245354e-03 -1.96210497e-04 -3.98770615e-01 -7.92045008e-01  6.05675129e-02]
319# solver= DynamicSolverType.TrapezoidalIndex2
320# sol= [ 2.30115800e-02 -2.49245878e-03 -1.96210497e-04 -3.98770615e-01 -7.92045009e-01  6.05675131e-02]
321# solver= DynamicSolverType.ExplicitEuler
322# sol= [ 2.30288928e-02  1.37441597e-02 -1.96200025e-04 -1.89273830e+00  1.44896879e-01  1.93767400e-01]
323# solver= DynamicSolverType.VelocityVerlet
324# sol= [ 2.30373557e-02  9.72852328e-04 -1.96201848e-04 -1.84638483e+00  1.42049290e-01  1.98406334e-01]