You can view and download this file on Github: generalContactFrictionTests.py
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# This is an EXUDYN example
#
# Details: test friction of spheres and spheres-trigs
#
# Author: Johannes Gerstmayr
# Date: 2021-12-06
#
# 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 *
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()
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0]))
#isPerformanceTest = exudynTestGlobals.isPerformanceTest
#useGraphics = False
# isPerformanceTest = True
# tEnd = 0.1
# if isPerformanceTest: tEnd *= 0.5
#%%+++++++++++++++++++++++++++++++++
#sphere-sphere with coordinate constraints, prestressed; fixed torque on one side, linear increasing torque on other side
#sphere on ground, rolling
#cube on ground, sliding (f=[f, f*mu*t, 0]), tangential force changing
#cube on ground with initial velocity
#cube-cube contact (meshed)
L = 1 #surrounding
a = 0.1 #base dimention of objects
r = 0.5*a #radius
t = 0.25*a #thickness
#contact coefficients:
mu = 0.8 #dry friction
m = 0.025 #mass
k = 1e3 #(linear) normal contact stiffness
d = 1e-4*k #(linear) contact damping
gFact = 10
g = [0,0,-gFact]
gContact = mbs.AddGeneralContact()
gContact.verboseMode = 1
#gContact.sphereSphereContact = False
gContact.frictionProportionalZone = 1e-3
#gContact.excludeDuplicatedTrigSphereContactPoints = False
fricMat = mu*np.eye(2)
fricMat[0,1] = 0.2
fricMat[1,0] = 0.2
gContact.SetFrictionPairings(fricMat)
gContact.SetSearchTreeCellSize(numberOfCells=[4,4,4])
#%% ground
p0 = np.array([0,0,-0.5*t])
color4wall = [0.9,0.9,0.7,0.5]
addNormals = False
gFloor = GraphicsDataOrthoCubePoint(p0,[L,L,t],color4steelblue,addNormals)
gFloorAdd = GraphicsDataOrthoCubePoint(p0+[-0.5*L,0,a],[t,L,2*a],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gFloorAdd = GraphicsDataOrthoCubePoint(p0+[ 0.5*L,0,a],[t,L,2*a],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gFloorAdd = GraphicsDataOrthoCubePoint(p0+[0,-0.5*L,a],[L,t,2*a],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gFloorAdd = GraphicsDataOrthoCubePoint(p0+[0, 0.5*L,a],[L,t,2*a],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
bb = 0.75*a
bh = 0.25*a
p1 = np.array([0.5*L,0.5*L,0])
gFloorAdd = GraphicsDataOrthoCubePoint(p1+[-bb*3, -bb, 0.5*bh],[bb,bb,bh],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gFloorAdd = GraphicsDataOrthoCubePoint(p1+[-bb*2, -bb, 1.5*bh],[bb,bb,bh],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gFloorAdd = GraphicsDataOrthoCubePoint(p1+[-bb*1, -bb, 2.5*bh],[bb,bb,bh],color4wall,addNormals)
gFloor = MergeGraphicsDataTriangleList(gFloor, gFloorAdd)
gDataList = [gFloor]
nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0] ))
mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
mGroundC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
[meshPoints, meshTrigs] = GraphicsData2PointsAndTrigs(gFloor)
#[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
# [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)
if True: #looses color
gFloor = GraphicsDataFromPointsAndTrigs(meshPoints, meshTrigs, color=color4wall) #show refined mesh
gDataList = [gFloor]
evalNodes = [] #collect nodes that are evaluated for test
#%%++++++++++++++++++++++++++++++++++++++++++++
#free rolling sphere:
gList = [GraphicsDataSphere(point=[0,0,0], radius=r, color= color4red, nTiles=24)]
omega0 = -4.*np.array([5,1.,0.])
pRef = [-0.4*L,-0.4*L,r-0*m*gFact/k]
RBinertia = InertiaSphere(m, r)
[nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
velocity=-np.cross([0,0,r],omega0),
rotationMatrix=RotationMatrixX(0.),
angularVelocity=omega0,
#gravity=g,
graphicsDataList=gList,
)
nNode0 = nMass
mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
mbs.AddLoad(Force(markerNumber=mNode, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
exu.Print('expect u0z=',2*r*0.01)
gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
if useGraphics:
sNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0.txt',
outputVariableType=exu.OutputVariableType.Displacement))
vNode0 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
outputVariableType=exu.OutputVariableType.Velocity))
evalNodes += [nMass]
#%%++++++++++++++++++++++++++++++++++++++++++++
#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
gList = [GraphicsDataSphere(point=[0,0,0], radius=r, color= color4yellow, nTiles=24)]
omega0 = -1e-12*np.array([1,0.1,0.])
pRef = [1e-15,-1e-14,r-2*m*gFact/k]
RBinertia = InertiaSphere(m, r)
[nMass, oMass] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
velocity=-np.cross([0,0,r],omega0),
rotationMatrix=RotationMatrixX(0.),
angularVelocity=omega0,
gravity=g,
graphicsDataList=gList,
)
nNode1 = nMass
mNode1 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMass))
#mbs.AddLoad(Force(markerNumber=mNode1, loadVector= [0,0,-k*r*0.01])) #==> uz = 2*r*0.01
#exu.Print('expect u0z=',2*r*0.01)
gContact.AddSphereWithMarker(mNode1, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
sNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode1, storeInternal=True, #fileName='solution/contactNode1.txt',
outputVariableType=exu.OutputVariableType.Displacement))
# vNode1 = mbs.AddSensor(SensorNode(nodeNumber=nNode0, storeInternal=True, #fileName='solution/contactNode0Vel.txt',
# outputVariableType=exu.OutputVariableType.Velocity))
#%%++++++++++++++++++++++++++++++++++++++++++++
#fixed pressure tests:
pf = np.array([-1.2*L,0,0])
nGroundF = mbs.AddNode(NodePointGround(referenceCoordinates=pf ))
mNode = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF))
gContact.AddSphereWithMarker(mNode, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
gDataList += [GraphicsDataSphere(point=pf, radius=r, color= color4grey, nTiles=24)]
gList = [GraphicsDataSphere(point=[0,0,0], radius=r, color= color4lightgreen, nTiles=24)]
pRef = pf+[0,2*r,0] #[-0.4*L,-0.4*L,r-m*gFact/k]
RBinertia = InertiaSphere(m, r)
[nMassF, oMassF] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
rotationMatrix=RotationMatrixX(0.),
#gravity=g,
graphicsDataList=gList,
)
mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassF, coordinate=0))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
mNodeF = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassF))
mbs.AddLoad(Force(markerNumber=mNodeF, loadVector= [0,-k*r*0.1,0])) #==> u = k*r*0.1/(0.5*k) = 2*r*0.1
exu.Print('expect uFy=',2*r*0.1)
gContact.AddSphereWithMarker(mNodeF, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
if useGraphics:
sNodeF = mbs.AddSensor(SensorNode(nodeNumber=nMassF, storeInternal=True, #fileName='solution/contactNodeF.txt',
outputVariableType=exu.OutputVariableType.Displacement))
evalNodes += [nMassF]
#%%++++++++++++++++++++++++++++++++++++++++++++
# sliding between spheres:
pr = np.array([-1.2*L,0.5*L,0])
nGroundF2 = mbs.AddNode(NodePointGround(referenceCoordinates=pr ))
mNode2 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGroundF2))
gContact.AddSphereWithMarker(mNode2, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
gDataList += [GraphicsDataSphere(point=pr, radius=r, color= color4lightgrey, nTiles=24)]
gList = [GraphicsDataSphere(point=[0,0,0], radius=r, color= color4lightred, nTiles=24)]
dRol = r*0.01
pRef = pr+[0,2*r-2*dRol,0] #force=k*r*0.01
fRol = k*dRol
exu.Print('force rolling=', fRol, ', torque=', fRol*mu*r)
RBinertia = InertiaSphere(m, r)
[nMassR, oMassR] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
rotationMatrix=RotationMatrixX(0.),
#angularVelocity=[10,0,0],
#gravity=g,
graphicsDataList=gList,
)
mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=0))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=1))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
mC = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nMassR, coordinate=2))
mbs.AddObject(CoordinateConstraint(markerNumbers=[mGroundC, mC]))
mNodeR = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassR))
def UFtorque(mbs, t, loadVector):
torque = 10*t*fRol*mu*r
if t > 0.3:
torque = 0
return [torque,0,0]
mbs.AddLoad(Torque(markerNumber=mNodeR, loadVectorUserFunction=UFtorque,
loadVector= [1,0,0])) #==> u = k*r*0.1/(0.5*k) = 2*r*0.1
gContact.AddSphereWithMarker(mNodeR, radius=r, contactStiffness=k, contactDamping=d, frictionMaterialIndex=0)
if useGraphics:
sNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeR.txt',
outputVariableType=exu.OutputVariableType.Rotation))
vNodeR = mbs.AddSensor(SensorNode(nodeNumber=nMassR, storeInternal=True, #fileName='solution/contactNodeRvel.txt',
outputVariableType=exu.OutputVariableType.AngularVelocity))
evalNodes += [nMassR]
#%%++++++++++++++++++++++++++++++++++++++++++++
#sphere on stairs
#%%++++++++++++++++++++++++++++++++++++++++++++
#free rolling sphere at midpoint, many triangles in close contact; slowly go through critical points:
gList = [GraphicsDataSphere(point=[0,0,0], radius=0.5*r, color= color4yellow, nTiles=24)]
omega0 = np.array([-0.05,-5,0.])
pRef = [0.5*L-1.45*bb, 0.5*L-1.20*bb, 3*bh+0.5*r-2*m*gFact/k] #[0.5*L-1.45*bb, 0.5*L-1.40*bb, ..] goes to edge
RBinertia = InertiaSphere(m, 0.5*r)
[nMassStair, oMassStair] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
velocity=-np.cross([0,0,0.5*r],omega0),
rotationMatrix=RotationMatrixX(0.),
angularVelocity=omega0,
gravity=g,
graphicsDataList=gList,
)
nNode3 = nMassStair
mNode3 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassStair))
gContact.AddSphereWithMarker(mNode3, radius=0.5*r, contactStiffness=k, contactDamping=20*d, frictionMaterialIndex=0)
if useGraphics:
sNode3 = mbs.AddSensor(SensorNode(nodeNumber=nNode3, storeInternal=True, #fileName='solution/contactNode3.txt',
outputVariableType=exu.OutputVariableType.Displacement))
evalNodes += [nMassStair]
#%%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#contact of cube with ground
tTrig = 0.25*r #size of contact points on mesh ('thickness')
gCube = GraphicsDataOrthoCubePoint(size=[3*r,2*r,r], color= color4steelblue,addNormals=addNormals)
[meshPoints, meshTrigs] = GraphicsData2PointsAndTrigs(gCube)
#for tests, 1 refinement!
[meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
# exu.Print("n points=",len(meshPoints))
# [meshPoints, meshTrigs] = RefineMesh(meshPoints, meshTrigs) #just to have more triangles on floor
# exu.Print("==> n points refined=",len(meshPoints))
# refinements give 26,98,386 points!
[meshPoints2, meshTrigs2] = ShrinkMeshNormalToSurface(meshPoints, meshTrigs, tTrig)
#add mesh to visualization
gCube = GraphicsDataFromPointsAndTrigs(meshPoints, meshTrigs, color=color4steelblue) #show refined mesh
gList = [gCube]
#add points for contact to visualization (shrinked)
for p in meshPoints2:
gList += [GraphicsDataSphere(point=p, radius=tTrig, color=color4red)]
pRef = [0.5*L-2*r, 0.25*L, 0.5*r+1.5*tTrig]
v0 = np.array([-2,0,0])
RBinertia = InertiaCuboid(density=m/(r*2*r*3*r), sideLengths=[3*r,2*r,r])
[nMassCube0, oMassCube0] = AddRigidBody(mainSys=mbs, inertia=RBinertia,
nodeType=exu.NodeType.RotationRotationVector,
position=pRef,
velocity=v0,
#rotationMatrix=RotationMatrixZ(0.),
angularVelocity=[0,0,0],
gravity=g,
graphicsDataList=gList,
)
nCube0 = nMassCube0
mCube0 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nMassCube0))
gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mCube0, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1,
pointList=meshPoints, triangleList=meshTrigs)
for p in meshPoints2:
mPoint = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oMassCube0, localPosition=p))
gContact.AddSphereWithMarker(mPoint, radius=tTrig, contactStiffness=k, contactDamping=d, frictionMaterialIndex=1)
if useGraphics:
sCube0 = mbs.AddSensor(SensorNode(nodeNumber=nCube0, storeInternal=True, #fileName='solution/contactCube0.txt',
outputVariableType=exu.OutputVariableType.Displacement))
evalNodes += [nMassCube0]
#%%++++++++++++++++++++++++++++++++++++++++++++
#add as last because of transparency
oGround = mbs.AddObject(ObjectGround(visualization=VObjectGround(graphicsData=gDataList)))
#%%+++++++++++++++++++++++++++++++++
mbs.Assemble()
tEnd = 0.8 #tEnd = 0.8 for test suite
h= 0.0002 #h= 0.0002 for test suite
# h*=0.1
# tEnd*=3
simulationSettings = exu.SimulationSettings()
#simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
simulationSettings.solutionSettings.writeSolutionToFile = False
if useGraphics:
simulationSettings.solutionSettings.solutionWritePeriod = 0.001
simulationSettings.solutionSettings.writeSolutionToFile = True
simulationSettings.solutionSettings.coordinatesSolutionFileName = 'solution/coordinatesSolution.txt'
else:
simulationSettings.solutionSettings.exportAccelerations = False
simulationSettings.solutionSettings.exportVelocities = False
simulationSettings.solutionSettings.sensorsWritePeriod = h*10
simulationSettings.solutionSettings.outputPrecision = 8 #make files smaller
# simulationSettings.displayComputationTime = True
# simulationSettings.displayGlobalTimers = True
#simulationSettings.displayStatistics = True
simulationSettings.timeIntegration.verboseMode = 1
simulationSettings.parallel.numberOfThreads = 1
simulationSettings.timeIntegration.newton.numericalDifferentiation.forODE2 = False
simulationSettings.timeIntegration.newton.useModifiedNewton = False
SC.visualizationSettings.general.graphicsUpdateInterval=0.05
# SC.visualizationSettings.general.drawWorldBasis = True
SC.visualizationSettings.general.circleTiling=200
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.showBasis =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.openGL.drawFaceNormals = False
SC.visualizationSettings.openGL.multiSampling = 4
SC.visualizationSettings.openGL.shadow = 0.25
SC.visualizationSettings.openGL.light0position = [-3,3,10,0]
if useGraphics:
SC.visualizationSettings.general.autoFitScene = False
exu.StartRenderer()
if 'renderState' in exu.sys:
SC.SetRenderState(exu.sys['renderState'])
mbs.WaitForUserToContinue()
simulationSettings.timeIntegration.numberOfSteps = int(tEnd/h)
simulationSettings.timeIntegration.endTime = tEnd
simulationSettings.timeIntegration.explicitIntegration.computeEndOfStepAccelerations = False #increase performance, accelerations less accurate
mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ExplicitEuler)
# mbs.SolveDynamic(simulationSettings, solverType=exu.DynamicSolverType.ODE23)
#compute error:
uSum=0
for node in evalNodes:
u = mbs.GetNodeOutput(node, exu.OutputVariableType.Coordinates)
exu.Print('coords node'+str(node)+' =',u)
for c in u:
uSum += abs(c) #add up all coordinates for comparison
exu.Print('solution of generalContactFrictionTest=',uSum)
exudynTestGlobals.testError = uSum - (10.132106712933348 )
exudynTestGlobals.testResult = uSum
if useGraphics:
SC.WaitForRenderEngineStopFlag()
if True:
SC.visualizationSettings.general.autoFitScene = False
SC.visualizationSettings.general.graphicsUpdateInterval=0.02
sol = LoadSolutionFile('solution/coordinatesSolution.txt', safeMode=True)#, maxRows=100)
print('start SolutionViewer')
mbs.SolutionViewer(sol)
exu.StopRenderer() #safely close rendering window!
if useGraphics:
mbs.PlotSensor([], closeAll=True)
mbs.PlotSensor([sNode3]*3, [0,1,2], figureName='node stair')