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