createSphereQuadContact.py
You can view and download this file on Github: createSphereQuadContact.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for CreateSphereQuadContact, combining sphere-sphere and 2 sphere-triangle contacts
5#
6# Author: Johannes Gerstmayr
7# Date: 2025-06-29
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 *
15from math import sin, cos
16import exudyn.graphics as graphics
17import numpy as np
18
19testIsActive = exu.sys.get('testIsActive', False)
20#testIsActive = True
21
22testSolution = 0
23
24SC = exu.SystemContainer()
25mbs = SC.AddSystem()
26
27methodList = ['trigs','quad'] #0=trigs, 1=quad; both tests should gave same results!
28listSolutions = []
29
30for methodNum, method in enumerate(methodList):
31 mbs.Reset()
32
33 exu.Print('\n\n***********************************')
34 exu.Print('*** Test contact create method: '+method+' ***')
35 exu.Print('***********************************\n')
36
37 radius=0.1
38 mass = 0.2 #mass in kg
39 contactStiffness = 2e4 #stiffness of spring-damper in N/m
40 contactDamping = 0*5e-4*contactStiffness #damping constant in N/(m/s)
41 dynamicFriction = 0.2
42 restitutionCoefficient = 0.75
43 impactModel = 2
44
45 isExplicitSolver = False
46 tEnd = 0.65 #end time of simulation
47 stepSize = 2e-4 #*10
48
49 g = 9.81
50
51 size = 1
52
53
54 oGround = mbs.CreateGround(graphicsDataList=[graphics.CheckerBoard(point=[0,0,0],size=2*size,
55 color=graphics.color.lightgrey[0:3]+[graphics.material.indexChrome],
56 alternatingColor=graphics.color.lightgrey2[0:3]+[graphics.material.indexChrome],
57 ), ])
58 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
59 mGround2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0,0,-radius]))
60
61 listMasses=[]
62 sPosList=[]
63
64 ny = 4
65 cnt = -1
66 for jy in range(ny):
67
68 for ix in range(max(1,jy)):
69 cnt+=1
70 x = (ix-(jy-1)*0.5)*2*radius
71 y = -4*radius + jy*radius*np.sqrt(3)
72
73 vy = 2*sin(cnt*pi/3)
74 vx = 2*cos(cnt*pi/3)
75
76 massFact = 1
77
78 #for explicit solver, we need Lie group node:
79 nodeType=exu.NodeType.RotationRotationVector if isExplicitSolver else exu.NodeType.RotationEulerParameters
80
81 oMass = mbs.CreateRigidBody(referencePosition=[x,y,radius],
82 initialVelocity=[vx,vy,0],
83 initialAngularVelocity=[0,0,0],
84 nodeType=nodeType,
85 inertia=InertiaSphere(mass=massFact*mass, radius=radius),
86 gravity = [0,0,-g],
87 graphicsDataList=[graphics.Sphere(radius=radius,
88 color=graphics.colorList[cnt][0:3]+[graphics.material.indexDefault],
89 nTiles=24)],
90 )
91 listMasses.append(oMass)
92 mMass = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass))
93
94 if False: #object mode
95 for oMass2 in listMasses[:-1]:
96 mMass2 = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass2))
97 nData1 = mbs.AddNode(NodeGenericData(initialCoordinates=[0.1,0,0,0],
98 numberOfDataCoordinates=4))
99 oSSC = mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mMass, mMass2],
100 nodeNumber=nData1,
101 spheresRadii=[radius, radius],
102 contactStiffness = contactStiffness,
103 dynamicFriction=dynamicFriction,
104 impactModel = impactModel,
105 restitutionCoefficient = restitutionCoefficient,
106 visualization=VObjectContactSphereSphere(show=True),
107 ))
108 else:
109 #create contact between each sphere:
110 for oMass2 in listMasses[:-1]:
111 oSSC = mbs.CreateSphereSphereContact(itemNumbers=[oMass, oMass2],
112 spheresRadii=[radius, radius],
113 contactStiffness = contactStiffness,
114 dynamicFriction = dynamicFriction,
115 impactModel = impactModel,
116 restitutionCoefficient = restitutionCoefficient,
117 show=True,
118 )
119
120 if method=='trigs': #create 2 triangles
121 trianglePoints0 = exu.Vector3DList([[-size,-size,0],[size,-size,0],[-size,size,0]])
122 trianglePoints1 = exu.Vector3DList([[size,-size,0],[size,size,0],[-size,size,0]])
123 includeEdgesList = [5,3]
124
125 trigList = [trianglePoints0,trianglePoints1]
126 for k, trianglePoints in enumerate(trigList):
127 oSSC = mbs.CreateSphereTriangleContact(itemNumbers=[oMass, oGround],
128 trianglePoints=trianglePoints,
129 includeEdges=includeEdgesList[k],
130 sphereRadius=radius,
131 contactStiffness = contactStiffness,
132 dynamicFriction=dynamicFriction,
133 impactModel = impactModel,
134 restitutionCoefficient = restitutionCoefficient,
135 show = True
136 )
137 else: #alternative with quads
138 quadPoints = exu.Vector3DList([[-size,-size,0],[size,-size,0],[size,size,0],[-size,size,0]])
139 oSSC = mbs.CreateSphereQuadContact(itemNumbers=[oMass, oGround],
140 quadPoints=quadPoints,
141 includeEdges=15, #all edges
142 sphereRadius=radius,
143 contactStiffness = contactStiffness,
144 dynamicFriction=dynamicFriction,
145 impactModel = impactModel,
146 restitutionCoefficient = restitutionCoefficient,
147 show = True
148 )
149
150
151 sPos=mbs.AddSensor(SensorBody(bodyNumber=oMass, storeInternal=True,
152 outputVariableType=exu.OutputVariableType.Position))
153 sPosList.append(sPos)
154
155 #exu.Print(mbs)
156 mbs.Assemble()
157
158 simulationSettings = exu.SimulationSettings()
159 simulationSettings.solution.file.write = (not testIsActive) #only the SolutionViewer reads it (#2492)
160 simulationSettings.solution.file.writePeriod = 0.005
161 simulationSettings.solution.sensors.writePeriod = 0.001 #output interval
162 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
163 simulationSettings.timeIntegration.endTime = tEnd
164 #simulationSettings.timeIntegration.realtime.active = True
165 simulationSettings.timeIntegration.newton.absoluteTolerance = 1e-6
166 simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
167
168 simulationSettings.timeIntegration.stepInformation = 3 #remove flag 64 which shows step reduction warnings
169
170
171 simulationSettings.timeIntegration.newton.useModifiedNewton = True
172 simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
173
174 simulationSettings.show.statistics = True
175 simulationSettings.timeIntegration.verboseMode = 1
176 SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
177 SC.visualizationSettings.general.showSolverInformation = False
178
179
180 if not testIsActive:
181 SC.renderer.Start() #start graphics visualization
182 if methodNum == 0:
183 SC.renderer.DoIdleTasks() #wait for pressing SPACE bar to continue
184
185 mbs.SolveDynamic(simulationSettings)
186
187 if not testIsActive:
188 #SC.renderer.DoIdleTasks()
189 SC.renderer.Stop() #safely close rendering window!
190
191 if False:
192 mbs.PlotSensor(sPosList, components=[2]*len(sPosList))
193
194 #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
195 ode2 = mbs.systemData.GetODE2Coordinates()
196 listSolutions.append(ode2)
197
198 testSolution += 0.1*np.linalg.norm(ode2)
199
200diff = listSolutions[0] - listSolutions[1]
201exu.Print('solution diff=',np.linalg.norm(diff))
202
203for i, sol in enumerate(listSolutions):
204 exu.Print('solver=',str(methodList[i]),'\nsol=',sol[0:6])
205
206#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
207exu.Print('solution of createSphereQuadContact=',testSolution)
208exu.sys['testResult'] = testSolution
209#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
210
211
212if not testIsActive:
213 mbs.SolutionViewer()