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