particlesTest.py

You can view and download this file on Github: particlesTest.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
 30L = 1
 31n = 32000 #*8*4 #32*8*8
 32a = 0.2*L*0.5
 33radius = 0.35*a
 34m = 0.05
 35k = 4e3*10 #4e3 needs h=1e-4
 36d = 0.001*k*4*0.5
 37markerList = []
 38radiusList = []
 39gDataList = []
 40
 41
 42rb = 30*L
 43H = 8*L
 44pos0 = [0,-rb-0.5*H,0]
 45pos1 = [-rb-H,0,0]
 46pos2 = [ rb+H,0,0]
 47posList=[pos0,pos1,pos2]
 48for pos in posList:
 49    #gDataList += [{'type':'Circle','position':pos,'radius':rb, 'color':color4grey}]
 50    gDataList += [GraphicsDataCylinder(pAxis=pos, vAxis=[0,0,0.1], radius=rb, color= color4grey, nTiles=200)]
 51    #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,10*H,color=color4red)]
 52    nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos))
 53    #oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
 54    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 55    markerList += [mThis]
 56    radiusList += [rb]
 57
 58    oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 59                                       visualization=VObjectGround(graphicsData=gDataList)))
 60
 61ns = 20
 62gDataSphere = []
 63for i in range(ns):
 64    gRad = radius*(0.75+0.4*(i/ns))
 65    gSphere = GraphicsDataCylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=color4blue, nTiles=12)
 66    gSphere2 = GraphicsDataCylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=color4steelblue, nTiles=10)
 67    gDataSphere += [[gSphere, gSphere2]]
 68
 69
 70print("start create: number of masses =",n)
 71for i in range(n):
 72    if (i%20000 == 0): print("create mass",i)
 73    offy = 0
 74    row = 80*2
 75    offy = -2*L-a*9+int(i/row)*a+a*0.2*np.random.random(1)[0]
 76
 77    valueRand = np.random.random(1)[0]
 78    gRad = radius*(0.75+0.4*valueRand)
 79    #gSphere = GraphicsDataCylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.25], radius=gRad, color= color4steelblue, nTiles=16)
 80    #gSphere2 = GraphicsDataCylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.3], radius=0.8*gRad, color= color4blue, nTiles=12)
 81    nMass = mbs.AddNode(NodePoint(referenceCoordinates=[-0.6*a-H + (i%row+1)*a+0.2*a*np.random.random(1)[0],offy,0],
 82                                  initialVelocities=[0,-5,0]))
 83    oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
 84                                    #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
 85                                    visualization=VMassPoint(graphicsData=gDataSphere[int(valueRand*ns)])
 86                                    ))
 87    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 88    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
 89    markerList += [mThis]
 90    radiusList += [gRad]
 91    #if (i==n-1):
 92    #    mbs.AddLoad(Force(markerNumber = mThis, loadVector = [5, -20, 0]))
 93
 94    #mbs.AddObject(CartesianSpringDamper(markerNumbers=[mLast, mThis],
 95    #                                    stiffness = [k,k,k], damping=[d,d,d], offset=[a,0,0],
 96    #                                    visualization = VCartesianSpringDamper(drawSize = 0.1*a)))
 97
 98    mLast = mThis
 99print("finish create")
100
101if True:
102    gContact = mbs.AddGeneralContact()
103    gContact.verboseMode = 1
104
105    for i in range(len(markerList)):
106        m = markerList[i]
107        r = radiusList[i]
108        gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
109
110    f=n/32000
111    ssx = 100 #search tree size
112    ssy = int(500*f) #search tree size
113    # mbs.Assemble()
114    # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,1]), frictionPairingsInit=np.eye(1),
115    #                          searchTreeBoxMin=np.array([-1.2*H,-0.75*H,0]), searchTreeBoxMax=np.array([1.2*H,10*H*f,1])
116    #                          )
117    gContact.SetFrictionPairings(np.eye(1))
118    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,1])
119    gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-0.75*H,0]), pMax=np.array([1.2*H,10*H*f,1]))
120
121mbs.Assemble()
122print("finish gContact")
123
124tEnd = 20
125h= 0.0001*0.25
126simulationSettings = exu.SimulationSettings()
127simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
128#simulationSettings.solutionSettings.writeSolutionToFile = True
129simulationSettings.solutionSettings.writeSolutionToFile = True
130simulationSettings.solutionSettings.solutionWritePeriod = 0.02
131
132simulationSettings.displayComputationTime = True
133#simulationSettings.displayStatistics = True
134simulationSettings.timeIntegration.verboseMode = 1
135simulationSettings.parallel.numberOfThreads = 8
136
137simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
138simulationSettings.timeIntegration.newton.useModifiedNewton = False
139
140SC.visualizationSettings.general.graphicsUpdateInterval=10
141SC.visualizationSettings.general.circleTiling=200
142SC.visualizationSettings.general.drawCoordinateSystem=False
143SC.visualizationSettings.loads.show=False
144SC.visualizationSettings.window.renderWindowSize=[1600,1200]
145#SC.visualizationSettings.window.renderWindowSize=[1024,1400]
146SC.visualizationSettings.openGL.multiSampling = 4
147#improved OpenGL rendering
148
149SC.visualizationSettings.exportImages.saveImageFileName = "animation/frame"
150if False:
151    simulationSettings.solutionSettings.recordImagesInterval = 0.025
152    SC.visualizationSettings.general.graphicsUpdateInterval=2
153
154
155simulate=True
156if simulate:
157    useGraphics = True
158    if useGraphics:
159        exu.StartRenderer()
160        mbs.WaitForUserToContinue()
161
162    #initial gContact statistics
163    #simulationSettings.timeIntegration.numberOfSteps = 1
164    #simulationSettings.timeIntegration.endTime = h
165    #mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
166    #print(gContact)
167
168    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
169    simulationSettings.timeIntegration.endTime = tEnd
170    simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
171    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.RK67)
172    print(gContact)
173    #p = mbs.GetNodeOutput(n, variableType=exu.OutputVariableType.Position)
174    #print("pEnd =", p[0], p[1])
175    print(gContact)
176
177    if useGraphics:
178        SC.WaitForRenderEngineStopFlag()
179        exu.StopRenderer() #safely close rendering window!
180else:
181    SC.visualizationSettings.general.autoFitScene = False
182    SC.visualizationSettings.general.graphicsUpdateInterval=0.5
183
184    sol = LoadSolutionFile('particles.txt')
185    mbs.SolutionViewer(sol)