Skip to content

Commit 5a7a3c4

Browse files
Merge pull request #76 from ModiaSim/mo_touching_contact
Contact hysteresis strategy slightly changed (and docu adapted + requ…
2 parents c688973 + ecc5d2e commit 5a7a3c4

File tree

4 files changed

+55
-37
lines changed

4 files changed

+55
-37
lines changed

docs/src/internal/ContactDetection.md

Lines changed: 42 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -23,59 +23,72 @@ otherwise a zero crossing cannot be detected.
2323

2424
The distance between two objects computed either in the narrow phase or in an approximate
2525
way in the broad phase (= distance between Axis Aligned Bounding Boxes) is called `distanceOrg`.
26-
If `distanceOrg < 0`, then the two objects are penetrating each other.
26+
If `distanceOrg < 0`, then the two objects are penetrating each other. The crossing functions are
27+
formulated with the help of a hysteresis defined as:
28+
29+
```julia
30+
contact_eps = max(1e-13, 10*mprTolerance)
31+
```
32+
33+
So it is in the order of ``10^{-13} m``, but at least a factor of 10 larger as the tolerance
34+
used for the distance computation with the MPR algorithm.
35+
36+
If the two shapes are treated in contact, the distance actually used for the elastic
37+
response calculation is
38+
39+
```
40+
distance = distanceOrg + contact_eps
41+
```
42+
43+
This means that penetration is assumed to occur at `distanceOrg < -contact_eps`, in order
44+
that there is by sure penetration, although `distanceOrg` is computed with some error.
2745

2846
Modia3D uses the dictionary `contactDict` to keep track of the
2947
contact situation. Every pair of objects is identified by a unique
3048
Integer value called *PairID* that is used as key in `conctactDict`.
3149
A dictionary value is an instance of [`Modia3D.ContactPair`](@ref)
3250

3351
At an **event instant** (including initialization), dictionary `contactDict` is emptied
34-
and all object pairs are stored newly in `contactDict` that have `distanceOrg <= 0`, so are
35-
either touching or penetrating.
52+
and all object pairs are stored newly in `contactDict` that have `distanceOrg < -2*contact_eps`, so are
53+
penetrating already a bit. Note, this is useful, for example in case of a gripper, where it is natural that two shapes
54+
have a distance of zero and no contact situation is present.
3655

3756
During **continuous integration**, two zero crossing functions are used:
3857

39-
* ``z_1(t)``: The maximum of `distanceOrg` ``- 10^{-16}`` of all object pairs that are in `contactDict`.
58+
* ``z_1(t)``: The maximum of `distanceOrg+contact_eps` of all object pairs that are in `contactDict`.
4059
``z_1`` monitors if an object pair that has been in contact, looses contact.
41-
Since `distanceOrg` values in `contactDict` are in the order of ``10^{-3} .. 10^{-6} ~m``
60+
Since penetration depth is in the order of ``10^{-3} .. 10^{-6} ~m``
4261
and `eps(Float64)` ``\approx 10^{-16}``, a hysteresis epsilon must be larger as ``10^{-19} .. 10^{-24}``.
43-
Actually, a hysteresis epsilon of ``10^{-16}`` is used. Note, that in `contactDict` all `distanceOrg` values
44-
are zero or negative at event restart. It is therefore guaranteed that ``z_1 \le - 10^{-16}`` at event restart
62+
Actually, a hysteresis epsilon of ``10^{-13}`` is used. Note, it is guaranteed that
63+
``z_1 \le`` `- contact_eps` at event restart
4564
(when no contact pair is present, a dummy value is used).
4665

47-
* ``z_2(t)``: The minimum `distanceOrg` of all object pairs that are **not** in `contactDict`.
66+
* ``z_2(t)``: The minimum of `distanceOrg+3*contact_eps` of all object pairs that are **not** in `contactDict`.
4867
``z_2`` monitors if two objects that have not been in contact to each other and start to penetrate each other.
49-
At an event restart (including the start of the integration after the initialization), it is guaranteed
50-
that all `distanceOrg` not in `contactDict` are positive (otherwise, they would be
51-
in `contactDict`). Therefore, it is guaranteed that ``z_2 > 0`` at event restart (without a hysteresis epsilon).
52-
53-
The distance between two objects is only computed up to a certain precision. Therefore, `distanceOrg` might be
54-
nearly zero but still positive, if two objects are placed initially in touching position.
55-
For this reason the actually used distance is `distanceOrg - 1e-16`,
56-
in order that the probability is larger that objects initially in touching position
57-
are treated to be in contact initially.
68+
At an event restart (including the start of the integration after the initialization),
69+
all contact pairs not in `contactDict` have ``z_2 > 0`` (otherwise, they would be
70+
in `contactDict`). Therefore, it is guaranteed that ``z_2 > 0`` at event restart.
5871

5972
To summarize, the following equations are actually used:
6073

6174
```
62-
distance = distanceOrg - 1e-16
63-
contact = isEvent ? distance <= 0 : <contact from previous event>
64-
z[1] = max(<distance that has contact=true>) - 1e-16
65-
z[2] = min(<distance that has contact=false>)
75+
contact = isEvent ? distanceOrg < -2*contact_eps : <contact from last event>
76+
distance = contact ? distanceOrg + contact_eps : distanceOrg + 3*contact_eps
77+
z[1] = max(<distanceOrg+contact_eps that has contact=true>)
78+
z[2] = min(<distanceOrg+3*contact_eps that has contact=false>)
6679
```
6780

6881
Typically, the `simulate!(..., log=true, ...)` produces the following output (here for a sphere boucing on ground):
6982

7083
```
71-
State event (zero-crossing) at time = 2.3924142776549155 s (due to z[2])
72-
distance(ground,sphere) = -3.815173674537923e-18 became <= 0
73-
contact normal = [0.0, 1.0, 0.0], contact position = [0.0, -0.1, 0.0], c_res = 1.1e11 d_res = 47.5
74-
restart = Restart
75-
76-
State event (zero-crossing) at time = 2.3924413895015424 s (due to z[1])
77-
distance(ground,sphere) = 1.0362163563391329e-16 became > 0
78-
restart = Restart
84+
State event (zero-crossing) at time = 2.3850903177319287 s (due to z[2])
85+
distance(ground,sphere) = -2.0001204970840633e-13 became <= 0
86+
contact normal = [0.0, 1.0, 0.0], contact position = [0.0, -0.1, 0.0], c_res = 1.1e11 d_res = 22.0
87+
restart = Restart
88+
89+
State event (zero-crossing) at time = 2.3851135168512942 s (due to z[1])
90+
distance(ground,sphere) = 3.1358549064074168e-18 became > 0
91+
restart = Restart
7992
```
8093

8194
Whenever the integrator requires the values of the zero crossing functions, the values of

src/contactDetection/ContactDetectionMPR/handler.jl

Lines changed: 11 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -76,11 +76,11 @@ function selectContactPairs!(sim, scene::Composition.Scene{F}, ch::Composition.C
7676
simh.z[scene.zStartIndex] = F(-42.0)
7777
else
7878
(pair, key) = findmax(ch.contactDict)
79-
simh.z[scene.zStartIndex] = pair.distanceWithHysteresis - ch.contact_eps
79+
simh.z[scene.zStartIndex] = pair.distanceWithHysteresis
8080
end
8181
# z[2] ... zero crossing function from no contact to contact
8282
# min. distance of inactive contact pairs
83-
simh.z[1+scene.zStartIndex] = ch.noContactMinVal
83+
simh.z[1+scene.zStartIndex] = ch.noContactMinVal #+ 2*ch.contact_eps
8484
end
8585
ch.distanceComputed = true
8686
return nothing
@@ -189,9 +189,15 @@ function storeDistancesForSolver!(world::Composition.Object3D{F}, pairID::Compos
189189
# println("AABB not overlapping")
190190
end
191191

192-
distanceWithHysteresis = distanceOrg - ch.contact_eps
193-
contact = hasEvent ? distanceWithHysteresis <= F(0) : hasContact
194-
192+
#distanceWithHysteresis = distanceOrg + ch.contact_eps
193+
#contact = hasEvent ? distanceWithHysteresis <= F(0) : hasContact
194+
195+
contact = hasEvent ? distanceOrg < -2*ch.contact_eps : hasContact
196+
distanceWithHysteresis = contact ? distanceOrg + ch.contact_eps : distanceOrg + 3*ch.contact_eps
197+
198+
#contact = hasEvent ? distanceOrg < -ch.contact_eps : hasContact
199+
#distanceWithHysteresis = contact ? distanceOrg : distanceOrg + 2*ch.contact_eps
200+
195201
if hasEvent
196202
# At event instant
197203
pushCollisionPair!(ch, contact, distanceWithHysteresis, pairID,

test/Collision/BouncingSphere2.jl

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,6 @@ BouncingSphere = Model(
1616
defaultFrameLength=0.2,
1717
visualizeBoundingBox = true,
1818
enableContactDetection=true,
19-
maximumContactDamping=1e3,
2019
visualizeContactPoints=false)),
2120
worldFrame = Object3D(parent=:world, feature=Visual(shape=CoordinateSystem(length=0.5))),
2221
ground = Object3D(parent=:world,

test/Robot/YouBotWithSphere.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -429,8 +429,8 @@ youbot = @instantiateModel(youbotModel, unitless=true, logCode=false, log=false)
429429

430430
stopTime = 5.0
431431
tolerance = 1e-6
432-
requiredFinalStates = [0.3872493507580361, -0.00016295756396963518, -0.3452201016630769, 0.15545045709830504, -3.11821309683207e-5, 1.9330433770085904e-7, 1.5711492875200297, -0.000498564591859475, -21.347006443600556, 0.0029361171926609074, 0.0008107479063527912, -6.245138200548544, -9.383806859286387e-8, 9.390066582743576e-8, -1.8339186184170656e-6, 1.8341287555337681e-6, -3.5545167503965664e-6, 3.5549976354276426e-6, -2.181811006068393e-6, 2.1821737128984775e-6, 8.488411980182589e-9, -8.489823786317174e-9, -0.0007168100048350477, 0.2276185370576126, -0.018209579211904647, -0.008528817908688886, 2.9526174065185742e-5]
433-
simulate!(youbot, stopTime=stopTime, tolerance=tolerance, requiredFinalStates_atol=0.001, log=true, logStates=false, requiredFinalStates=requiredFinalStates)
432+
requiredFinalStates = [0.38481650062588163, -0.00016295433833872128, -0.3452201051978677, 0.1536726214508672, -3.110346807218785e-5, 2.172037533153316e-7, 1.571174013759277, -0.00046190239792032596, -21.251971315603686, 0.0028911425459974627, 0.0007934643317099368, -6.173710969663145, -9.385568936457577e-8, 9.391828914937129e-8, -1.8339879373678193e-6, 1.8341980832035228e-6, -3.554672609027041e-6, 3.5551535140628727e-6, -2.181898901956027e-6, 2.1822616237942205e-6, 8.488749686546869e-9, -8.490161551079439e-9, -0.0007169446781250749, 0.22761800725155104, -0.018210342819184065, -0.008529123645210713, 2.9527348738343718e-5]
433+
simulate!(youbot, stopTime=stopTime, tolerance=tolerance, requiredFinalStates_atol=0.01, log=true, logStates=false, requiredFinalStates=requiredFinalStates)
434434

435435
@usingModiaPlot
436436
plot(youbot, ["free.rot"], figure=1)

0 commit comments

Comments
 (0)