[PPT] - Stable, Robust, and Versatile Multibody Dynamics Animation Kenny PowerPoint Presentation

SLIDE 1

Stable, Robust, and Versatile Multibody Dynamics Animation

Kenny Erleben Department of Computer Science University of Copenhagen

c

K. Erleben

– p. 1/25

SLIDE 2

The Problem

Large scale (today ≥ 1000 objects). Dense random stacking Dense structured stacking Sparse Stacking Dense ⇒ more difficult. Random looks OK if top-most objects OK. More sparse ⇒ less dependent constraints ≡ faster convergence. Notice: Dense ⇒ more dependent constraints. Stable ⇒ Handling Physical Unstable and Ill-posed. Robust ⇒ Handling Degenerate and Faulty Cases. The Attack of The Malicious User!

c

K. Erleben

– p. 2/25

SLIDE 3

My Solution

Ingredients to achieve high-performance: Velocity-based complementarity formulations. Leap-frog like time-stepping scheme. Iterative LCP solver. Error-correction by projection. Shock-propagation, followed by a smoothing correction phase.

c

K. Erleben

– p. 3/25

SLIDE 4

Complementarity Formulation

From physics Normal force repulsive (λ1 ≥ 0). Normal force zero at separation (λ1 = 0). If bodies are moving apart w1 = Jrow1 u > 0 ⇒ λ1 = 0,

(1)

If bodies are resting λ1 > 0 ⇒ w1 = Jrow1 u = 0.

(2) w

1

=0 w

1

>0 f

1

>0 f

1

=0

Either λ1 is zero and Jrow1 u is non-zero or vice-versa. w1 = Jrow1 u ≥ 0 compl. λ1 ≥ 0

(3)

c

K. Erleben

– p. 4/25

SLIDE 5

Complementarity Formulation

Coulomb’s Friction Model.

J

row1u>0

f

1

f

2

friction cone

Dynamic friction: w2 = Jrow2 u > 0 ⇒ λ2 = −µλ1,

(4)

w2 = Jrow2 u < 0 ⇒ λ2 = µλ1,

(5)

µ is the coefficient of friction.

J

row1u=0

f

1

f

2

=?

friction cone

Static friction: w2 = Jrow2 u = 0 ⇒ λ2 <| µλ1 |,

(6)

Similar constraints for w3 = Jrow3 u and λ3.

c

K. Erleben

– p. 5/25

SLIDE 6

Complementarity Formulation

Combine discritization of M˙

u =

fext − J T λ with w = J u

w = JM−1J T

| {z }

A

∆t λ + J

u t + ∆tM−1

fext

|

{z }

b
w = A

λ + b

(7)

Notice ∆t is moved into the Lagrange multiplier vector.

w linear in

λ. Further λi = λloi ⇒ wi ≥ 0

(8)

λi = λhii ⇒ wi ≤ 0

(9)

λloi < λi <

λhii ⇒ wi = 0

(10)

A box linear complementarity problem (LCP) problem.

c

K. Erleben

– p. 6/25

SLIDE 7

The Time-Stepping Method

Dynamics explicit time-step scheme, dynamics (∆t)

(11)

given by

b = J
u t + ∆tM−1

fext

(12a)

A = JM−1J T

(12b)

λ = lcp
A,

b

,

(12c)

u t+1 =

u t + M−1J T λ + ∆tM−1 Fext,

(12d)

s t+1 =

s t + ∆tS u t+1.

(12e)

But how do we solve for the Lagrange multipliers? Brute Force...?

c

K. Erleben

– p. 7/25

SLIDE 8

Iterative LCP solver

The splitting A = L + D + U Iterative Gauss-Seidel A λ = − b

(13)

(L + D + U) λ = − b

(14)

λ k+1 = D−1
L

λ k − U λ k − b

(15)

Loop over all variables i ∈ [1..3K]

λ k+1

i

=

−

Pi−1

j=0 Li,j

λ k+1

j

−

Pn−1

j=i+1 Ui,j

λ k

j −

bi

Di,i

(16)

= −Lrowi λ k+1 − Urowi λ k − bi Di,i .

(17)

Superscript = iteration number.

c

K. Erleben

– p. 8/25

SLIDE 9

Iterative LCP solver

Upper and lower limits are updated if the i’th variable is a friction constraint, (r = i mod 3) = 0 ⇒ λloi = − λhii = µ λi−r,

(18)

projection step

λ k+1

i

= max

min
λloi,

λ k+1

i

,

λhii

.

(19)

Now Each row of the Jacobian have exactly 12 non-zero elements. M is 3-by-3 block diagonal matrix. Using a sparse matrix representation for M, J, and A, Linear time to computing A and b.

c

K. Erleben

– p. 9/25

SLIDE 10

LCP solver Results

Convergence Results Initial position. 10 iterations after 4 secs. 100 iterations after 4 secs. 10 iterations after 4 secs. 100 iterations after 4 secs. 10 iterations after 4 secs. 100 iterations after 4 secs.

c

K. Erleben

– p. 10/25

SLIDE 11

Convergence Rate

Convergence rate ≈ O(e−kn)

100 200 300 0.1 0.2 0.3 0.4 0.5

Convergence All

iterations θ(λi) θ(λi) = H(λi)T H(λi) / 2 Large mass Box Stack Ball Grid Wall Tower 100 200 300 10

−30

10

−20

10

−10

10

Log Convergence All

iterations log( θ(λi)) θ(λi) = H(λi)T H(λi) / 2 Large mass Box Stack Ball Grid Wall Tower

More structure (in sense of all-pair dependency) means increasing k, less structure means decreasing k. Value of θ has nothing to do with accuracy, when θ is “flat” we got enough accuracy.

c

K. Erleben

– p. 11/25

SLIDE 12

Error Correction

Penetrations are unavoidable Curved surfaces.

g

Undetected contacts.

g time t time t+dt g

c

K. Erleben

– p. 12/25

SLIDE 13

Error Correction

A first order world simulation: The equations of motion

F = m

a and

τ = d

L dt

(20)

is replaced by

F = m

v and

τ = I

ω.

(21)

Used for error correction, yields the scheme, correction ()

(22)

given by A = JM−1J T ,

(23a)

λ = lcp
A,

dpenetration

,

(23b)

s t+1 =

s t + SM−1J T λ,

(23c)

c

K. Erleben

– p. 13/25

SLIDE 14

Shock-Propagation

The general idea Velocity-based-Shock-Propagation(f,dt) collision detection at time t

dynamics(f*dt)

shock-propagration(

dynamics((1-f)*dt),
correction()

) collision detection at time t + dt

correction()

t = t + dt Velocity Correction Smoothing Position Correction

c

K. Erleben

– p. 14/25

SLIDE 15

Shock-Propagation

Stack Height Definition. Bodies are assigned a stack-height number. The stack height: ♯bodies on the closest path to a fixed body. Free bodies are assigned the maximum stack-height.

1 1 2 2 3 3 4 4 4 4

Notice: Stack-height is the path-cost of a breadth-first traversal starting at the fixed bodies.

c

K. Erleben

– p. 15/25

SLIDE 16

Shock-Propagation

Stack-Layer Definition: Stack layer i: Bodies with stack height i and i + 1, contact points between body-pairs with stack height i and i + 1, with stack height i + 1 and i + 1.

1 1 2 2 3 3 4 4 4 4

Stack Layer 1

Edges are given a stack layer number: If stack-heights are equal stack layer number is stack height minus one. Otherwise stack layer number is minimum.

c

K. Erleben

– p. 16/25

SLIDE 17

Shock-Propagation

The Shock-Propagation algorithm: Apply an algorithm sequentially to all stack layers. Treat stack layers in bottom-to-top order. Set bottom-most bodies to fixed before applying, afterwards unfix. shock-propagation(algorithm A) compute contact graph for each stack layer in bottom up order fixate bottom-most objects of layer apply algorithm A to layer un-fixate bottom-most objects of layer next layer

c

K. Erleben

– p. 17/25

SLIDE 18

Results

1000 balls. 2.0 secs. 4.0 secs. 6.0 secs. 8.0 secs. 250 cows. 1.0 secs. 2.0 secs. 3.0 secs. 4.0 secs. Roof and canon-ball. 2.0 secs. 4.0 secs. 6.0 secs. 8.0 secs.

c

K. Erleben

– p. 18/25

SLIDE 19

Results

Worst case frame-times: Engraving 0.18 secs Cow pile 1.2 secs Roof 0.25 secs Silo 1.1 secs Time-step = 0.01 secs. Pentium M, 1.7 GHz, 1GB RAM. Time-Stepping is linear in the number of contact points. Different slopes are a result of the stacking topology.

c

K. Erleben

– p. 19/25

SLIDE 20

Comparisons

Novodex Default

c

K. Erleben

– p. 20/25

SLIDE 21

Comparisons

Novodex 300 Iterations low separation distance

c

K. Erleben

– p. 21/25

SLIDE 22

Comparisons

Novodex 30 Iterations, Tweaked Mass, gravity,...

c

K. Erleben

– p. 22/25

SLIDE 23

Comparisons

Guendelman et. al. (Siggraph 2003)

c

K. Erleben

– p. 23/25

SLIDE 24

Comparisons

Maya Fatal error...

c

K. Erleben

– p. 24/25

SLIDE 25

Comparisons

200 400 600 800 1000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Iteration Time (secs)

Performance Comparison

Erleben Guendelman et. al. Novodex Default Novodex 30 Novodex 300

c

K. Erleben

– p. 25/25