Skip to content

Commit 293ae8f

Browse files
Merge pull request #80 from ModiaSim/gh_revertContactHysteresisModification
Revert contact hysteresis modification of #76
2 parents 16f346b + 4d5d129 commit 293ae8f

File tree

3 files changed

+36
-55
lines changed

3 files changed

+36
-55
lines changed

docs/src/internal/ContactDetection.md

Lines changed: 29 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -23,72 +23,59 @@ 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. 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.
26+
If `distanceOrg < 0`, then the two objects are penetrating each other.
4527

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

5133
At an **event instant** (including initialization), dictionary `contactDict` is emptied
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.
34+
and all object pairs are stored newly in `contactDict` that have `distanceOrg <= 0`, so are
35+
either touching or penetrating.
5536

5637
During **continuous integration**, two zero crossing functions are used:
5738

58-
* ``z_1(t)``: The maximum of `distanceOrg+contact_eps` of all object pairs that are in `contactDict`.
39+
* ``z_1(t)``: The maximum of `distanceOrg` ``- 10^{-16}`` of all object pairs that are in `contactDict`.
5940
``z_1`` monitors if an object pair that has been in contact, looses contact.
60-
Since penetration depth is in the order of ``10^{-3} .. 10^{-6} ~m``
41+
Since `distanceOrg` values in `contactDict` are in the order of ``10^{-3} .. 10^{-6} ~m``
6142
and `eps(Float64)` ``\approx 10^{-16}``, a hysteresis epsilon must be larger as ``10^{-19} .. 10^{-24}``.
62-
Actually, a hysteresis epsilon of ``10^{-13}`` is used. Note, it is guaranteed that
63-
``z_1 \le`` `- contact_eps` at event restart
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
6445
(when no contact pair is present, a dummy value is used).
6546

66-
* ``z_2(t)``: The minimum of `distanceOrg+3*contact_eps` of all object pairs that are **not** in `contactDict`.
47+
* ``z_2(t)``: The minimum `distanceOrg` of all object pairs that are **not** in `contactDict`.
6748
``z_2`` monitors if two objects that have not been in contact to each other and start to penetrate each other.
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.
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.
7158

7259
To summarize, the following equations are actually used:
7360

7461
```
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>)
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>)
7966
```
8067

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

8370
```
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
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
9279
```
9380

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

src/contactDetection/ContactDetectionMPR/handler.jl

Lines changed: 5 additions & 11 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
79+
simh.z[scene.zStartIndex] = pair.distanceWithHysteresis - ch.contact_eps
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 #+ 2*ch.contact_eps
83+
simh.z[1+scene.zStartIndex] = ch.noContactMinVal
8484
end
8585
ch.distanceComputed = true
8686
return nothing
@@ -189,15 +189,9 @@ 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-
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-
192+
distanceWithHysteresis = distanceOrg - ch.contact_eps
193+
contact = hasEvent ? distanceWithHysteresis <= F(0) : hasContact
194+
201195
if hasEvent
202196
# At event instant
203197
pushCollisionPair!(ch, contact, distanceWithHysteresis, pairID,

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.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)
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)
434434

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

0 commit comments

Comments
 (0)