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