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