particlesTest3D.py

You can view and download this file on Github: particlesTest3D.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  test with parallel computation and particles
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2021-11-01
  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
 20
 21SC = exu.SystemContainer()
 22mbs = SC.AddSystem()
 23
 24#create an environment for mini example
 25
 26nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 27#mLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nGround))
 28
 29np.random.seed(1) #always get same results
 30
 31useGraphics = True
 32
 33L = 1
 34n = 8000 #up to 256000; *8*4 #32*8*8
 35# n = 5000
 36a = 0.2*L*0.5*10*0.5
 37radius = 0.35*a
 38m = 0.05
 39k = 4e3*10*0.5 #4e3 needs h=1e-4
 40d = 0.001*k*4*0.5
 41markerList = []
 42radiusList = []
 43gDataList = []
 44
 45
 46rb = 30*L
 47H = 8*L
 48pos0 = [0,-rb-0.5*H,0]
 49pos1 = [-rb-H,0,0]
 50pos2 = [ rb+H,0,0]
 51pos3 = [ 0,0,rb+H]
 52pos4 = [ 0,0,-rb-H]
 53posList=[pos0,pos1,pos2,pos3,pos4]
 54for pos in posList:
 55    #gDataList += [{'type':'Circle','position':pos,'radius':rb, 'color':graphics.color.grey}]
 56    #gDataList += [graphics.Cylinder(pAxis=pos, vAxis=[0,0,0.1], radius=rb, color= graphics.color.grey, nTiles=200)]
 57    colBG = graphics.color.grey
 58    colBG[3] = 0.05
 59    gDataList += [graphics.Sphere(point=pos, radius=rb, color= colBG, nTiles=100)]
 60    #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,16*H,color=graphics.color.red)]
 61    nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos,
 62                        visualization=VNodePointGround(show=False)))
 63    #oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
 64    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 65    markerList += [mThis]
 66    radiusList += [rb]
 67
 68
 69ns = 20
 70gDataSphere = []
 71for i in range(ns):
 72    gRad = radius*(0.75+0.4*(i/ns))
 73    # gSphere = graphics.Cylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=graphics.color.blue, nTiles=12)
 74    # gSphere2 = graphics.Cylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=graphics.color.steelblue, nTiles=10)
 75    gSphere = graphics.Sphere(point=[0,0,0], radius=gRad, color=graphics.color.blue, nTiles=8)
 76    gDataSphere += [[gSphere]]
 77
 78gDataSphere = []
 79
 80color4node = graphics.color.blue
 81print("start create: number of masses =",n)
 82for i in range(n):
 83
 84    kk = int(i/12800)
 85    color4node = graphics.colorList[min(kk%12,11)]
 86    # if (i%10000 == 0):
 87        # gDataSphere = []
 88        # for i in range(ns):
 89        #     gRad = radius*(0.75+0.4*(i/ns))
 90        #     # gSphere = graphics.Cylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=graphics.color.blue, nTiles=12)
 91        #     # gSphere2 = graphics.Cylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=graphics.color.steelblue, nTiles=10)
 92        #     gSphere = graphics.Sphere(point=[0,0,0], radius=gRad, color=graphics.colorList[min(k%12,11)], nTiles=8)
 93        #     gDataSphere += [[gSphere]]
 94
 95
 96    if (i%20000 == 0): print("create mass",i)
 97    offy = 0
 98    row = 8*2 #160
 99    offy = -0.25*H-1.5*a+int(i/(row*row))*a+a*0.2*np.random.random(1)[0]
100
101    offx = -0.6*a-H*0.5 + (i%row+1)*a+0.2*a*np.random.random(1)[0]
102    offz = -0.6*a-H*0.5 + (int(i/row)%row+1)*a+0.2*a*np.random.random(1)[0]
103
104    valueRand = np.random.random(1)[0]
105    gRad = radius*(0.75+0.4*valueRand)
106    #gSphere = graphics.Cylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.25], radius=gRad, color= graphics.color.steelblue, nTiles=16)
107    #gSphere2 = graphics.Cylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.3], radius=0.8*gRad, color= graphics.color.blue, nTiles=12)
108    nMass = mbs.AddNode(NodePoint(referenceCoordinates=[offx,offy,offz],
109                                  initialVelocities=[0,-20,0],
110                                  visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
111    # gData = gDataSphere[int(valueRand*ns)]
112    # if not useGraphics:
113    #     gData = []
114    # if i%2 != 0:
115    #     gData = []
116
117    oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
118                                    #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
119                                    # visualization=VMassPoint(graphicsData=gData)
120                                    ))
121    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
122    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
123    markerList += [mThis]
124    radiusList += [gRad]
125    #if (i==n-1):
126    #    mbs.AddLoad(Force(markerNumber = mThis, loadVector = [5, -20, 0]))
127
128    #mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLast, mThis],
129    #                                    stiffness = [k,k,k], damping=[d,d,d], offset=[a,0,0],
130    #                                    visualization = VCartesianSpringDamper(drawSize = 0.1*a)))
131
132    mLast = mThis
133print("finish create")
134#put here, such that it is transparent in background
135oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
136                                   visualization=VObjectGround(graphicsData=gDataList)))
137
138if True:
139    gContact = mbs.AddGeneralContact()
140    gContact.verboseMode = 1
141
142    for i in range(len(markerList)):
143        m = markerList[i]
144        r = radiusList[i]
145        gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
146
147    # f=n/32000
148    ssx = 20 #search tree size
149    #ssy = int(500*f) #search tree size
150    ssy = 200
151    # mbs.Assemble()
152    # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,ssx]), frictionPairingsInit=np.eye(1),
153    #                          searchTreeBoxMin=np.array([-1.2*H,-0.75*H,-1.2*H]),
154    #                          searchTreeBoxMax=np.array([1.2*H,4*16*H,1.2*H])
155    #                          )
156    gContact.SetFrictionPairings(np.eye(1))
157    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssx])
158    gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-0.75*H,-1.2*H]), pMax=np.array([1.2*H,4*16*H,1.2*H]))
159    print('treesize=',ssx*ssx*ssy)
160
161mbs.Assemble()
162print("finish gContact")
163
164tEnd = 10
165h= 0.0001*0.25
166simulationSettings = exu.SimulationSettings()
167simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
168#simulationSettings.solutionSettings.writeSolutionToFile = True
169simulationSettings.solutionSettings.writeSolutionToFile = True
170simulationSettings.solutionSettings.solutionWritePeriod = 0.02
171simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
172simulationSettings.solutionSettings.exportAccelerations = False
173simulationSettings.solutionSettings.exportVelocities = False
174#simulationSettings.solutionSettings.coordinatesSolutionFileName = 'particles3D.txt'
175simulationSettings.displayComputationTime = True
176#simulationSettings.displayStatistics = True
177simulationSettings.timeIntegration.verboseMode = 1
178simulationSettings.parallel.numberOfThreads = 4
179
180simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
181simulationSettings.timeIntegration.newton.useModifiedNewton = False
182
183SC.visualizationSettings.general.graphicsUpdateInterval=0.5
184SC.visualizationSettings.general.circleTiling=200
185SC.visualizationSettings.general.drawCoordinateSystem=False
186SC.visualizationSettings.loads.show=False
187SC.visualizationSettings.bodies.show=False
188SC.visualizationSettings.markers.show=False
189
190SC.visualizationSettings.nodes.show=True
191SC.visualizationSettings.nodes.drawNodesAsPoint = False
192SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
193SC.visualizationSettings.nodes.tiling = 4
194
195SC.visualizationSettings.window.renderWindowSize=[1200,1200]
196#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
197SC.visualizationSettings.openGL.multiSampling = 4
198#improved OpenGL rendering
199
200SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
201SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
202if False:
203    simulationSettings.solutionSettings.recordImagesInterval = 0.025
204    SC.visualizationSettings.general.graphicsUpdateInterval=2
205
206
207simulate=False
208if simulate:
209    if useGraphics:
210        SC.visualizationSettings.general.autoFitScene = False
211        exu.StartRenderer()
212        if 'renderState' in exu.sys:
213            SC.SetRenderState(exu.sys['renderState'])
214        mbs.WaitForUserToContinue()
215
216    #initial gContact statistics
217    #simulationSettings.timeIntegration.numberOfSteps = 1
218    #simulationSettings.timeIntegration.endTime = h
219    #mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
220    #print(gContact)
221
222    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
223    simulationSettings.timeIntegration.endTime = tEnd
224    simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
225    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
226    print(gContact)
227    #p = mbs.GetNodeOutput(n, variableType=exu.OutputVariableType.Position)
228    #print("pEnd =", p[0], p[1])
229    print(gContact)
230
231    if useGraphics:
232        SC.WaitForRenderEngineStopFlag()
233        exu.StopRenderer() #safely close rendering window!
234else:
235    SC.visualizationSettings.general.autoFitScene = False
236    SC.visualizationSettings.general.graphicsUpdateInterval=0.5
237
238    #load previously computed solution
239    # print('load solution file')
240    # sol = LoadSolutionFile('particles3DX.txt', safeMode=True)
241    print('start SolutionViewer')
242    # mbs.SolutionViewer(sol)
243    mbs.SolutionViewer()