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# Revised:  2024-03-24
  9#
 10# 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.
 11#
 12#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 13
 14import exudyn as exu
 15from exudyn.itemInterface import *
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17import exudyn.graphics as graphics #only import if it does not conflict
 18from exudyn.graphicsDataUtilities import *
 19
 20import numpy as np
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25#create an environment for mini example
 26
 27nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 28#mLast = mbs.AddMarker(MarkerNodePosition(nodeNumber=nGround))
 29
 30np.random.seed(1) #always get same results
 31
 32L = 1
 33n = 4000 #*8*4 #32*8*8
 34a = L*0.25
 35radius = 0.35*a
 36m = 0.05
 37k = 4e4 #4e3 needs h=1e-4
 38d = 0.00005*k
 39markerList = []
 40radiusList = []
 41gDataList = []
 42
 43
 44rb = 30*L
 45H = 8*L
 46pos0 = [0,-rb-0.5*H,0]
 47pos1 = [-rb-H,0,0]
 48pos2 = [ rb+H,0,0]
 49posList=[pos0,pos1,pos2]
 50for pos in posList:
 51    #gDataList += [{'type':'Circle','position':pos,'radius':rb, 'color':graphics.color.grey}]
 52    gDataList += [graphics.Cylinder(pAxis=pos, vAxis=[0,0,0.1], radius=rb, color= graphics.color.grey, nTiles=200)]
 53    #gDataList += [GraphicsDataRectangle(-1.2*H,-H*0.75,1.2*H,10*H,color=graphics.color.red)]
 54    nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos))
 55    #oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass))
 56    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 57    markerList += [mThis]
 58    radiusList += [rb]
 59
 60    oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
 61                                       visualization=VObjectGround(graphicsData=gDataList)))
 62
 63ns = 20
 64gDataSphere = []
 65for i in range(ns):
 66    gRad = radius*(0.75+0.4*(i/ns))
 67    gSphere = graphics.Cylinder(pAxis=[0,0,-0.25], vAxis=[0,0,0.5], radius=gRad, color=graphics.color.blue, nTiles=12)
 68    gSphere2 = graphics.Cylinder(pAxis=[0,0,-0.3], vAxis=[0,0,0.6], radius=0.8*gRad, color=graphics.color.steelblue, nTiles=10)
 69    gDataSphere += [[gSphere, gSphere2]]
 70
 71
 72print("start create: number of masses =",n)
 73for i in range(n):
 74    if (i%20000 == 0): print("create mass",i)
 75    offy = 0
 76    row = 32*2
 77    offy = -2*L-a*0+int(i/row)*a+a*0.2*np.random.random(1)[0]
 78
 79    valueRand = np.random.random(1)[0]
 80    gRad = radius*(0.75+0.4*valueRand)
 81
 82    nMass = mbs.AddNode(NodePoint(referenceCoordinates=[-0.6*a-H + (i%row+1)*a+0.2*a*np.random.random(1)[0],offy,0],
 83                                  initialVelocities=[0,-5,0]))
 84    oMass = mbs.AddObject(MassPoint(physicsMass=m, nodeNumber=nMass,
 85                                    #visualization=VMassPoint(graphicsData=[gSphere,gSphere2])
 86                                    visualization=VMassPoint(graphicsData=gDataSphere[int(valueRand*ns)])
 87                                    ))
 88    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 89    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
 90    markerList += [mThis]
 91    radiusList += [gRad]
 92
 93    mLast = mThis
 94print("finish create")
 95
 96if True:
 97    gContact = mbs.AddGeneralContact()
 98    gContact.verboseMode = 1
 99
100    for i in range(len(markerList)):
101        m = markerList[i]
102        r = radiusList[i]
103        gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d,
104                                     frictionMaterialIndex=-1)
105
106    ssx = int(sqrt(n)+1) #search tree size
107    ssy = ssx #search tree size
108
109    gContact.SetFrictionPairings(np.eye(1))
110    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,1])
111    gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-0.75*H,0]), pMax=np.array([1.2*H,10*H,1]))
112
113mbs.Assemble()
114print("finish gContact")
115
116tEnd = 20
117stepSize = 0.0005
118simulationSettings = exu.SimulationSettings()
119simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
120#simulationSettings.solutionSettings.writeSolutionToFile = True
121simulationSettings.solutionSettings.writeSolutionToFile = True
122simulationSettings.solutionSettings.solutionWritePeriod = 0.04
123
124simulationSettings.displayComputationTime = True
125#simulationSettings.displayStatistics = True
126simulationSettings.timeIntegration.verboseMode = 1
127simulationSettings.parallel.numberOfThreads = 8
128
129simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
130simulationSettings.timeIntegration.newton.useModifiedNewton = False
131
132SC.visualizationSettings.general.graphicsUpdateInterval=0.1
133SC.visualizationSettings.general.circleTiling=200
134SC.visualizationSettings.general.drawCoordinateSystem=False
135SC.visualizationSettings.loads.show=False
136SC.visualizationSettings.window.renderWindowSize=[1600,1200]
137SC.visualizationSettings.openGL.multiSampling = 4
138
139
140simulate=True
141if simulate:
142    useGraphics = True
143    if useGraphics:
144        exu.StartRenderer()
145        if 'renderState' in exu.sys:
146            SC.SetRenderState(exu.sys['renderState'])
147        # mbs.WaitForUserToContinue()
148
149    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
150    simulationSettings.timeIntegration.endTime = tEnd
151    simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
152    mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitMidpoint)
153    # print(gContact)
154
155    if useGraphics:
156        SC.WaitForRenderEngineStopFlag()
157        exu.StopRenderer() #safely close rendering window!
158else:
159    SC.visualizationSettings.general.autoFitScene = False
160    SC.visualizationSettings.general.graphicsUpdateInterval=0.5
161
162    sol = LoadSolutionFile('particles.txt')
163    mbs.SolutionViewer(sol)