You can view and download this file on Github: generalContactCylinderTest.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: model with GeneralContact of a cylinder modelled by spheres rolling on triangle ground mesh
#
# Author: Johannes Gerstmayr
# Date: 2024-03-16
#
# 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.
#
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
import exudyn as exu
from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
import exudyn.graphics as graphics #only import if it does not conflict
import numpy as np
useGraphics = True #without test
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
try: #only if called from test suite
from modelUnitTests import exudynTestGlobals #for globally storing test results
useGraphics = exudynTestGlobals.useGraphics
except:
class ExudynTestGlobals:
pass
exudynTestGlobals = ExudynTestGlobals()
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SC = exu.SystemContainer()
mbs = SC.AddSystem()
#create an environment for mini example
gravity = [0,0,-9.81]
planeL = 4
p0 = [0.5*planeL,0,0]
nPhi = 80 #spheres around circumference
nThick = 1 #spheres along cylinder axis (+1)
rCyl = 0.25
tCyl = 0.1
rMarker = 0.02 #radius of spheres for contact points
sRad = 0.02
k = 1e4
d = 0.0005*k
stepSize = 1e-3
frictionCoeff = 0.2
ss = 10
markerList = []
radiusList = []
gDataList = []
gContact = mbs.AddGeneralContact()
#gContact.verboseMode = 1
gContact.SetFrictionPairings(frictionCoeff*np.eye(1))
gContact.SetSearchTreeCellSize(numberOfCells=[ss,ss,ss])
# gContact.SetSearchTreeBox(pMin=np.array([-0.5*planeL,-0.5*planeL,-0.1]),
# pMax=np.array([0.5*planeL,0.5*planeL,0.1]))
#print('treesize=',ssx*ssx*ssy)
# gContact.sphereSphereContact = False
#%% ground
gFloor = graphics.CheckerBoard(p0,size=planeL, nTiles=10)
gDataList = [gFloor]
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
[meshPoints, meshTrigs] = graphics.ToPointsAndTrigs(gFloor)
#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mGround, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0,
pointList=meshPoints, triangleList=meshTrigs)
#rigid body containing sphere markers:
inertia = InertiaCylinder(1000, tCyl, rCyl, 1)
# print(inertia)
omegaY = 12
bCyl=mbs.CreateRigidBody(referencePosition=[0,0,rCyl*1.0],
initialVelocity=[omegaY*rCyl,0,0],
initialAngularVelocity=[0,omegaY*0.5,0],
initialRotationMatrix=RotationMatrixX(0.1),
inertia=inertia,
gravity = gravity,
nodeType = exu.NodeType.RotationRotationVector,
graphicsDataList=[graphics.Cylinder(pAxis=[0,-0.5*tCyl,0],vAxis=[0,tCyl,0], radius=rCyl,
color=[0.3,0.3,0.3,1],
alternatingColor=graphics.color.lightgrey,nTiles=64)]
)
nCyl = mbs.GetObject(bCyl)['nodeNumber']
#print(mbs.GetNode(nCyl))
sPos = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
outputVariableType=exu.OutputVariableType.Position))
sRot = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
outputVariableType=exu.OutputVariableType.Rotation))
sVel = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
outputVariableType=exu.OutputVariableType.Velocity))
sOmega = mbs.AddSensor(SensorBody(bodyNumber=bCyl, storeInternal=True,
outputVariableType=exu.OutputVariableType.AngularVelocity))
for phiI in range(nPhi):
phi = phiI/nPhi*2*pi
for j in range(nThick+1):
rCylMod = rCyl-rMarker
#compute local coordinates for markers
y = (j/nThick-0.5)*tCyl
x = rCylMod*sin(phi)
z = rCylMod*cos(phi)
m = mbs.AddMarker(MarkerBodyRigid(bodyNumber=bCyl, localPosition=[x,y,z]))
gContact.AddSphereWithMarker(m, radius=rMarker, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
#put ground here, such that it is transparent in background
oGround=mbs.AddObject(ObjectGround(referencePosition= [0,0,0],
visualization=VObjectGround(graphicsData=gDataList)))
mbs.Assemble()
items=gContact.GetItemsInBox(pMin=[-4,-4,0], pMax=[4,4,20])
#print('n spheres=',len(items['MarkerBasedSpheres']))
tEnd = 2
#tEnd = h*100
simulationSettings = exu.SimulationSettings()
simulationSettings.solutionSettings.writeSolutionToFile = False
#simulationSettings.displayComputationTime = True
#simulationSettings.displayStatistics = True
simulationSettings.timeIntegration.verboseMode = 1
# simulationSettings.parallel.numberOfThreads = 4
# simulationSettings.timeIntegration.simulateInRealtime = True
SC.visualizationSettings.general.graphicsUpdateInterval=0.02
SC.visualizationSettings.general.drawCoordinateSystem=True
SC.visualizationSettings.loads.show=False
SC.visualizationSettings.bodies.show=True
SC.visualizationSettings.markers.show=False
SC.visualizationSettings.nodes.show=True
SC.visualizationSettings.nodes.drawNodesAsPoint = False
SC.visualizationSettings.nodes.defaultSize = 0 #must not be -1, otherwise uses autocomputed size
SC.visualizationSettings.nodes.tiling = 4
SC.visualizationSettings.window.renderWindowSize=[2000,1200]
SC.visualizationSettings.openGL.multiSampling = 4
#improved OpenGL rendering
SC.visualizationSettings.contact.showSpheres = True
SC.visualizationSettings.general.autoFitScene = False
if useGraphics:
exu.StartRenderer()
if 'renderState' in exu.sys:
SC.SetRenderState(exu.sys['renderState'])
# mbs.WaitForUserToContinue()
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize)
simulationSettings.timeIntegration.endTime = tEnd
mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
if useGraphics:
#SC.WaitForRenderEngineStopFlag()
exu.StopRenderer() #safely close rendering window!
#%%+++++++++++++++++++
q = mbs.GetSensorValues(sPos)
q += mbs.GetSensorValues(sVel)
q += mbs.GetSensorValues(sOmega)
q += mbs.GetSensorValues(sRot)
#print('q=', q)
u = NormL2(q)
exu.Print('solution of generalContactCylinderTest =',u)
exudynTestGlobals.testError = u - (12.42377622187738 )
exudynTestGlobals.testResult = u