generalContactCylinderTest.py

You can view and download this file on Github: generalContactCylinderTest.py

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  model with GeneralContact of a cylinder modelled by spheres rolling on triangle ground mesh
  5#
  6# Author:   Johannes Gerstmayr
  7# Date:     2024-03-16
  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
 16import numpy as np
 17
 18useGraphics = True #without test
 19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 20#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
 21try: #only if called from test suite
 22    from modelUnitTests import exudynTestGlobals #for globally storing test results
 23    useGraphics = exudynTestGlobals.useGraphics
 24except:
 25    class ExudynTestGlobals:
 26        pass
 27    exudynTestGlobals = ExudynTestGlobals()
 28#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 29
 30SC = exu.SystemContainer()
 31mbs = SC.AddSystem()
 32
 33#create an environment for mini example
 34
 35gravity = [0,0,-9.81]
 36
 37planeL = 4
 38p0 = [0.5*planeL,0,0]
 39
 40nPhi = 80  #spheres around circumference
 41nThick = 1 #spheres along cylinder axis (+1)
 42
 43rCyl = 0.25
 44tCyl = 0.1
 45rMarker = 0.02 #radius of spheres for contact points
 46
 47sRad = 0.02
 48k = 1e4
 49d = 0.0005*k
 50stepSize = 1e-3
 51
 52frictionCoeff = 0.2
 53ss = 10
 54
 55markerList = []
 56radiusList = []
 57gDataList = []
 58
 59
 60gContact = mbs.AddGeneralContact()
 61#gContact.verboseMode = 1
 62gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
 63gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
 64# gContact.SetSearchTreeBox(pMin=np.array([-0.5*planeL,-0.5*planeL,-0.1]),
 65#                           pMax=np.array([0.5*planeL,0.5*planeL,0.1]))
 66#print('treesize=',ssx*ssx*ssy)
 67# gContact.sphereSphereContact = False
 68
 69#%% ground
 70gFloor = graphics.CheckerBoard(p0,size=planeL, nTiles=10)
 71gDataList = [gFloor]
 72
 73
 74nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
 75mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
 76
 77[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
 78#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
 79gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
 80    pointList=meshPoints,  triangleList=meshTrigs)
 81
 82
 83#rigid body containing sphere markers:
 84inertia = InertiaCylinder(1000, tCyl, rCyl, 1)
 85# print(inertia)
 86
 87omegaY = 12
 88bCyl=mbs.CreateRigidBody(referencePosition=[0,0,rCyl*1.0],
 89                    initialVelocity=[omegaY*rCyl,0,0],
 90                    initialAngularVelocity=[0,omegaY*0.5,0],
 91                    initialRotationMatrix=RotationMatrixX(0.1),
 92                    inertia=inertia,
 93                    gravity = gravity,
 94                    nodeType = exu.NodeType.RotationRotationVector,
 95                    graphicsDataList=[graphics.Cylinder(pAxis=[0,-0.5*tCyl,0],vAxis=[0,tCyl,0], radius=rCyl,
 96                                                           color=[0.3,0.3,0.3,1],
 97                                                           alternatingColor=graphics.color.lightgrey,nTiles=64)]
 98                    )
 99nCyl = mbs.GetObject(bCyl)['nodeNumber']
100#print(mbs.GetNode(nCyl))
101sPos = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
102                                outputVariableType=exu.OutputVariableType.Position))
103sRot = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
104                                outputVariableType=exu.OutputVariableType.Rotation))
105sVel = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
106                                outputVariableType=exu.OutputVariableType.Velocity))
107sOmega = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
108                                outputVariableType=exu.OutputVariableType.AngularVelocity))
109
110for phiI in range(nPhi):
111    phi = phiI/nPhi*2*pi
112    for j in range(nThick+1):
113        rCylMod = rCyl-rMarker
114        #compute local coordinates for markers
115        y = (j/nThick-0.5)*tCyl
116        x = rCylMod*sin(phi)
117        z = rCylMod*cos(phi)
118
119        m = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bCyl, localPosition=[x,y,z]))
120        gContact.AddSphereWithMarker(m, radius=rMarker, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
121
122
123
124#put ground here, such that it is transparent in background
125oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
126                                    visualization=VObjectGround(graphicsData=gDataList)))
127
128
129mbs.Assemble()
130
131items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
132#print('n spheres=',len(items['MarkerBasedSpheres']))
133
134
135tEnd = 2
136#tEnd = h*100
137simulationSettings = exu.SimulationSettings()
138simulationSettings.solutionSettings.writeSolutionToFile = False
139#simulationSettings.displayComputationTime = True
140#simulationSettings.displayStatistics = True
141simulationSettings.timeIntegration.verboseMode = 1
142#simulationSettings.parallel.numberOfThreads = 4
143
144SC.visualizationSettings.general.graphicsUpdateInterval=0.02
145SC.visualizationSettings.general.drawCoordinateSystem=True
146SC.visualizationSettings.loads.show=False
147SC.visualizationSettings.bodies.show=True
148SC.visualizationSettings.markers.show=False
149
150SC.visualizationSettings.nodes.show=True
151SC.visualizationSettings.nodes.drawNodesAsPoint = False
152SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
153SC.visualizationSettings.nodes.tiling = 4
154
155SC.visualizationSettings.window.renderWindowSize=[2000,1200]
156SC.visualizationSettings.openGL.multiSampling = 4
157#improved OpenGL rendering
158
159SC.visualizationSettings.contact.showSpheres = True
160
161SC.visualizationSettings.general.autoFitScene = False
162
163if useGraphics:
164    exu.StartRenderer()
165    if 'renderState' in exu.sys:
166        SC.SetRenderState(exu.sys['renderState'])
167    # mbs.WaitForUserToContinue()
168
169
170simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
171simulationSettings.timeIntegration.endTime = tEnd
172mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
173
174if useGraphics:
175    #SC.WaitForRenderEngineStopFlag()
176    exu.StopRenderer() #safely close rendering window!
177
178#%%+++++++++++++++++++
179q = mbs.GetSensorValues(sPos)
180q += mbs.GetSensorValues(sVel)
181q += mbs.GetSensorValues(sOmega)
182q += mbs.GetSensorValues(sRot)
183#print('q=', q)
184
185u = NormL2(q)
186exu.Print('solution of generalContactCylinderTest =',u)
187
188exudynTestGlobals.testError = u - (12.42377622187738 )
189exudynTestGlobals.testResult = u