generalContactSpheresTest.py

You can view and download this file on Github: generalContactSpheresTest.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.utilities import * #includes itemInterface and rigidBodyUtilities
 15import exudyn.graphics as graphics
 16
 17import numpy as np
 18import time
 19
 20testIsActive = exu.sys.get('testIsActive', False)
 21
 22SC = exu.SystemContainer()
 23mbs = SC.AddSystem()
 24
 25nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
 26
 27np.random.seed(1) #always get same results
 28
 29
 30# testIsActive = True
 31
 32L = 1
 33n = 500
 34row = 8 #number of spheres in plane
 35
 36a = L
 37vInit= -20
 38
 39radius = 0.5*a
 40m = 0.05
 41k = 4e4
 42d = 0.001*k*4*0.5*0.2 *0.25
 43markerList = []
 44radiusList = []
 45gDataList = []
 46
 47
 48rb = 30*L
 49H = 8*L
 50Hy=3*L
 51pos0 = [0,-rb-0.5*H,0]
 52pos1 = [-rb-H,-Hy,0]
 53pos2 = [ rb+H,-Hy,0]
 54pos3 = [ 0,-Hy,rb+H]
 55pos4 = [ 0,-Hy,-rb-H]
 56posList=[pos0,pos1,pos2,pos3,pos4]
 57for pos in posList:
 58    colBG = graphics.color.grey
 59    colBG[3] = 0.05
 60    gDataList += [graphics.Sphere(point=pos, radius=rb, color= colBG, nTiles=25)]
 61    nMass = mbs.AddNode(NodePointGround(referenceCoordinates=pos,
 62                        visualization=VNodePointGround(show=False)))
 63
 64    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
 65    markerList += [mThis]
 66    radiusList += [rb]
 67
 68
 69ns = 20
 70gDataSphere = []
 71for i in range(ns):
 72    gRad = radius*(0.75+0.4*(i/ns))
 73    gSphere = graphics.Sphere(point=[0,0,0], radius=gRad, color=graphics.color.blue, nTiles=5)
 74    gDataSphere += [[gSphere]]
 75
 76gDataSphere = []
 77
 78timeCreateStart= -time.time()
 79
 80color4node = graphics.color.blue
 81maxY = 0
 82for i in range(n):
 83
 84    kk = int(i/int(n/16))
 85    color4node = graphics.colorList[min(kk%9,9)]
 86
 87    if (i%20000 == 0): exu.Print("create mass",i)
 88    offy = 0
 89    
 90    iy = int(i/(row*row))
 91    ix = i%row
 92    iz = int(i/row)%row
 93
 94    if iy % 2 == 1:
 95        ix+=0.5
 96        iz+=0.5
 97
 98    offy = -0.25*H-1.5*a+iy*a*0.74 #0.70x is limit value!
 99    offx = -0.6*a-H*0.5 + (ix+1)*a
100    offz = -0.6*a-H*0.5 + (iz+1)*a
101
102    valueRand = np.random.random(1)[0]
103    rFact = 0.2 #random part
104    gRad = radius*(1-rFact+rFact*valueRand)
105    pRef = np.array([offx,offy,offz])
106    maxY = max(maxY,offy)
107    nMass = mbs.AddNode(NodePoint(referenceCoordinates=pRef,
108                                  initialVelocities=[0,vInit,0],
109                                  visualization=VNodePoint(show=True,drawSize=2*gRad, color=color4node)))
110    if i==row*int(row/4)-int(row/2):
111        sNodeNum = nMass
112        if not testIsActive:
113            sNode=mbs.AddSensor(SensorNode(nodeNumber=nMass, fileName='solution/generalContactSpheres.txt',
114                                     outputVariableType=exu.OutputVariableType.Position))
115        
116    oMass = mbs.AddObject(MassPoint(mass=m, nodeNumber=nMass))
117    mThis = mbs.AddMarker(MarkerNodePosition(nodeNumber=nMass))
118    mbs.AddLoad(Force(markerNumber=mThis, loadVector= [0,-m*9.81,0]))
119    markerList += [mThis]
120    radiusList += [gRad]
121
122    mLast = mThis
123
124#put here, such that it is transparent in background
125oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
126                                   visualization=VObjectGround(graphicsData=gDataList)))
127
128exu.Print('generalContactSpheresTest: create bodies:',timeCreateStart+time.time(),'seconds')
129timeCreateStart= -time.time()
130
131if True:
132    gContact = mbs.AddGeneralContact()
133    gContact.verboseMode = 1
134
135    for i in range(len(markerList)):
136        m = markerList[i]
137        r = radiusList[i]
138        gContact.AddSphereWithMarker(m, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
139
140    ssx = 20 #search tree size
141    ssy = 20
142    # mbs.Assemble()
143    # gContact.FinalizeContact(mbs, searchTreeSize=np.array([ssx,ssy,ssx]), frictionPairingsInit=np.eye(1), 
144    #                          searchTreeBoxMin=np.array([-1.2*H,-H,-1.2*H]), searchTreeBoxMax=np.array([1.2*H,14*H,1.2*H]) #80000 particles
145    #                          )
146
147    gContact.SetFrictionPairings(0.*np.eye(1))
148    gContact.SetSearchTreeCellSize(numberOfCells=[ssx,ssy,ssx])
149    gContact.SetSearchTreeBox(pMin=np.array([-1.2*H,-H,-1.2*H]), pMax=np.array([1.2*H,maxY,1.2*H]))
150
151    exu.Print('treesize=',ssx*ssx*ssy)
152
153exu.Print('generalContactSpheresTest: gContact:',timeCreateStart+time.time(),'seconds')
154
155mbs.Assemble()
156exu.Print("finish gContact")
157
158tEnd = 0.1
159h= 0.0002
160simulationSettings = exu.SimulationSettings()
161simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
162simulationSettings.solution.file.write = True
163simulationSettings.solution.file.writePeriod = 0.02
164simulationSettings.solution.sensors.writePeriod = h*10
165simulationSettings.solution.precision = 5 #make files smaller
166simulationSettings.solution.file.export.accelerations = False
167simulationSettings.solution.file.export.velocities = False
168simulationSettings.solution.file.name = 'solution/test.txt'
169simulationSettings.show.computationTime = True
170#simulationSettings.show.statistics = True
171simulationSettings.timeIntegration.verboseMode = 1
172simulationSettings.parallel.numberOfThreads = 1 #use 1 thread to create reproducible results (due to round off errors in sparse vector?)
173
174simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
175simulationSettings.timeIntegration.newton.useModifiedNewton = False
176
177SC.visualizationSettings.general.graphicsUpdateInterval=0.5
178SC.visualizationSettings.general.circleTiling=200
179SC.visualizationSettings.view0.scene.drawCoordinateSystem=False
180SC.visualizationSettings.loads.show=False
181SC.visualizationSettings.bodies.show=True
182SC.visualizationSettings.markers.show=False
183
184SC.visualizationSettings.nodes.show=True
185SC.visualizationSettings.nodes.drawNodesAsPoint = False
186SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
187SC.visualizationSettings.nodes.tiling = 4
188
189SC.visualizationSettings.view0.window.renderWindowSize=[800,800]
190#SC.visualizationSettings.view0.window.renderWindowSize=[1024,1400]
191SC.visualizationSettings.openGL.multiSampling = 4
192#improved OpenGL rendering
193
194
195if not testIsActive:
196    SC.visualizationSettings.general.autoFitScene = False
197    SC.renderer.Start()
198    SC.renderer.RestoreSavedState()
199    SC.renderer.DoIdleTasks()
200
201simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
202simulationSettings.timeIntegration.endTime = tEnd
203simulationSettings.timeIntegration.explicit.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
204simulationSettings.timeIntegration.explicit.computeMassMatrixInversePerBody = True ##2022-12-16: increase performance for multi-threading, Newton increment faster by factor 6 for 8 threads
205
206#one run with a single thread, which is what makes this test reproducible. The scaling over
207#thread counts is measured by generalContactSpheresPerf.py (#2513).
208mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
209
210u = mbs.GetNodeOutput(sNodeNum, exu.OutputVariableType.Coordinates)
211uSum = u[0] + u[1] + u[2]
212exu.Print("u =", u)
213exu.Print('solution of generalContactSpheresTest=',uSum)
214
215
216exu.sys['testResult'] = uSum
217
218
219if not testIsActive:
220    SC.renderer.DoIdleTasks()
221    SC.renderer.Stop() #safely close rendering window!
222
223if not testIsActive:
224    
225    # mbs.PlotSensor([sNode], [2])
226    mbs.PlotSensor([sNode,sNode], [0,1])