particlesSilo.py
You can view and download this file on Github: particlesSilo.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 * #includes itemInterface and rigidBodyUtilities
15import exudyn.graphics as graphics #only import if it does not conflict
16from exudyn.graphicsDataUtilities import *
17
18import numpy as np
19
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23
24#create an environment for mini example
25
26nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
27#mLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nGround))
28
29np.random.seed(1) #always get same results
30
31useGraphics = True
32
33
34useRigidBody = True
35L = 1
36n = 4000*25 #test: 4000
37n = 8000 #fast simulation for testing
38row = 8*2
39a = L*0.5*0.5
40h= 0.0001
41m = 0.05
42ss=16*2
43holeRad = 3*a
44
45if n >= 4000*8:
46 a*=0.4
47 row=40
48 ss = 16*3
49 holeRad *= 1.4 #better 1.4 !
50
51if n >= 4000*64:
52 a*=0.5
53 row*=2
54 h *= 0.5
55 m *=0.125
56 ss*=2
57
58radius = 0.5*a
59t = 0.5*a
60k = 8e4*2 #4e3 needs h=1e-4
61d = 0.001*k*4*0.5*0.2
62
63frictionCoeff = 0.5
64if not useRigidBody:
65 frictionCoeff = 0
66
67markerList = []
68radiusList = []
69gDataList = []
70
71# rb = 30*L
72H = 8*L
73Hy=3*L
74
75gContact = mbs.AddGeneralContact()
76gContact.verboseMode = 1
77gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
78gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
79#gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-H,-1.2*H]), pMax=np.array([1.2*H,14*H,1.2*H]))
80#print('treesize=',ssx*ssx*ssy)
81
82
83#%% ground
84LL=6*L
85p0 = np.array([0,0,-0.5*t])
86color4wall = [0.6,0.6,0.6,0.5]
87addNormals = False
88hw=10*a
89gFloor = graphics.Brick(p0,[LL,LL,t],graphics.color.steelblue,addNormals)
90gFloorAdd = graphics.Brick(p0+[-0.5*LL,0,0.5*hw],[t,LL,hw],color4wall,addNormals)
91gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
92gFloorAdd = graphics.Brick(p0+[ 0.5*LL,0,0.5*hw],[t,LL,hw],color4wall,addNormals)
93gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
94gFloorAdd = graphics.Brick(p0+[0,-0.5*LL,0.5*hw],[LL,t,hw],color4wall,addNormals)
95gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
96gFloorAdd = graphics.Brick(p0+[0, 0.5*LL,0.5*hw],[LL,t,hw],color4wall,addNormals)
97gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
98
99gDataList = [gFloor]
100
101
102nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
103mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
104#mGroundC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
105
106[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
107#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
108# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
109gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
110 pointList=meshPoints, triangleList=meshTrigs)
111
112if True: #looses color
113 gFloor = graphics.FromPointsAndTrigs(meshPoints, meshTrigs, color=color4wall) #show refined mesh
114 gDataList = [gFloor]
115
116
117color4node = graphics.color.blue
118print("start create: number of masses =",n)
119for i in range(n):
120
121 kk = int(i/int(n/8))
122 color4node = graphics.colorList[min(kk%9,9)]
123
124 if (i%20000 == 0): print("create mass",i)
125 offy = 0
126
127 iz = int(i/(row*row))
128 ix = i%row
129 iy = int(i/row)%row
130
131 if iz % 2 == 1:
132 ix+=0.5
133 iy+=0.5
134
135 offz = 5*L+0.5*a+iz*a*0.74 #0.70x is limit value!
136 offx = -0.6*a-row*0.5*a + (ix+1)*a
137 offy = -0.6*a-row*0.5*a + (iy+1)*a
138
139 valueRand = np.random.random(1)[0]
140 rFact = 0.2 #random part
141 gRad = radius*(1-rFact+rFact*valueRand)
142 v0 = [0,0,-2]
143 pRef = [offx,offy,offz]
144
145 if not useRigidBody:
146 nMass = mbs.AddNode(NodePoint(referenceCoordinates=pRef,
147 initialVelocities=v0,
148 visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
149
150 oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
151 #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
152 # visualization=VMassPoint(graphicsData=gData)
153 ))
154 mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
155 else:
156 RBinertia = InertiaSphere(m, radius)
157 [nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
158 nodeType=exu.NodeType.RotationRotationVector,
159 position=pRef, velocity=v0,
160 #graphicsDataList=gList,
161 )
162 mbs.SetNodeParameter(nMass, 'VdrawSize', 2*gRad)
163 mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
164 mThis = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
165
166 mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,0,-m*9.81]))
167
168 gContact.AddSphereWithMarker(mThis, radius=gRad, contactStiffness=k, contactDamping=d,
169 frictionMaterialIndex=0)
170
171
172
173if True: #add Silo
174 SR = 3.1*L
175 SH = 2*L
176 SH2 = 1*L #hole
177 SR2 = holeRad #hole
178 ST = 0.25*L
179 #contour=8*np.array([[0,0.2],[0.3,0.2],[0.5,0.3],[0.7,0.4],[1,0.4],[1,0.]])
180 contour=np.array([[0,SR2],[0,SR2+ST],[SH2-ST,SR2+ST],[2*SH2-ST,SR+ST],[2*SH2+SH,SR+ST],
181 [2*SH2+SH,SR],[2*SH2,SR],[SH2,SR2],[0,SR2]])
182 contour = list(contour)
183 contour.reverse()
184 gSilo = graphics.SolidOfRevolution(pAxis=[0,0,3*L], vAxis=[0,0,1],
185 contour=contour, color=[0.8,0.1,0.1,0.5], nTiles = 64)
186
187 [meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gSilo)
188 gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
189 pointList=meshPoints, triangleList=meshTrigs)
190
191
192#put here, such that it is transparent in background
193oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
194 visualization=VObjectGround(graphicsData=[gSilo]+gDataList)))
195
196
197mbs.Assemble()
198print("finish gContact")
199
200items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
201print('n spheres=',len(items['MarkerBasedSpheres']))
202
203
204tEnd = 50
205#tEnd = h*100
206simulationSettings = exu.SimulationSettings()
207simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
208#simulationSettings.solutionSettings.writeSolutionToFile = True
209simulationSettings.solutionSettings.writeSolutionToFile = True
210simulationSettings.solutionSettings.solutionWritePeriod = 0.01
211simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
212simulationSettings.solutionSettings.exportAccelerations = False
213simulationSettings.solutionSettings.exportVelocities = False
214simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
215simulationSettings.displayComputationTime = True
216#simulationSettings.displayStatistics = True
217simulationSettings.timeIntegration.verboseMode = 1
218simulationSettings.parallel.numberOfThreads = 8
219
220simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
221simulationSettings.timeIntegration.newton.useModifiedNewton = False
222
223SC.visualizationSettings.general.graphicsUpdateInterval=0.5*4
224SC.visualizationSettings.general.circleTiling=200
225SC.visualizationSettings.general.drawCoordinateSystem=True
226SC.visualizationSettings.loads.show=False
227SC.visualizationSettings.bodies.show=True
228SC.visualizationSettings.markers.show=False
229
230SC.visualizationSettings.nodes.show=True
231SC.visualizationSettings.nodes.drawNodesAsPoint = False
232SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
233SC.visualizationSettings.nodes.tiling = 4
234
235SC.visualizationSettings.window.renderWindowSize=[1200,1200]
236#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
237SC.visualizationSettings.openGL.multiSampling = 4
238#improved OpenGL rendering
239
240SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
241SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
242
243if False:
244 simulationSettings.solutionSettings.recordImagesInterval = 0.005
245 SC.visualizationSettings.general.graphicsUpdateInterval=2
246
247
248simulate=True
249if simulate:
250 if useGraphics:
251 SC.visualizationSettings.general.autoFitScene = False
252 exu.StartRenderer()
253 if 'renderState' in exu.sys:
254 SC.SetRenderState(exu.sys['renderState'])
255 mbs.WaitForUserToContinue()
256
257 #initial gContact statistics
258 #simulationSettings.timeIntegration.numberOfSteps = 1
259 #simulationSettings.timeIntegration.endTime = h
260 #mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
261 #print(gContact)
262
263 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
264 simulationSettings.timeIntegration.endTime = tEnd
265 mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
266 #print(gContact)
267 #p = mbs.GetNodeOutput(n, variableType=exu.OutputVariableType.Position)
268 #print("pEnd =", p[0], p[1])
269 #print(gContact)
270
271 if useGraphics:
272 SC.WaitForRenderEngineStopFlag()
273 exu.StopRenderer() #safely close rendering window!
274
275if not simulate:
276 SC.visualizationSettings.general.autoFitScene = False
277 SC.visualizationSettings.general.graphicsUpdateInterval=0.5
278
279 print('load solution file')
280 #sol = LoadSolutionFile('solution/test2.txt', safeMode=False)
281 sol = LoadSolutionFile('solution/test.txt', safeMode=True, verbose = True)#, maxRows=100)
282 print('start SolutionViewer')
283 mbs.SolutionViewer(sol)