contactComparisonTest.py
You can view and download this file on Github: contactComparisonTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Three contact objects, one drop: a ball of radius r falls onto the ground with the same
5# linear penalty law (contact stiffness k, damping d). ObjectContactCoordinate sees it as
6# one coordinate, ObjectContactSphereSphere as a hit on a large ground sphere,
7# ObjectContactSphereTriangle as a hit on a ground triangle - the motion is the same, so the
8# three give the same height z(t) (#2749), a cross-check of three implementations.
9# Measured: all three agree to 1e-12 over the whole drop - they recommend the same step
10# size where the contact begins and where it ends (#2109); the deepest point 0.0963898
11# against 0.0963927 of an independent RK4 integration of the same law.
12# The same for a nonlinear law with the impact model of Gonthier et al. / Carvalho-Martins,
13# which ObjectContactCoordinate has since #2750.
14# A fourth case is GeneralContact, a sphere on a ground triangle, implicit (#1848): its
15# stiffnesses act in series and its dampings in parallel, so each gets 2k and d/2. Without
16# damping it agrees with the others to 1e-15; with damping it releases the ball earlier and
17# it rebounds higher, because its implicit solver keeps a contact only while the contact
18# force presses - the contact objects, and GeneralContact with an explicit solver, keep it
19# while the ball penetrates and let the damping pull: identical until the release, then a
20# rebound to 0.1391 instead of 0.1352. With keepContactWhilePenetrating = True (#2848) it does
21# what the contact objects do, the fifth case: identical until the release, a rebound to 0.13511
22# instead of 0.13523, and both converge to 0.1351859 with smaller steps (equal to 2e-15 with
23# 80000 steps) - they differ only in the step where the contact ends.
24# Then a ball that slides with friction (#1947): ObjectContactSphereTriangle and
25# GeneralContact, implicit, against the rigid-body solution - the ball rolls after
26# t = 2 v0/(7 mu g) with 5/7 of its initial velocity. Measured: both to 3e-7 while it
27# slides, to 6e-4 when it rolls (the regularized friction and the penetration).
28#
29# Author: Johannes Gerstmayr
30# Date: 2026-09-29
31#
32# 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.
33#
34#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
35
36import exudyn as exu
37from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
38
39import numpy as np
40
41testIsActive = exu.sys.get('testIsActive', False)
42
43r = 0.1 #radius of the ball
44m = 1 #mass
45k = 1e5 #contact stiffness
46d = 100 #contact damping
47g = 9.81
48z0 = 0.2 #initial height of the center: falls 0.1 before contact
49tEnd = 0.3 #one contact and the rebound
50nSteps = 20000
51
52
53trianglePoints = [[-1, -1, 0], [2, -1, 0], [0, 2, 0]] #the ground triangle, around the path of the sliding ball
54
55
56def Drop(kind, law, z0=z0, vx=0, mu=0, frictionProportionalZone=1e-3, sensorType=exu.OutputVariableType.Position,
57 component=3):
58 """a component of a sensor of the ball over time - the height of its center by default -, for one of the four
59 contact cases, a contact law (a dict of the law's parameters), and a ball that starts at height z0 with the
60 velocity vx and slides with the friction coefficient mu"""
61 nImpact = 1 if law.get('impactModel', 0) != 0 else 0
62 SC = exu.SystemContainer()
63 mbs = SC.AddSystem()
64 oGround = mbs.AddObject(ObjectGround())
65 if kind == 'coordinate':
66 #the coordinate is the height; gap = q - offset
67 node = mbs.AddNode(Node1D(referenceCoordinates=[0], initialCoordinates=[z0]))
68 mbs.AddObject(ObjectMass1D(nodeNumber=node, mass=m))
69 mBall = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=node, coordinate=0))
70 mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
71 nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=1+nImpact, initialCoordinates=[z0 - r]+[0]*nImpact))
72 mbs.AddObject(ObjectContactCoordinate(markerNumbers=[mGround, mBall], nodeNumber=nData,
73 offset=r, **law))
74 mbs.AddLoad(LoadCoordinate(markerNumber=mBall, load=-m*g))
75 sensor = mbs.AddSensor(SensorNode(nodeNumber=node, outputVariableType=exu.OutputVariableType.Coordinates,
76 storeInternal=True, writeToFile=False))
77 component = 1
78 else:
79 node = mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[0.1, 0.1, z0] + eulerParameters0,
80 initialVelocities=[vx, 0, 0, 0, 0, 0, 0]))
81 ball = mbs.AddObject(ObjectRigidBody(nodeNumber=node, mass=m,
82 inertia=InertiaSphere(mass=m, radius=r).GetInertia6D()))
83 mBall = mbs.AddMarker(MarkerBodyRigid(bodyNumber=ball))
84 mbs.AddLoad(LoadForceVector(markerNumber=mBall, loadVector=[0, 0, -m*g]))
85 friction = {'dynamicFriction': mu, 'frictionProportionalZone': frictionProportionalZone} if mu != 0 else {}
86 nData = mbs.AddNode(NodeGenericData(numberOfDataCoordinates=4, initialCoordinates=[0, 0, 0, 0]))
87 if kind in ['general', 'generalPenetrating']:
88 #the linear law only; stiffnesses in series, dampings in parallel
89 gContact = mbs.AddGeneralContact()
90 gContact.keepContactWhilePenetrating = (kind == 'generalPenetrating')
91 gContact.frictionProportionalZone = frictionProportionalZone
92 gContact.SetFrictionPairings(mu*np.eye(1))
93 gContact.SetSearchTreeCellSize(numberOfCells=[1, 1, 1])
94 gContact.AddTrianglesRigidBodyBased(rigidBodyMarkerIndex=mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround)),
95 contactStiffness=2*law['contactStiffness'],
96 contactDamping=0.5*law['contactDamping'], frictionMaterialIndex=0,
97 pointList=trianglePoints, triangleList=[[0, 1, 2]])
98 gContact.AddSphereWithMarker(mbs.AddMarker(MarkerNodeRigid(nodeNumber=node)), radius=r,
99 contactStiffness=2*law['contactStiffness'],
100 contactDamping=0.5*law['contactDamping'], frictionMaterialIndex=0)
101 elif kind == 'sphere':
102 R = 1 #a ground sphere whose top is at z=0, right below the ball
103 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround, localPosition=[0.1, 0.1, -R]))
104 mbs.AddObject(ObjectContactSphereSphere(markerNumbers=[mGround, mBall], nodeNumber=nData,
105 spheresRadii=[R, r], **law, **friction))
106 else:
107 #the sphere is marker 0, the triangle marker 1
108 mGround = mbs.AddMarker(MarkerBodyRigid(bodyNumber=oGround))
109 mbs.AddObject(ObjectContactSphereTriangle(markerNumbers=[mBall, mGround], nodeNumber=nData, sphereRadius=r,
110 trianglePoints=exu.Vector3DList(trianglePoints),
111 **law, **friction))
112 sensor = mbs.AddSensor(SensorNode(nodeNumber=node, outputVariableType=sensorType,
113 storeInternal=True, writeToFile=False))
114 mbs.Assemble()
115 simulationSettings = exu.SimulationSettings()
116 simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
117 simulationSettings.timeIntegration.numberOfSteps = nSteps
118 simulationSettings.timeIntegration.endTime = tEnd
119 simulationSettings.timeIntegration.verboseMode = 0
120 simulationSettings.solution.sensors.writePeriod = tEnd/nSteps
121 simulationSettings.solution.file.write = False
122 mbs.SolveDynamic(simulationSettings)
123 data = mbs.GetSensorStoredData(sensor)
124 #on one time grid: the adaptive step control puts the sensor times of the release step at different places
125 return np.interp(np.linspace(0, tEnd, nSteps+1), data[:, 0], data[:, component])
126
127
128#the linear law, and a nonlinear law with the impact model of Gonthier et al. / Carvalho-Martins (#2750)
129laws = {'linear': {'contactStiffness': k, 'contactDamping': d},
130 'impact': {'contactStiffness': k, 'contactDamping': 0, 'contactStiffnessExponent': 1.5,
131 'impactModel': 2, 'restitutionCoefficient': 0.5}}
132t = np.linspace(0, tEnd, nSteps+1)
133testResult = 0
134for (lawName, law) in laws.items():
135 kinds = ['coordinate', 'sphere', 'triangle'] + (['general', 'generalPenetrating'] if lawName == 'linear' else [])
136 z = dict((kind, Drop(kind, law)) for kind in kinds)
137 release = np.argmax((t > 0.14) & (z['coordinate'] > r)) #the ball leaves the ground again
138 exu.Print(lawName + ':')
139 exu.Print(' deepest point :', z['coordinate'].min(), z['sphere'].min(), z['triangle'].min())
140 exu.Print(' sphere-triangle :', np.abs(z['sphere'] - z['triangle']).max())
141 if 'general' in z:
142 exu.Print(' GeneralContact : until release', np.abs(z['general'][:release] - z['sphere'][:release]).max(),
143 ', rebound to', z['general'][release:].max(), 'instead of', z['sphere'][release:].max())
144 exu.Print(' GeneralContact, keepContactWhilePenetrating: until release',
145 np.abs(z['generalPenetrating'][:release] - z['sphere'][:release]).max(), ', rebound to',
146 z['generalPenetrating'][release:].max(), ', largest difference', np.abs(z['generalPenetrating'] - z['sphere']).max())
147 exu.Print(' until release :', np.abs(z['coordinate'][:release] - z['sphere'][:release]).max())
148 exu.Print(' after release :', np.abs(z['coordinate'] - z['sphere']).max())
149 #the ratio of the velocities leaving and hitting the ground, from the rebound height; gravity acts during the contact
150 exu.Print(' rebound ratio :', np.sqrt(2*g*(z['sphere'][release:].max() - r))/np.sqrt(2*g*(z0 - r)))
151 #the deepest points of the three and the heights at the end - moves if any of the three changes
152 testResult += sum(z[kind].min() + z[kind][-1] for kind in z)
153
154#sliding with friction: the ball starts on the ground with the velocity v0 and no rotation; it slides until
155#t = 2 v0/(7 mu g) and then rolls with 5/7 v0
156v0 = 1
157mu = 0.2
158tRoll = 2*v0/(7*mu*g)
159zRest = r - m*g/k #the static penetration
160for kind in ['triangle', 'general']:
161 vx = Drop(kind, laws['linear'], z0=zRest, vx=v0, mu=mu, sensorType=exu.OutputVariableType.Velocity, component=1)
162 iHalf = np.argmax(t >= 0.5*tRoll)
163 exu.Print(kind + ': velocity while sliding', vx[iHalf], 'instead of', v0 - mu*g*t[iHalf],
164 ', when rolling', vx[-1], 'instead of', 5/7*v0)
165 testResult += vx[iHalf] + vx[-1]
166
167exu.Print('solution of contactComparisonTest=', testResult)
168exu.sys['testResult'] = testResult