Skip to content

Commit bef2493

Browse files
Merge pull request #109 from ModiaSim/an_mprImprovements
An mpr improvements
2 parents a36390e + 4260d5a commit bef2493

File tree

9 files changed

+266
-72
lines changed

9 files changed

+266
-72
lines changed

src/Composition/scene.jl

Lines changed: 10 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -239,7 +239,7 @@ struct SceneOptions{F <: Modia3D.VarFloatType}
239239
gap = 0.001,
240240
enableVisualization = true,
241241
animationFile = nothing,
242-
provideAnimationHistory = false,
242+
provideAnimationHistory = false,
243243
visualizeFrames = false,
244244
visualizeBoundingBox = false,
245245
visualizeContactPoints = false,
@@ -323,6 +323,7 @@ Defines global properties of the system, such as the gravity field. Exactly one
323323
| `elasticContactReductionFactor` | 1.0 (> 0.0, <= 1.0) |
324324
| `maximumContactDamping` | 2000.0 |
325325
| `mprTolerance` | 1.0e-20 |
326+
| `mprIterMax` | 120 |
326327
| `visualizeFrames` | false |
327328
| `visualizeBoundingBox` | false |
328329
| `visualizeContactPoints` | false |
@@ -364,6 +365,8 @@ Defines global properties of the system, such as the gravity field. Exactly one
364365
- `mprTolerance::1.0e-20`: Local tolerance used for terminating the mpr algorithm (that computes the distances between shapes). Changing this value might improve speed.
365366
For integrators with step-size control, a value is needed that is much smaller as the relative tolerance used for the integration.
366367
368+
- `mprIterMax::120`: Local maximum amount of iterations used for mpr algorithm. If more iterations are needed a message is printed.
369+
367370
- `visualizeFrames::Bool`: = true, to visualize the coordinate system of every [Object3D](@ref) that is not explicitly switched off.
368371
369372
- `visualizeBoundingBox::Bool`: Flag enabled for visualizing Axis Aligned Bounding Box (AABB) for all solid shapes allowed to collide
@@ -427,11 +430,11 @@ mutable struct Scene{F <: Modia3D.VarFloatType} <: Modia3D.AbstractScene
427430
AABB::Vector{Vector{Basics.BoundingBox{F}}} # Bounding boxes of elements that can collide
428431
zStartIndex::Int # start index of collision zero crossing functions
429432
forceElements::Vector{Modia3D.AbstractForceElement}
430-
provideAnimationData::Bool # = true, if animation data shall be provided
433+
provideAnimationData::Bool # = true, if animation data shall be provided
431434
exportAnimation::Bool # animation file export is enabled
432435
animation::Vector{animationStep} # animation data of visible Object3Ds
433436
outputCounter::Int64 # animation/visualization output step counter
434-
437+
435438
# Data specific to a particular joint type
436439
revolute::Vector{Revolute{F}}
437440
prismatic::Vector{Prismatic{F}}
@@ -442,6 +445,7 @@ mutable struct Scene{F <: Modia3D.VarFloatType} <: Modia3D.AbstractScene
442445
useOptimizedStructure = true,
443446
enableContactDetection = true,
444447
mprTolerance = 1.0e-20,
448+
mprIterMax = 120,
445449
elasticContactReductionFactor = F(1.0),
446450
maximumContactDamping = F(2000),
447451
gap = 0.001,
@@ -465,7 +469,7 @@ mutable struct Scene{F <: Modia3D.VarFloatType} <: Modia3D.AbstractScene
465469

466470
sceneOptions = SceneOptions{F}(gravityField = gravityField,
467471
useOptimizedStructure = useOptimizedStructure,
468-
contactDetection = ContactDetectionMPR_handler{Modia3D.MPRFloatType, F}(tol_rel = mprTolerance),
472+
contactDetection = ContactDetectionMPR_handler{Modia3D.MPRFloatType, F}(tol_rel = mprTolerance, niter_max = mprIterMax),
469473
nVisualContSupPoints = nVisualContSupPoints,
470474
gap = gap,
471475
enableContactDetection = enableContactDetection,
@@ -476,7 +480,7 @@ mutable struct Scene{F <: Modia3D.VarFloatType} <: Modia3D.AbstractScene
476480
defaultFrameLength = defaultFrameLength,
477481
enableVisualization = enableVisualization,
478482
animationFile = animationFile,
479-
provideAnimationHistory = provideAnimationHistory,
483+
provideAnimationHistory = provideAnimationHistory,
480484
visualizeFrames = visualizeFrames,
481485
visualizeBoundingBox = visualizeBoundingBox,
482486
visualizeContactPoints = visualizeContactPoints,
@@ -524,7 +528,7 @@ mutable struct Scene{F <: Modia3D.VarFloatType} <: Modia3D.AbstractScene
524528
Vector{Vector{Basics.BoundingBox{F}}}[],
525529
1,
526530
Vector{Modia3D.AbstractForceElement}[],
527-
provideAnimationData,
531+
provideAnimationData,
528532
exportAnimation,
529533
Vector{animationStep}[],
530534
0,

src/contactDetection/ContactDetectionMPR/ContactDetectionMPR_handler.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -106,7 +106,7 @@ mutable struct ContactDetectionMPR_handler{T,F} <: Modia3D.AbstractContactDetect
106106
defaultContactSphereDiameter::Float64
107107

108108
function ContactDetectionMPR_handler{T,F}(;tol_rel = 1.0e-20,
109-
niter_max = 100) where {T,F}
109+
niter_max = 120) where {T,F}
110110
@assert(tol_rel > 0.0)
111111
@assert(niter_max > 0)
112112
contact_eps = F(max(1e-13, 10*tol_rel))

src/contactDetection/ContactDetectionMPR/mpr.jl

Lines changed: 47 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ struct SupportPoint{T}
1515
end
1616

1717

18-
function getSupportPoint(shapeA::Modia3D.Composition.Object3D, shapeB::Composition.Object3D, n::SVector{3,T}; scale::T=T(1.0)) where {T}
18+
function getSupportPoint(shapeA::Modia3D.Composition.Object3D{F}, shapeB::Composition.Object3D{F}, n::SVector{3,T}; scale::T=T(1.0)) where {T,F}
1919
a = Modia3D.supportPoint(shapeA, n)
2020
b = Modia3D.supportPoint(shapeB, -n)
2121
return SupportPoint{T}((a-b).*scale,n,a,b)
@@ -66,15 +66,15 @@ end
6666

6767
# checks if centers of shapeA and shapeB are overlapping
6868
# belongs to construction of r0
69-
function checkCentersOfShapesOverlapp(r0::SupportPoint{T}, shapeA::Composition.Object3D, shapeB::Composition.Object3D) where {T}
69+
function checkCentersOfShapesOverlapp(r0::SupportPoint{T}, shapeA::Composition.Object3D{F}, shapeB::Composition.Object3D{F}) where {T,F}
7070
if norm(r0.p) <= Modia3D.nepsType(T)
7171
error("MPR: Too large penetration (prerequisite of MPR violated). Centers are overlapping. Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
7272
end
7373
end
7474

7575

76-
function checkIfShapesArePlanar(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,n2::SVector{3,T},
77-
shapeA::Composition.Object3D,shapeB::Composition.Object3D) where {T}
76+
function checkIfShapesArePlanar(r0::SupportPoint{T},r1::SupportPoint{T},r2::SupportPoint{T},n2::SVector{3,T},
77+
shapeA::Composition.Object3D{F},shapeB::Composition.Object3D{F}) where {T,F}
7878
r3 = SupportPoint{T}
7979

8080
# r3 is in the direction of plane normal that contains triangle r0-r1-r2
@@ -86,14 +86,13 @@ function checkIfShapesArePlanar(r0::SupportPoint,r1::SupportPoint,r2::SupportPoi
8686
# because we are still interested in distances if shapes are not intersecting
8787
n2 = -n2
8888
r2 = getSupportPoint(shapeA, shapeB, n2)
89-
if abs(dot((r2.p-r1.p),n2)) <= neps
89+
n3 = cross(r1.p-r0.p, r2.p-r0.p) # new normal to the triangle plane (r0-r1-r2_new)
90+
if norm(n3) <= neps
9091
# Shape is purely planar. Computing the shortest distance for a planar shape
9192
# requires an MPR 2D algorithm (using lines instead of triangles as portals).
9293
# However, this is not implemented and therefore the shortest distance cannot be computed
93-
error("MPR: Shapes are planar and MPR2D is not supported. abs(dot((r2.p-r1.p),n2)). Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
94+
error("MPR: Shapes are planar and MPR2D is not supported. norm(cross(r1.p-r0.p, r2.p-r0.p)). Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
9495
end
95-
# new normal to the triangle plane (r0-r1-r2_new)
96-
n3 = cross(r1.p-r0.p, r2.p-r0.p) # |n3| > 0 guaranteed, due to construction
9796
end
9897

9998
if dot(n3,r0.p) >= 0.0
@@ -106,11 +105,11 @@ function checkIfShapesArePlanar(r0::SupportPoint,r1::SupportPoint,r2::SupportPoi
106105
if norm(n3b) <= neps
107106
# change search direction for r3
108107
r3 = getSupportPoint(shapeA, shapeB, -r3.n)
109-
if abs(dot((r3.p-r1.p),r3.n)) <= neps
108+
if norm(cross(r2.p-r1.p, r3.p-r1.p)) <= neps
110109
# Shape is purely planar. Computing the shortest distance for a planar shape
111110
# requires an MPR 2D algorithm (using lines instead of triangles as portals).
112111
# However, this is not implemented and therefore the shortest distance cannot be computed
113-
error("MPR: Shapes are planar and MPR2D is not supported. r1, r2, r3 are on the same ray. abs(dot((r3.p-r1.p),r3.n)) <= neps. Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
112+
error("MPR: Shapes are planar and MPR2D is not supported. r1, r2, r3 are on the same ray. norm(cross(r2.p-r1.p, r3.p-r1.p)) <= neps. Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
114113
end
115114
end
116115

@@ -122,10 +121,10 @@ end
122121
# loop around to "ensure" the tetrahedron r0,r1,r2 and r3 encloses the origin
123122
# stimmt so nicht wirklich, muss ich nochmal nachlesen!!!
124123
# Der Ursprung muss nicht enthalten sein!!!
125-
function tetrahedronEncloseOrigin(r0::SupportPoint, r1::SupportPoint,
126-
r2::SupportPoint, r3::SupportPoint,
124+
function tetrahedronEncloseOrigin(r0::SupportPoint{T}, r1::SupportPoint{T},
125+
r2::SupportPoint{T}, r3::SupportPoint{T},
127126
niter_max::Int64,
128-
shapeA::Composition.Object3D, shapeB::Composition.Object3D, scale::T) where {T}
127+
shapeA::Composition.Object3D{F}, shapeB::Composition.Object3D{F}, scale::T) where {T,F}
129128
r1org = r1
130129
r2org = r2
131130
r3org = r3
@@ -134,13 +133,13 @@ function tetrahedronEncloseOrigin(r0::SupportPoint, r1::SupportPoint,
134133
success = false
135134
for i in 1:niter_max
136135
aux = cross(r1.p-r0.p,r3.p-r0.p)
137-
if dot(aux,r0.n) < -neps
136+
if dot(aux,r0.p) > 0.0
138137
r2 = r3
139138
r3 = getSupportPoint(shapeA,shapeB,Basics.normalizeVector(aux), scale=scale)
140139
continue
141140
end
142141
aux = cross(r3.p-r0.p,r2.p-r0.p)
143-
if dot(aux,r0.n) < -neps
142+
if dot(aux,r0.p) > 0.0
144143
r1 = r3
145144
r3 = getSupportPoint(shapeA,shapeB,Basics.normalizeVector(aux), scale=scale)
146145
continue
@@ -149,33 +148,28 @@ function tetrahedronEncloseOrigin(r0::SupportPoint, r1::SupportPoint,
149148
break
150149
end
151150
if success != true
152-
if niter_max < 100
153-
@warn("MPR (phase 2): Max. number of iterations (= $niter_max) is reached. niter_max increased locally by 10 and phase 2 is rerun. Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
154-
tetrahedronEncloseOrigin(r0, r1org, r2org, r3org, niter_max + 10, shapeA, shapeB, scale)
155-
else
156-
error("MPR (phase 2): Max. number of iterations (= $niter_max) is reached and $niter_max > 100, look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
157-
end
151+
error("MPR (phase 2): Max. number of iterations (mprIterMax = $niter_max) is reached. Please, increase mprIterMax. Look at shapes: $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
158152
end
159153
return (r1, r2, r3)
160154
end
161155

162156

163157
########### Phase 3, Minkowski Portal Refinement ###################
164158
# construction of r4
165-
function constructR4(r0::SupportPoint,r1::SupportPoint,r2::SupportPoint,r3::SupportPoint,
166-
shapeA::Composition.Object3D,shapeB::Composition.Object3D, scale::T) where {T}
159+
function constructR4(r0::SupportPoint{T},r1::SupportPoint{T},r2::SupportPoint{T},r3::SupportPoint{T},
160+
shapeA::Composition.Object3D{F},shapeB::Composition.Object3D{F}, scale::T) where {T,F}
167161
r4 = SupportPoint{T}
168162
n4 = cross(r2.p-r1.p, r3.p-r1.p)
169163
neps = Modia3D.nepsType(T)
170164
if norm(n4) <= neps
171165
r3 = getSupportPoint(shapeA, shapeB, -r3.n, scale=scale) # change search direction
172-
if abs(dot((r3.p-r1.p),r3.n)) <= neps
166+
n4 = cross(r2.p-r1.p, r3.p-r1.p)
167+
if norm(n4) <= neps
173168
# Shape is purely planar. Computing the shortest distance for a planar shape
174169
# requires an MPR 2D algorithm (using lines instead of triangles as portals).
175170
# However, this is not implemented and therefore the shortest distance cannot be computed
176-
error("MPR: Shapes are planar and MPR2D is not supported. abs(dot((r3.p-r1.p),r3.n)). Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
171+
error("MPR: Shapes are planar and MPR2D is not supported. norm(n4). Look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
177172
end
178-
n4 = cross(r2.p-r1.p, r3.p-r1.p) # |n4| > 0 guaranteed, due to construction
179173
end
180174
if dot(n4,r0.p) >= 0.0
181175
n4 = -n4
@@ -196,7 +190,7 @@ isNextPortal(r0::SupportPoint, r1::SupportPoint, r2::SupportPoint, r4::SupportPo
196190
# r1r2r4
197191
# r2r3r4
198192
# r3r1r4
199-
function createBabyTetrahedrons(r0::SupportPoint, r1::SupportPoint, r2::SupportPoint, r3::SupportPoint, r4::SupportPoint)
193+
function createBabyTetrahedrons(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoint{T}, r3::SupportPoint{T}, r4::SupportPoint{T}) where {T}
200194
nextPortal = true
201195
if isNextPortal(r0,r1,r2,r4)
202196
r3 = r4
@@ -224,7 +218,7 @@ function skalarization(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoin
224218
return (r0, r1, r2, r3, scale)
225219
end
226220

227-
function finalTC2(r1::SupportPoint, r2::SupportPoint, r3::SupportPoint, r4::SupportPoint)
221+
function finalTC2(r1::SupportPoint{T}, r2::SupportPoint{T}, r3::SupportPoint{T}, r4::SupportPoint{T}) where {T}
228222
#println("TC 2")
229223
#if !analyzeFinalPortal(r1.p, r2.p, r3.p, r4.p)
230224
# error("shapeA = ", shapeA, " shapeB = ", shapeB)
@@ -239,9 +233,6 @@ function finalTC3(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoint{T},
239233
#println("TC 3")
240234

241235
#doesRayIntersectPortal(r1.p,r2.p,r3.p, r4.p)
242-
#println("r1.p = ", r1.p , " r2.p = ", r2.p ," r3.p = ", r3.p)
243-
#println("r4.p = ", r4.p)
244-
#println(" ")
245236
#if !analyzeFinalPortal(r1.p, r2.p, r3.p, r4.p)
246237
# error("shapeA = ", shapeA, " shapeB = ", shapeB)
247238
#end
@@ -253,7 +244,20 @@ function finalTC3(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoint{T},
253244
end
254245

255246

256-
function phase3(r0::SupportPoint, r1::SupportPoint, r2::SupportPoint, r3::SupportPoint, niter_max::Int64, tol_rel::T, shapeA::Composition.Object3D, shapeB::Composition.Object3D, scale::T) where {T}
247+
function terminateMPR(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoint{T},
248+
r3::SupportPoint{T}, r4::SupportPoint{T}, isTC2::Bool, isTC3::Bool)::Tuple{T, SupportPoint{T}, SupportPoint{T}, SupportPoint{T}, SupportPoint{T}} where {T}
249+
if isTC2
250+
(distance,r1,r2,r3,r4) = finalTC2(r1, r2, r3, r4)
251+
return distance, r1, r2, r3, r4
252+
end
253+
if isTC3
254+
(distance,r1,r2,r3,r4) = finalTC3(r0, r1, r2, r3, r4)
255+
return distance, r1, r2, r3, r4
256+
end
257+
end
258+
259+
260+
function phase3(r0::SupportPoint{T}, r1::SupportPoint{T}, r2::SupportPoint{T}, r3::SupportPoint{T}, niter_max::Int64, tol_rel::T, shapeA::Composition.Object3D{F}, shapeB::Composition.Object3D{F}, scale::T) where {T,F}
257261
r1org = r1
258262
r2org = r2
259263
r3org = r3
@@ -277,13 +281,11 @@ function phase3(r0::SupportPoint, r1::SupportPoint, r2::SupportPoint, r3::Suppor
277281
TC3 = abs(dot(r4.p-r1.p, r4.n)) # TC3
278282
## TERMINATION CONDITION 2 ##
279283
if TC2 < tol_rel
280-
(distance,r1,r2,r3,r4) = finalTC2(r1,r2,r3,r4)
281-
return distance, r1, r2, r3, r4
284+
return finalTC2(r1,r2,r3,r4)
282285

283286
## TERMINATION CONDITION 3 ##
284287
elseif TC3 < tol_rel
285-
(distance,r1,r2,r3,r4) = finalTC3(r0, r1, r2, r3, r4)
286-
return distance, r1, r2, r3, r4
288+
return finalTC3(r0, r1, r2, r3, r4)
287289
else
288290
if TC2 < new_tol
289291
new_tol = TC2
@@ -311,32 +313,11 @@ function phase3(r0::SupportPoint, r1::SupportPoint, r2::SupportPoint, r3::Suppor
311313

312314
if !nextPortal # createBabyTetrahedrons failed
313315
@warn("MPR (phase 3): Numerical issues with distance computation between $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)). tol_rel increased locally for this computation to $new_tol.")
314-
if isTC2
315-
(distance,r1,r2,r3,r4) = finalTC2(r1_new,r2_new,r3_new,r4_new)
316-
return distance, r1, r2, r3, r4
317-
end
318-
if isTC3
319-
(distance,r1,r2,r3,r4) = finalTC3(r0, r1_new, r2_new, r3_new, r4_new)
320-
return distance, r1, r2, r3, r4
321-
end
316+
return terminateMPR(r0, r1_new, r2_new, r3_new, r4_new, isTC2, isTC3)
322317
end
323318
end
324-
if niter_max < 100
325-
@warn("MPR (phase 3): Numerical issues with distance computation between $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)). Max. number of iterations (= $niter_max) is reached. niter_max increased locally by 10 and phase 3 is rerun.")
326-
phase3(r0, r1org, r2org, r3org, niter_max + 10, tol_rel, shapeA, shapeB, scale)
327-
else
328-
@warn("MPR (phase 3): Max. number of iterations (= $niter_max) is reached and $niter_max > 100, look at $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)). tol_rel increased locally for this computation to $new_tol.")
329-
if isTC2
330-
(distance,r1,r2,r3,r4) = finalTC2(r1_new,r2_new,r3_new,r4_new)
331-
return distance, r1, r2, r3, r4
332-
end
333-
if isTC3
334-
(distance,r1,r2,r3,r4) = finalTC3(r0, r1_new, r2_new, r3_new, r4_new)
335-
return distance, r1, r2, r3, r4
336-
end
337-
end
338-
@error("passiert das?!?!?")
339-
return distance, r1, r2, r3, r4 # needed for a unique return type
319+
@warn("MPR (phase 3): Max. number of iterations (mprIterMax = $niter_max) is reached. Please, increase mprIterMax. tol_rel increased locally for this computation to $new_tol. Look at shapes $(Modia3D.fullName(shapeA)) and $(Modia3D.fullName(shapeB)).")
320+
return terminateMPR(r0, r1_new, r2_new, r3_new, r4_new, isTC2, isTC3)
340321
end
341322

342323
# MPR - Minkowski Portal Refinement algorithm
@@ -363,7 +344,7 @@ end
363344
# Termination Condition 2
364345
# Termination Condition 3
365346
# Phase 3.3: construct baby tetrahedrons with r1,r2,r3,r4 and create a new portal
366-
function mprGeneral(ch::Composition.ContactDetectionMPR_handler{T,F}, shapeA::Composition.Object3D, shapeB::Modia3D.Composition.Object3D) where {T,F}
347+
function mprGeneral(ch::Composition.ContactDetectionMPR_handler{T,F}, shapeA::Composition.Object3D{F}, shapeB::Modia3D.Composition.Object3D{F}) where {T,F}
367348
tol_rel = ch.tol_rel
368349
niter_max = ch.niter_max
369350
neps = Modia3D.nepsType(T)
@@ -379,6 +360,11 @@ function mprGeneral(ch::Composition.ContactDetectionMPR_handler{T,F}, shapeA::Co
379360
# the direction of the origin ray r0 is -r0.p
380361
centroidA = getCentroid(shapeA)
381362
centroidB = getCentroid(shapeB)
363+
if isnan(centroidA[1]) || isnan(centroidB[1]) ||
364+
isnan(centroidA[2]) || isnan(centroidB[2]) ||
365+
isnan(centroidA[3]) || isnan(centroidB[3])
366+
error("MPR: One of the absolute position or translation is NaN. Look at $(Modia3D.fullName(shapeA)): r_abs = $(shapeA.r_abs), R_abs = $(shapeA.R_abs) and $(Modia3D.fullName(shapeB)): r_abs = $(shapeB.r_abs), R_abs = $(shapeB.R_abs).")
367+
end
382368
r0 = SupportPoint{T}(centroidA-centroidB, -(centroidA-centroidB), SVector{3,T}(0.0,0.0,0.0), SVector{3,T}(0.0,0.0,0.0))
383369
# check if centers of shapes are overlapping
384370
checkCentersOfShapesOverlapp(r0, shapeA, shapeB)

0 commit comments

Comments
 (0)