particleClusters.py
You can view and download this file on Github: particleClusters.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: test with clusters of spheres attached to rigid body
5#
6# Author: Johannes Gerstmayr
7# Date: 2021-12-23
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
19from math import sqrt
20
21SC = exu.SystemContainer()
22mbs = SC.AddSystem()
23
24np.random.seed(1) #always get same results
25
26L = 1
27# n = 16000 #*8*4 #32*8*8
28a = 0.2*L
29# radius = 0.5*0.125
30# rows = 14
31# rows2 = 40*2
32# mu = 0.1
33m = 0.0125
34k = 4e3 #2e4
35d = 0.0004*k
36# ssx = 32 #search tree size
37# ssy = 16 #search tree size
38# ssz = 12 #search tree size
39LL=12*L
40hWall = 8*L
41wWall = 4*L
42
43useClusters = True
44useNodeMarker = not useClusters
45n = 200#0*3
46rows = 7
47rows2 = 12
48radius = 0.5*0.25
49mu = 0.4*3
50
51ssx = 16 #search tree size
52ssy = 10 #search tree size
53ssz = 8 #search tree size
54hWall = 12*L
55drx = radius*0.75 #for cluster
56drz = radius*0.5 #for cluster
57
58markerList = []
59radiusList = []
60p0 = np.array([0.45*LL,-4*radius,0])
61color4wall = [0.9,0.9,0.7,0.25]
62gFloor = GraphicsDataOrthoCubePoint(p0,[LL,a,wWall],color4steelblue)
63gFloorAdd = GraphicsDataOrthoCubePoint(p0+[-0.5*LL,0.5*hWall,0],[a,hWall,wWall],color4wall)
64gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
65gFloorAdd = GraphicsDataOrthoCubePoint(p0+[ 0.5*LL,0.5*hWall,0],[a,hWall,wWall],color4wall)
66gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
67gFloorAdd = GraphicsDataOrthoCubePoint(p0+[ 0,0.5*hWall,-0.5*wWall],[LL,hWall,a],color4wall)
68gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
69gFloorAdd = GraphicsDataOrthoCubePoint(p0+[ 0,0.5*hWall, 0.5*wWall],[LL,hWall,a],color4wall)
70gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
71
72[meshPoints, meshTrigs] = GraphicsData2PointsAndTrigs(gFloor)
73
74gDataList = [gFloor]
75
76
77rb = radius
78H = 8*L
79RR = 2*radius*sqrt(0.5)-1e-12
80pos0 = [0.1,-RR*2,0]
81pos1 = [5*radius,-RR*2,0]
82#posList=[pos0, pos1]
83#posList=[pos0]
84posList=[]
85for pos in posList:
86 #gDataList += [{'type':'Circle','position':pos,'radius':rb, 'color':color4grey}]
87 gDataList += [GraphicsDataSphere(point=pos, radius=rb, color= color4grey, nTiles=40)]
88 #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,10*H,color=color4red)]
89 nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos))
90 #oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
91 mThis = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
92 markerList += [mThis]
93 radiusList += [rb]
94
95
96
97#ns = 20
98minP = np.array([1e10,1e10,1e10])
99maxP = -minP
100gContact = mbs.AddGeneralContact()
101
102lastColor=[0,0,0,0]
103for i in range(n):
104
105 color4node = color4list[min(9, int((int(i/rows)%rows2 * 10)/rows2))]
106 if lastColor != color4node:
107 lastColor = color4node
108 gList = [GraphicsDataSphere(point=[0,0,0], radius=radius, color= color4node, nTiles=8)]
109 if useClusters:
110 gList = []
111 gList += [GraphicsDataSphere(point=[ 0,0,-drz], radius=radius, color= color4node, nTiles=8)]
112 gList += [GraphicsDataSphere(point=[ 0,0, drz], radius=radius, color= color4node, nTiles=8)]
113 gList += [GraphicsDataSphere(point=[-drx,0,-drz], radius=radius, color= color4node, nTiles=8)]
114 gList += [GraphicsDataSphere(point=[ drx,0,-drz], radius=radius, color= color4node, nTiles=8)]
115 gList += [GraphicsDataSphere(point=[ drx,0, drz], radius=radius, color= color4node, nTiles=8)]
116 gList += [GraphicsDataSphere(point=[-drx,0, drz], radius=radius, color= color4node, nTiles=8)]
117 #gList = [GraphicsDataOrthoCubePoint(centerPoint=[0,0,0], size=[2*radius+2*drx,2*radius,2*radius+2*drz], color= color4node)]
118
119
120 dd = 2.5*radius
121 ox=(int(i/rows)%2)*0.5*dd#rows
122 pRef = [(i%rows)*(dd+2*drx)+ox, (int(i/rows)%rows2)*0.8*dd, -wWall*0.5+dd+ox+int(i/(rows*rows2))*(dd+drz)]
123
124 minP = np.minimum(minP, pRef)
125 maxP = np.maximum(maxP, pRef)
126 v0 = [0,-10*0.1,0]
127 rot0 = np.eye(3)
128 forceY = -m*9.81*0.25
129
130 RBinertia = InertiaSphere(m, radius*1)
131 isCube = (i == n-1)
132 cubeX = 8*radius
133 cubeY = 8*radius
134 cubeZ = 8*radius
135 if isCube:
136 pRef = [10*L, 2*L,-wWall*0.5+cubeZ]
137 rot0=RotationMatrixZ(0.2) @ RotationMatrixX(0.1)
138 RBinertia = InertiaSphere(10*m, radius*8)
139 forceY *= 10
140 gCube=GraphicsDataOrthoCubePoint(centerPoint=[0,0,0], size=[2*radius+cubeX,2*radius+cubeY,2*radius+cubeZ], color= color4steelblue)
141 gList = [gCube]
142
143 if isCube:
144 for ix in range(2):
145 for iy in range(2):
146 for iz in range(2):
147 gList += [GraphicsDataSphere(point=[ cubeX*(ix-0.5), cubeY*(iy-0.5), cubeZ*(iz-0.5) ], radius=radius, color= color4node, nTiles=8)]
148
149 [nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
150 #nodeType=exu.NodeType.RotationRxyz,
151 nodeType=exu.NodeType.RotationRotationVector,
152 position=pRef, velocity=v0,
153 rotationMatrix=rot0,
154 #angularVelocity=omega0,
155 #gravity=[0.,-9.81,0.],
156 graphicsDataList=gList,
157 )
158 mbs.SetNodeParameter(nMass, 'VdrawSize', radius*2)
159 mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
160
161 if i == n-2:
162 nMassOutput = nMass #for solution output
163
164
165 mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
166 if useNodeMarker:
167 markerList += [mNode]
168 elif not useClusters:
169 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[0,0,0]))
170 markerList += [mBody]
171 else:
172 if isCube:
173 [meshPointsCube, meshTrigsCube] = GraphicsData2PointsAndTrigs(gCube)
174 mCube = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[0,0,0]))
175 gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCube, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
176 pointList=meshPointsCube, triangleList=meshTrigsCube)
177 for ix in range(2):
178 for iy in range(2):
179 for iz in range(2):
180 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ cubeX*(ix-0.5), cubeY*(iy-0.5), cubeZ*(iz-0.5)]))
181 markerList += [mBody]
182 else:
183 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ 0, 0,-drz]))
184 markerList += [mBody]
185 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ 0, 0, drz]))
186 markerList += [mBody]
187 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[-drx, 0,-drz]))
188 markerList += [mBody]
189 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ drx, 0,-drz]))
190 markerList += [mBody]
191 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ drx, 0, drz]))
192 markerList += [mBody]
193 mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[-drx, 0, drz]))
194 markerList += [mBody]
195
196
197 #mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
198 #if (i < 2):
199 mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,forceY,0]))
200
201oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
202 visualization=VObjectGround(graphicsData=gDataList)))
203mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
204
205gContact.verboseMode = 1
206#gContact.sphereSphereContact = False
207gContact.frictionProportionalZone = 1e-6
208#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs)
209gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
210 pointList=meshPoints, triangleList=meshTrigs)
211
212for i in range(len(markerList)):
213 m = markerList[i]
214 gContact.AddSphereWithMarker(m, radius=radius, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
215
216#gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,1]), frictionPairingsInit=np.eye(1),
217# searchTreeBoxMin=np.array([-H,-H,-H]), searchTreeBoxMax=np.array([H,H,H])
218# )
219gContact.SetFrictionPairings(mu*np.eye(1))
220gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
221#gContact.SetSearchTreeBox(pMin=np.array([-H,-H,-0.2*H]), pMax=np.array([H,H,0.2*H]))
222
223sNode = mbs.AddSensor(SensorNode(nodeNumber=nMassOutput, fileName='solution/particlesNode.txt',
224 outputVariableType=exu.OutputVariableType.Position))
225
226mbs.Assemble()
227
228tEnd = 1
229h= 0.0001
230
231
232simulationSettings = exu.SimulationSettings()
233simulationSettings.solutionSettings.writeSolutionToFile = True
234simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
235
236simulationSettings.displayComputationTime = True
237#simulationSettings.displayStatistics = True
238simulationSettings.timeIntegration.verboseMode = 1
239simulationSettings.parallel.numberOfThreads = 8
240
241#SPEEDUPs:
242simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
243simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False
244simulationSettings.solutionSettings.solutionWritePeriod = 0.02
245simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
246simulationSettings.solutionSettings.exportAccelerations = False
247simulationSettings.solutionSettings.exportVelocities = False
248
249SC.visualizationSettings.general.graphicsUpdateInterval=1
250SC.visualizationSettings.general.circleTiling=200
251SC.visualizationSettings.general.drawCoordinateSystem=True
252SC.visualizationSettings.general.drawCoordinateSystem=False
253SC.visualizationSettings.loads.show=False
254SC.visualizationSettings.nodes.showBasis = False
255SC.visualizationSettings.nodes.basisSize = radius*2
256
257SC.visualizationSettings.window.renderWindowSize=[1600,1200]
258SC.visualizationSettings.openGL.multiSampling = 4
259#improved OpenGL rendering
260
261SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
262if False:
263 simulationSettings.solutionSettings.recordImagesInterval = 0.025
264 SC.visualizationSettings.general.graphicsUpdateInterval=2
265
266
267simulate=True
268if simulate:
269 useGraphics = True
270 if useGraphics:
271 exu.StartRenderer()
272 mbs.WaitForUserToContinue()
273
274 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
275 simulationSettings.timeIntegration.endTime = tEnd
276 mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
277
278 p = mbs.GetSensorValues(sNode)
279 #after one second (200 particles, h=0.0001):
280 #p= [0.854546923638082, 0.4801722062275968, -0.7822195683611741]
281 print("p=", list(p))
282
283 if useGraphics:
284 SC.WaitForRenderEngineStopFlag()
285 exu.StopRenderer() #safely close rendering window!
286
287 if True:
288
289 mbs.PlotSensor(sensorNumbers=[sNode,sNode,sNode], components=[0,1,2], closeAll=True)
290else:
291 SC.visualizationSettings.general.autoFitScene = False
292 SC.visualizationSettings.general.graphicsUpdateInterval=0.1
293
294 sol = LoadSolutionFile('particles.txt')
295 mbs.SolutionViewer(sol)