generalContactSpheresTest.py
You can view and download this file on Github: generalContactSpheresTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: test with parallel computation and particles
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-11-01
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 *
15
16import numpy as np
17import time
18
19useGraphics = True #without test
20#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
21#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
22try: #only if called from test suite
23 from modelUnitTests import exudynTestGlobals #for globally storing test results
24 useGraphics = exudynTestGlobals.useGraphics
25except:
26 class ExudynTestGlobals:
27 pass
28 exudynTestGlobals = ExudynTestGlobals()
29 exudynTestGlobals.isPerformanceTest = False
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31
32SC = exu.SystemContainer()
33mbs = SC.AddSystem()
34
35nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
36
37np.random.seed(1) #always get same results
38
39
40isPerformanceTest = exudynTestGlobals.isPerformanceTest
41# useGraphics = False
42# isPerformanceTest = True
43
44L = 1
45n = 500
46row = 8 #number of spheres in plane
47
48a = L
49vInit= -20
50
51if isPerformanceTest:
52 n *= 100
53 a *= 0.5
54 row *= 2
55 vInit *= 10
56
57radius = 0.5*a
58m = 0.05
59k = 4e4
60d = 0.001*k*4*0.5*0.2 *0.25
61markerList = []
62radiusList = []
63gDataList = []
64
65
66rb = 30*L
67H = 8*L
68Hy=3*L
69pos0 = [0,-rb-0.5*H,0]
70pos1 = [-rb-H,-Hy,0]
71pos2 = [ rb+H,-Hy,0]
72pos3 = [ 0,-Hy,rb+H]
73pos4 = [ 0,-Hy,-rb-H]
74posList=[pos0,pos1,pos2,pos3,pos4]
75for pos in posList:
76 colBG = color4grey
77 colBG[3] = 0.05
78 gDataList += [GraphicsDataSphere(point=pos, radius=rb, color= colBG, nTiles=50)]
79 nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos,
80 visualization=VNodePointGround(show=False)))
81
82 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
83 markerList += [mThis]
84 radiusList += [rb]
85
86
87ns = 20
88gDataSphere = []
89for i in range(ns):
90 gRad = radius*(0.75+0.4*(i/ns))
91 gSphere = GraphicsDataSphere(point=[0,0,0], radius=gRad, color=color4blue, nTiles=5)
92 gDataSphere += [[gSphere]]
93
94gDataSphere = []
95
96timeCreateStart= -time.time()
97
98color4node = color4blue
99maxY = 0
100for i in range(n):
101
102 kk = int(i/int(n/16))
103 color4node = color4list[min(kk%9,9)]
104
105 if (i%20000 == 0): exu.Print("create mass",i)
106 offy = 0
107
108 iy = int(i/(row*row))
109 ix = i%row
110 iz = int(i/row)%row
111
112 if iy % 2 == 1:
113 ix+=0.5
114 iz+=0.5
115
116 offy = -0.25*H-1.5*a+iy*a*0.74 #0.70x is limit value!
117 offx = -0.6*a-H*0.5 + (ix+1)*a
118 offz = -0.6*a-H*0.5 + (iz+1)*a
119
120 valueRand = np.random.random(1)[0]
121 rFact = 0.2 #random part
122 gRad = radius*(1-rFact+rFact*valueRand)
123 pRef = np.array([offx,offy,offz])
124 maxY = max(maxY,offy)
125 nMass = mbs.AddNode(NodePoint(referenceCoordinates=pRef,
126 initialVelocities=[0,vInit,0],
127 visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
128 if i==row*int(row/4)-int(row/2):
129 sNodeNum = nMass
130 if useGraphics:
131 sNode=mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/generalContactSpheres.txt',
132 outputVariableType=exu.OutputVariableType.Position))
133
134 oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
135 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
136 mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
137 markerList += [mThis]
138 radiusList += [gRad]
139
140 mLast = mThis
141
142#put here, such that it is transparent in background
143oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
144 visualization=VObjectGround(graphicsData=gDataList)))
145
146exu.Print('generalContactSpheresTest: create bodies:',timeCreateStart+time.time(),'seconds')
147timeCreateStart= -time.time()
148
149if True:
150 gContact = mbs.AddGeneralContact()
151 gContact.verboseMode = 1
152
153 for i in range(len(markerList)):
154 m = markerList[i]
155 r = radiusList[i]
156 gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
157
158 ssx = 20 #search tree size
159 ssy = 20
160 if isPerformanceTest:
161 ssy*=4
162 # mbs.Assemble()
163 # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,ssx]), frictionPairingsInit=np.eye(1),
164 # searchTreeBoxMin=np.array([-1.2*H,-H,-1.2*H]), searchTreeBoxMax=np.array([1.2*H,14*H,1.2*H]) #80000 particles
165 # )
166
167 gContact.SetFrictionPairings(0.*np.eye(1))
168 gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssx])
169 gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-H,-1.2*H]), pMax=np.array([1.2*H,maxY,1.2*H]))
170
171 exu.Print('treesize=',ssx*ssx*ssy)
172
173exu.Print('generalContactSpheresTest: gContact:',timeCreateStart+time.time(),'seconds')
174
175mbs.Assemble()
176exu.Print("finish gContact")
177
178tEnd = 0.1
179if isPerformanceTest: tEnd *= 0.5
180h= 0.0002
181simulationSettings = exu.SimulationSettings()
182simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
183#simulationSettings.solutionSettings.writeSolutionToFile = True
184simulationSettings.solutionSettings.writeSolutionToFile = True
185simulationSettings.solutionSettings.solutionWritePeriod = 0.02
186simulationSettings.solutionSettings.sensorsWritePeriod = h*10
187simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
188simulationSettings.solutionSettings.exportAccelerations = False
189simulationSettings.solutionSettings.exportVelocities = False
190simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
191simulationSettings.displayComputationTime = True
192#simulationSettings.displayStatistics = True
193simulationSettings.timeIntegration.verboseMode = 1
194simulationSettings.parallel.numberOfThreads = 1 #use 1 thread to create reproducible results (due to round off errors in sparse vector?)
195if isPerformanceTest: simulationSettings.parallel.numberOfThreads = 8
196
197simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
198simulationSettings.timeIntegration.newton.useModifiedNewton = False
199
200SC.visualizationSettings.general.graphicsUpdateInterval=0.5
201SC.visualizationSettings.general.circleTiling=200
202SC.visualizationSettings.general.drawCoordinateSystem=False
203SC.visualizationSettings.loads.show=False
204SC.visualizationSettings.bodies.show=True
205SC.visualizationSettings.markers.show=False
206
207SC.visualizationSettings.nodes.show=True
208SC.visualizationSettings.nodes.drawNodesAsPoint = False
209SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
210SC.visualizationSettings.nodes.tiling = 4
211
212SC.visualizationSettings.window.renderWindowSize=[800,800]
213#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
214SC.visualizationSettings.openGL.multiSampling = 4
215#improved OpenGL rendering
216
217
218if useGraphics:
219 SC.visualizationSettings.general.autoFitScene = False
220 exu.StartRenderer()
221 if 'renderState' in exu.sys:
222 SC.SetRenderState(exu.sys['renderState'])
223 mbs.WaitForUserToContinue()
224
225simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
226simulationSettings.timeIntegration.endTime = tEnd
227simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
228simulationSettings.timeIntegration.explicitIntegration.computeMassMatrixInversePerBody = True ##2022-12-16: increase performance for multi-threading, Newton increment faster by factor 6 for 8 threads
229
230mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
231
232u = mbs.GetNodeOutput(sNodeNum, exu.OutputVariableType.Coordinates)
233uSum = u[0] + u[1] + u[2]
234exu.Print("u =", u)
235exu.Print('solution of generalContactSpheresTest=',uSum)
236
237if isPerformanceTest:
238 exudynTestGlobals.testError = uSum - (-5.946497644233068)
239else:
240 exudynTestGlobals.testError = uSum - (-1.0947542400425323)
241
242exudynTestGlobals.testResult = uSum
243
244
245if useGraphics:
246 SC.WaitForRenderEngineStopFlag()
247 exu.StopRenderer() #safely close rendering window!
248
249if useGraphics:
250
251 # mbs.PlotSensor([sNode], [2])
252 mbs.PlotSensor([sNode,sNode], [0,1])