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)