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
 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 = graphics.Brick(p0,[LL,a,wWall],graphics.color.steelblue)
 63gFloorAdd = graphics.Brick(p0+[-0.5*LL,0.5*hWall,0],[a,hWall,wWall],color4wall)
 64gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 65gFloorAdd = graphics.Brick(p0+[ 0.5*LL,0.5*hWall,0],[a,hWall,wWall],color4wall)
 66gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 67gFloorAdd = graphics.Brick(p0+[ 0,0.5*hWall,-0.5*wWall],[LL,hWall,a],color4wall)
 68gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 69gFloorAdd = graphics.Brick(p0+[ 0,0.5*hWall, 0.5*wWall],[LL,hWall,a],color4wall)
 70gFloor = graphics.MergeTriangleLists(gFloor, gFloorAdd)
 71
 72[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(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':graphics.color.grey}]
 87    gDataList += [graphics.Sphere(point=pos, radius=rb, color= graphics.color.grey, nTiles=40)]
 88    #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,10*H,color=graphics.color.red)]
 89    nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos))
 90    #oMass = mbs.AddObject(MassPoint(mass=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 = graphics.colorList[min(9, int((int(i/rows)%rows2 * 10)/rows2))]
106    if lastColor != color4node:
107        lastColor = color4node
108        gList = [graphics.Sphere(point=[0,0,0], radius=radius, color= color4node, nTiles=8)]
109        if useClusters:
110            gList = []
111            gList += [graphics.Sphere(point=[   0,0,-drz], radius=radius, color= color4node, nTiles=8)]
112            gList += [graphics.Sphere(point=[   0,0, drz], radius=radius, color= color4node, nTiles=8)]
113            gList += [graphics.Sphere(point=[-drx,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.Brick(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=graphics.Brick(centerPoint=[0,0,0], size=[2*radius+cubeX,2*radius+cubeY,2*radius+cubeZ], color= graphics.color.steelblue)
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 += [graphics.Sphere(point=[ cubeX*(ix-0.5), cubeY*(iy-0.5), cubeZ*(iz-0.5) ], radius=radius, color= color4node, nTiles=8)]
148
149    dictMass = mbs.CreateRigidBody(
150                  inertia=RBinertia, 
151                  nodeType=exu.NodeType.RotationRotationVector,
152                  referencePosition=pRef, 
153                  initialVelocity=v0,
154                  referenceRotationMatrix=rot0,
155                  graphicsDataList=gList,
156                  returnDict=True)
157    [nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
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.solution.file.write = True
235simulationSettings.solution.file.name = 'solution/test.txt'
236
237simulationSettings.show.computationTime = True
238#simulationSettings.show.statistics = True
239simulationSettings.timeIntegration.verboseMode = 1
240simulationSettings.parallel.numberOfThreads = 8
241
242#SPEEDUPs:
243simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
244simulationSettings.timeIntegration.explicit.computeEndOfStepAccelerations = False
245simulationSettings.solution.file.writePeriod = 0.02
246simulationSettings.solution.precision = 5 #make files smaller
247simulationSettings.solution.file.export.accelerations = False
248simulationSettings.solution.file.export.velocities = False
249
250SC.visualizationSettings.general.graphicsUpdateInterval=1
251SC.visualizationSettings.general.circleTiling=200
252SC.visualizationSettings.view0.scene.drawCoordinateSystem=True
253SC.visualizationSettings.view0.scene.drawCoordinateSystem=False
254SC.visualizationSettings.loads.show=False
255SC.visualizationSettings.nodes.showBasis = False
256SC.visualizationSettings.nodes.basisSize = radius*2
257
258SC.visualizationSettings.view0.window.renderWindowSize=[1600,1200]
259SC.visualizationSettings.openGL.multiSampling = 4
260#improved OpenGL rendering
261
262SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
263if False:
264    simulationSettings.solution.recordImagesInterval = 0.025
265    SC.visualizationSettings.general.graphicsUpdateInterval=2
266
267
268simulate=True
269if simulate:
270    useGraphics = True
271    if useGraphics:
272        SC.renderer.Start()
273        SC.renderer.DoIdleTasks()
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.renderer.DoIdleTasks()
286        SC.renderer.Stop() #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)