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