particlesTest3D2.py

You can view and download this file on Github: particlesTest3D2.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 = 50000
 35a = 0.2*L*0.5*10*0.5
 36radius = 0.5*a
 37m = 0.05
 38k = 8e4*2 #4e3 needs h=1e-4
 39d = 0.001*k*4*0.5*0.2
 40markerList = []
 41radiusList = []
 42gDataList = []
 43
 44
 45rb = 30*L
 46H = 8*L
 47Hy=3*L
 48pos0 = [0,-rb-0.5*H,0]
 49pos1 = [-rb-H,-Hy,0]
 50pos2 = [ rb+H,-Hy,0]
 51pos3 = [ 0,-Hy,rb+H]
 52pos4 = [ 0,-Hy,-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,1.2*H,14*H,color=graphics.color.red)]#80000 particles
 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
 69color4node = graphics.color.blue
 70print("start create: number of masses =",n)
 71for i in range(n):
 72
 73    kk = int(i/int(n/16))
 74    color4node = graphics.colorList[min(kk%9,9)]
 75
 76    if (i%20000 == 0): print("create mass",i)
 77    offy = 0
 78    row = 8*2 #160
 79
 80    iy = int(i/(row*row))
 81    ix = i%row
 82    iz = int(i/row)%row
 83
 84    if iy % 2 == 1:
 85        ix+=0.5
 86        iz+=0.5
 87
 88    offy = -0.25*H-3.5*a+iy*a*0.74 #0.70x is limit value!
 89    offx = -0.6*a-H*0.5 + (ix+1)*a
 90    offz = -0.6*a-H*0.5 + (iz+1)*a
 91
 92    valueRand = np.random.random(1)[0]
 93    rFact = 0.2 #random part
 94    gRad = radius*(1-rFact+rFact*valueRand)
 95    nMass = mbs.AddNode(NodePoint(referenceCoordinates=[offx,offy,offz],
 96                                  initialVelocities=[0,-20,0],
 97                                  visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
 98
 99    oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
100                                    #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
101                                    # visualization=VMassPoint(graphicsData=gData)
102                                    ))
103    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
104    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
105    markerList += [mThis]
106    radiusList += [gRad]
107
108    mLast = mThis
109print("finish create")
110#put here, such that it is transparent in background
111oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
112                                   visualization=VObjectGround(graphicsData=gDataList)))
113
114if True:
115    gContact = mbs.AddGeneralContact()
116    gContact.verboseMode = 1
117
118    for i in range(len(markerList)):
119        m = markerList[i]
120        r = radiusList[i]
121        gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
122
123    # f=n/32000
124    ssx = 20 #search tree size
125    #ssy = int(500*f) #search tree size
126    ssy = 200
127    # mbs.Assemble()
128    # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,ssx]), frictionPairingsInit=np.eye(1),
129    #                          searchTreeBoxMin=np.array([-1.2*H,-H,-1.2*H]), searchTreeBoxMax=np.array([1.2*H,14*H,1.2*H]) #80000 particles
130    #                          )
131    gContact.SetFrictionPairings(np.eye(1))
132    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssx])
133    gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-H,-1.2*H]), pMax=np.array([1.2*H,14*H,1.2*H]))
134    print('treesize=',ssx*ssx*ssy)
135
136mbs.Assemble()
137print("finish gContact")
138
139tEnd = 10
140h= 0.0001*0.25
141simulationSettings = exu.SimulationSettings()
142simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
143#simulationSettings.solutionSettings.writeSolutionToFile = True
144simulationSettings.solutionSettings.writeSolutionToFile = True
145simulationSettings.solutionSettings.solutionWritePeriod = 0.02
146simulationSettings.solutionSettings.outputPrecision = 5 #make files smaller
147simulationSettings.solutionSettings.exportAccelerations = False
148simulationSettings.solutionSettings.exportVelocities = False
149simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/test.txt'
150simulationSettings.displayComputationTime = True
151#simulationSettings.displayStatistics = True
152simulationSettings.timeIntegration.verboseMode = 1
153simulationSettings.parallel.numberOfThreads = 4
154
155simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
156simulationSettings.timeIntegration.newton.useModifiedNewton = False
157
158SC.visualizationSettings.general.graphicsUpdateInterval=0.5
159SC.visualizationSettings.general.circleTiling=200
160SC.visualizationSettings.general.drawCoordinateSystem=False
161SC.visualizationSettings.loads.show=False
162SC.visualizationSettings.bodies.show=True
163SC.visualizationSettings.markers.show=False
164
165SC.visualizationSettings.nodes.show=True
166SC.visualizationSettings.nodes.drawNodesAsPoint = False
167SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
168SC.visualizationSettings.nodes.tiling = 4
169
170SC.visualizationSettings.window.renderWindowSize=[1200,1200]
171#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
172SC.visualizationSettings.openGL.multiSampling = 4
173#improved OpenGL rendering
174
175SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
176SC.visualizationSettings.exportImages.saveImageTimeOut=10000 #5000 is too shot sometimes!
177if False:
178    simulationSettings.solutionSettings.recordImagesInterval = 0.025
179    SC.visualizationSettings.general.graphicsUpdateInterval=2
180
181
182simulate=True
183if simulate:
184    if useGraphics:
185        SC.visualizationSettings.general.autoFitScene = False
186        exu.StartRenderer()
187        if 'renderState' in exu.sys:
188            SC.SetRenderState(exu.sys['renderState'])
189        mbs.WaitForUserToContinue()
190
191    #initial gContact statistics
192    #simulationSettings.timeIntegration.numberOfSteps = 1
193    #simulationSettings.timeIntegration.endTime = h
194    #mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
195    #print(gContact)
196
197    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
198    simulationSettings.timeIntegration.endTime = tEnd
199    simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
200    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
201    #print(gContact)
202    #p = mbs.GetNodeOutput(n, variableType=exu.OutputVariableType.Position)
203    #print("pEnd =", p[0], p[1])
204    #print(gContact)
205
206    if useGraphics:
207        SC.WaitForRenderEngineStopFlag()
208        exu.StopRenderer() #safely close rendering window!
209else:
210    SC.visualizationSettings.general.autoFitScene = False
211    SC.visualizationSettings.general.graphicsUpdateInterval=0.5
212
213    print('load solution file')
214    sol = LoadSolutionFile('particles3Db.txt', safeMode=True)
215    #sol = LoadSolutionFile('coordinatesSolution2.txt')
216    print('start SolutionViewer')
217    mbs.SolutionViewer(sol)