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