particlesTest3D.py
You can view and download this file on Github: particlesTest3D.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.itemInterface import *
15from exudyn.utilities import *
16from exudyn.graphicsDataUtilities import *
17
18import numpy as np
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23#create an environment for mini example
24
25nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
26#mLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nGround))
27
28np.random.seed(1) #always get same results
29
30useGraphics = True
31
32L = 1
33n = 8000*8*4 #*8*4 #32*8*8
34# n = 5000
35a = 0.2*L*0.5*10*0.5
36radius = 0.35*a
37m = 0.05
38k = 4e3*10*0.5 #4e3 needs h=1e-4
39d = 0.001*k*4*0.5
40markerList = []
41radiusList = []
42gDataList = []
43
44
45rb = 30*L
46H = 8*L
47pos0 = [0,-rb-0.5*H,0]
48pos1 = [-rb-H,0,0]
49pos2 = [ rb+H,0,0]
50pos3 = [ 0,0,rb+H]
51pos4 = [ 0,0,-rb-H]
52posList=[pos0,pos1,pos2,pos3,pos4]
53for pos in posList:
54 #gDataList += [{'type':'Circle','position':pos,'radius':rb, 'color':color4grey}]
55 #gDataList += [GraphicsDataCylinder(pAxis=pos, vAxis=[0,0,0.1], radius=rb, color= color4grey, nTiles=200)]
56 colBG = color4grey
57 colBG[3] = 0.05
58 gDataList += [GraphicsDataSphere(point=pos, radius=rb, color= colBG, nTiles=100)]
59 #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,16*H,color=color4red)]
60 nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos,
61 visualization=VNodePointGround(show=False)))
62 #oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
63 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
64 markerList += [mThis]
65 radiusList += [rb]
66
67
68ns = 20
69gDataSphere = []
70for i in range(ns):
71 gRad = radius*(0.75+0.4*(i/ns))
72 # gSphere = GraphicsDataCylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=color4blue, nTiles=12)
73 # gSphere2 = GraphicsDataCylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=color4steelblue, nTiles=10)
74 gSphere = GraphicsDataSphere(point=[0,0,0], radius=gRad, color=color4blue, nTiles=8)
75 gDataSphere += [[gSphere]]
76
77gDataSphere = []
78
79color4node = color4blue
80print("start create: number of masses =",n)
81for i in range(n):
82
83 kk = int(i/12800)
84 color4node = color4list[min(kk%12,11)]
85 # if (i%10000 == 0):
86 # gDataSphere = []
87 # for i in range(ns):
88 # gRad = radius*(0.75+0.4*(i/ns))
89 # # gSphere = GraphicsDataCylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=color4blue, nTiles=12)
90 # # gSphere2 = GraphicsDataCylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=color4steelblue, nTiles=10)
91 # gSphere = GraphicsDataSphere(point=[0,0,0], radius=gRad, color=color4list[min(k%12,11)], nTiles=8)
92 # gDataSphere += [[gSphere]]
93
94
95 if (i%20000 == 0): print("create mass",i)
96 offy = 0
97 row = 8*2 #160
98 offy = -0.25*H-1.5*a+int(i/(row*row))*a+a*0.2*np.random.random(1)[0]
99
100 offx = -0.6*a-H*0.5 + (i%row+1)*a+0.2*a*np.random.random(1)[0]
101 offz = -0.6*a-H*0.5 + (int(i/row)%row+1)*a+0.2*a*np.random.random(1)[0]
102
103 valueRand = np.random.random(1)[0]
104 gRad = radius*(0.75+0.4*valueRand)
105 #gSphere = GraphicsDataCylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.25], radius=gRad, color= color4steelblue, nTiles=16)
106 #gSphere2 = GraphicsDataCylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.3], radius=0.8*gRad, color= color4blue, nTiles=12)
107 nMass = mbs.AddNode(NodePoint(referenceCoordinates=[offx,offy,offz],
108 initialVelocities=[0,-20,0],
109 visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
110 # gData = gDataSphere[int(valueRand*ns)]
111 # if not useGraphics:
112 # gData = []
113 # if i%2 != 0:
114 # gData = []
115
116 oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
117 #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
118 # visualization=VMassPoint(graphicsData=gData)
119 ))
120 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
121 mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
122 markerList += [mThis]
123 radiusList += [gRad]
124 #if (i==n-1):
125 # mbs.AddLoad(Force(markerNumber = mThis, loadVector = [5, -20, 0]))
126
127 #mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLast, mThis],
128 # stiffness = [k,k,k], damping=[d,d,d], offset=[a,0,0],
129 # visualization = VCartesianSpringDamper(drawSize = 0.1*a)))
130
131 mLast = mThis
132print("finish create")
133#put here, such that it is transparent in background
134oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
135 visualization=VObjectGround(graphicsData=gDataList)))
136
137if True:
138 gContact = mbs.AddGeneralContact()
139 gContact.verboseMode = 1
140
141 for i in range(len(markerList)):
142 m = markerList[i]
143 r = radiusList[i]
144 gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
145
146 # f=n/32000
147 ssx = 20 #search tree size
148 #ssy = int(500*f) #search tree size
149 ssy = 200
150 # mbs.Assemble()
151 # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,ssx]), frictionPairingsInit=np.eye(1),
152 # searchTreeBoxMin=np.array([-1.2*H,-0.75*H,-1.2*H]),
153 # searchTreeBoxMax=np.array([1.2*H,4*16*H,1.2*H])
154 # )
155 gContact.SetFrictionPairings(np.eye(1))
156 gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssx])
157 gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-0.75*H,-1.2*H]), pMax=np.array([1.2*H,4*16*H,1.2*H]))
158 print('treesize=',ssx*ssx*ssy)
159
160mbs.Assemble()
161print("finish gContact")
162
163tEnd = 10
164h= 0.0001*0.25
165simulationSettings = exu.SimulationSettings()
166simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
167#simulationSettings.solutionSettings.writeSolutionToFile = True
168simulationSettings.solutionSettings.writeSolutionToFile = True
169simulationSettings.solutionSettings.solutionWritePeriod = 0.02
170simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
171simulationSettings.solutionSettings.exportAccelerations = False
172simulationSettings.solutionSettings.exportVelocities = False
173#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'particles3D.txt'
174simulationSettings.displayComputationTime = True
175#simulationSettings.displayStatistics = True
176simulationSettings.timeIntegration.verboseMode = 1
177simulationSettings.parallel.numberOfThreads = 4
178
179simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
180simulationSettings.timeIntegration.newton.useModifiedNewton = False
181
182SC.visualizationSettings.general.graphicsUpdateInterval=0.5
183SC.visualizationSettings.general.circleTiling=200
184SC.visualizationSettings.general.drawCoordinateSystem=False
185SC.visualizationSettings.loads.show=False
186SC.visualizationSettings.bodies.show=False
187SC.visualizationSettings.markers.show=False
188
189SC.visualizationSettings.nodes.show=True
190SC.visualizationSettings.nodes.drawNodesAsPoint = False
191SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
192SC.visualizationSettings.nodes.tiling = 4
193
194SC.visualizationSettings.window.renderWindowSize=[1200,1200]
195#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
196SC.visualizationSettings.openGL.multiSampling = 4
197#improved OpenGL rendering
198
199SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
200SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
201if False:
202 simulationSettings.solutionSettings.recordImagesInterval = 0.025
203 SC.visualizationSettings.general.graphicsUpdateInterval=2
204
205
206simulate=False
207if simulate:
208 if useGraphics:
209 SC.visualizationSettings.general.autoFitScene = False
210 exu.StartRenderer()
211 if 'renderState' in exu.sys:
212 SC.SetRenderState(exu.sys['renderState'])
213 mbs.WaitForUserToContinue()
214
215 #initial gContact statistics
216 #simulationSettings.timeIntegration.numberOfSteps = 1
217 #simulationSettings.timeIntegration.endTime = h
218 #mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
219 #print(gContact)
220
221 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
222 simulationSettings.timeIntegration.endTime = tEnd
223 simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
224 mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
225 print(gContact)
226 #p = mbs.GetNodeOutput(n, variableType=exu.OutputVariableType.Position)
227 #print("pEnd =", p[0], p[1])
228 print(gContact)
229
230 if useGraphics:
231 SC.WaitForRenderEngineStopFlag()
232 exu.StopRenderer() #safely close rendering window!
233else:
234 SC.visualizationSettings.general.autoFitScene = False
235 SC.visualizationSettings.general.graphicsUpdateInterval=0.5
236
237 print('load solution file')
238 sol = LoadSolutionFile('particles3DX.txt', safeMode=True)
239 #sol = LoadSolutionFile('coordinatesSolution2.txt')
240 print('start SolutionViewer')
241 mbs.SolutionViewer(sol)