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    dictMass = mbs.CreateRigidBody(
151                  inertia=RBinertia,
152                  nodeType=exu.NodeType.RotationRotationVector,
153                  referencePosition=pRef,
154                  initialVelocity=v0,
155                  referenceRotationMatrix=rot0,
156                  graphicsDataList=gList,
157                  returnDict=True)
158    [nMass, oMass] = [dictMass['nodeNumber'], dictMass['bodyNumber']]
159
160    mbs.SetNodeParameter(nMass, 'VdrawSize', radius*2)
161    mbs.SetNodeParameter(nMass, 'Vcolor', color4node)
162
163    if i == n-2:
164        nMassOutput = nMass #for solution output
165
166
167    mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
168    if useNodeMarker:
169        markerList += [mNode]
170    elif not useClusters:
171        mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[0,0,0]))
172        markerList += [mBody]
173    else:
174        if isCube:
175            [meshPointsCube, meshTrigsCube] = graphics.ToPointsAndTrigs(gCube)
176            mCube = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[0,0,0]))
177            gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCube, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
178                                                pointList=meshPointsCube,  triangleList=meshTrigsCube)
179            for ix in range(2):
180                for iy in range(2):
181                    for iz in range(2):
182                        mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[ cubeX*(ix-0.5), cubeY*(iy-0.5), cubeZ*(iz-0.5)]))
183                        markerList += [mBody]
184        else:
185            mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[   0, 0,-drz]))
186            markerList += [mBody]
187            mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[   0, 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            mBody = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMass, localPosition=[-drx, 0, drz]))
196            markerList += [mBody]
197
198
199    #mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
200    #if (i < 2):
201    mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,forceY,0]))
202
203oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
204                                    visualization=VObjectGround(graphicsData=gDataList)))
205mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
206
207gContact.verboseMode = 1
208#gContact.sphereSphereContact = False
209gContact.frictionProportionalZone = 1e-6
210#[meshPoints,  meshTrigs] = RefineMesh(meshPoints,  meshTrigs)
211gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
212    pointList=meshPoints,  triangleList=meshTrigs)
213
214for i in range(len(markerList)):
215    m = markerList[i]
216    gContact.AddSphereWithMarker(m, radius=radius, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
217
218#gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,1]), frictionPairingsInit=np.eye(1),
219#                         searchTreeBoxMin=np.array([-H,-H,-H]), searchTreeBoxMax=np.array([H,H,H])
220#                         )
221gContact.SetFrictionPairings(mu*np.eye(1))
222gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssz])
223#gContact.SetSearchTreeBox(pMin=np.array([-H,-H,-0.2*H]), pMax=np.array([H,H,0.2*H]))
224
225sNode = mbs.AddSensor(SensorNode(nodeNumber=nMassOutput, fileName='solution/particlesNode.txt',
226                         outputVariableType=exu.OutputVariableType.Position))
227
228mbs.Assemble()
229
230tEnd = 1
231h= 0.0001
232
233
234simulationSettings = exu.SimulationSettings()
235simulationSettings.solutionSettings.writeSolutionToFile = True
236simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
237
238simulationSettings.displayComputationTime = True
239#simulationSettings.displayStatistics = True
240simulationSettings.timeIntegration.verboseMode = 1
241simulationSettings.parallel.numberOfThreads = 8
242
243#SPEEDUPs:
244simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
245simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False
246simulationSettings.solutionSettings.solutionWritePeriod = 0.02
247simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
248simulationSettings.solutionSettings.exportAccelerations = False
249simulationSettings.solutionSettings.exportVelocities = False
250
251SC.visualizationSettings.general.graphicsUpdateInterval=1
252SC.visualizationSettings.general.circleTiling=200
253SC.visualizationSettings.general.drawCoordinateSystem=True
254SC.visualizationSettings.general.drawCoordinateSystem=False
255SC.visualizationSettings.loads.show=False
256SC.visualizationSettings.nodes.showBasis = False
257SC.visualizationSettings.nodes.basisSize = radius*2
258
259SC.visualizationSettings.window.renderWindowSize=[1600,1200]
260SC.visualizationSettings.openGL.multiSampling = 4
261#improved OpenGL rendering
262
263SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
264if False:
265    simulationSettings.solutionSettings.recordImagesInterval = 0.025
266    SC.visualizationSettings.general.graphicsUpdateInterval=2
267
268
269simulate=True
270if simulate:
271    useGraphics = True
272    if useGraphics:
273        exu.StartRenderer()
274        mbs.WaitForUserToContinue()
275
276    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
277    simulationSettings.timeIntegration.endTime = tEnd
278    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
279
280    p = mbs.GetSensorValues(sNode)
281    #after one second (200 particles, h=0.0001):
282    #p= [0.854546923638082, 0.4801722062275968, -0.7822195683611741]
283    print("p=", list(p))
284
285    if useGraphics:
286        SC.WaitForRenderEngineStopFlag()
287        exu.StopRenderer() #safely close rendering window!
288
289    if True:
290
291        mbs.PlotSensor(sensorNumbers=[sNode,sNode,sNode], components=[0,1,2], closeAll=True)
292else:
293    SC.visualizationSettings.general.autoFitScene = False
294    SC.visualizationSettings.general.graphicsUpdateInterval=0.1
295
296    sol = LoadSolutionFile('particles.txt')
297    mbs.SolutionViewer(sol)