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