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

Stable, Robust, and Versatile Multibody Dynamics Animation Kenny Erleben Department of Computer Science University of Copenhagen � K. Erleben c – p. 1/25

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! � K. Erleben c – p. 2/25

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. � K. Erleben c – p. 3/25

Complementarity Formulation From physics Normal force repulsive ( λ 1 ≥ 0 ). Normal force zero at separation ( λ 1 = 0 ). If bodies are moving apart w 1 = J row 1 � u > 0 ⇒ λ 1 = 0 , (1) If bodies are resting λ 1 > 0 ⇒ w 1 = J row 1 � u = 0 . (2) w� >0� 1� f� >0� 1� w� =0� 1� f� =0� 1� Either λ 1 is zero and J row 1 � u is non-zero or vice-versa. w 1 = J row 1 � u ≥ 0 λ 1 ≥ 0 compl. (3) � K. Erleben c – p. 4/25

Complementarity Formulation Coulomb’s Friction Model. Dynamic friction: friction cone� f� 1� w 2 = J row 2 � u > 0 ⇒ λ 2 = − µλ 1 , (4) J� row1� u>0� w 2 = J row 2 � u < 0 ⇒ λ 2 = µλ 1 , (5) µ is the coefficient of friction. f� 2� friction cone� Static friction: f� 1� f� =?� 2� w 2 = J row 2 � u = 0 ⇒ λ 2 < | µλ 1 | , (6) J� row1� u=0� Similar constraints for w 3 = J row 3 � u and λ 3 . � K. Erleben c – p. 5/25

Complementarity Formulation � � | {z } | {z } Combine discritization of M ˙ u = � f ext − J T � w = J � � λ with � u u t + ∆ t M − 1 � ∆ t� w = J M − 1 J T λ + J � � f ext A � b w = A � λ + � � b (7) Notice ∆ t is moved into the Lagrange multiplier vector. w linear in � λ . � Further λ i = � λ lo i ⇒ � w i ≥ 0 (8) λ i = � λ hi i ⇒ � w i ≤ 0 (9) � λ lo i < λ i < � λ hi i ⇒ � w i = 0 (10) A box linear complementarity problem (LCP) problem. � K. Erleben c – p. 6/25

The Time-Stepping Method But how do we solve for the Lagrange multipliers? Dynamics explicit time-step scheme, � � dynamics (∆ t ) (11) � � given by u t + ∆ t M − 1 � � b = J � f ext (12a) A = J M − 1 J T (12b) λ = lcp � A ,� b , (12c) u t +1 = � u t + M − 1 J T � λ + ∆ t M − 1 � � F ext , (12d) s t +1 = � s t + ∆ t S � u t +1 . � (12e) Brute Force...? � K. Erleben c – p. 7/25

Iterative LCP solver The splitting A = L + D + U � � Iterative Gauss-Seidel A � λ = − � b (13) ( L + D + U ) � λ = − � b (14) � � P i − 1 P n − 1 λ k +1 = D − 1 λ k − U � λ k − � � L � b (15) Loop over all variables i ∈ [1 .. 3 K ] j =0 L i,j � λ k +1 j = i +1 U i,j � j − � λ k − − b i j λ k +1 � = (16) i D i,i λ k +1 − U row i � λ k − � = − L row i � b i . (17) D i,i Superscript = iteration number. � K. Erleben c – p. 8/25

Iterative LCP solver � � � � Upper and lower limits are updated if the i ’th variable is a friction constraint, ( r = i mod 3) � = 0 ⇒ � λ lo i = − � λ hi i = µ� λ i − r , (18) projection step λ k +1 λ k +1 � λ lo i ,� � ,� = max min λ hi i . (19) i i 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 . � K. Erleben c – p. 9/25

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. � K. Erleben c – p. 10/25

Convergence Rate Convergence rate ≈ O ( e − kn ) Convergence All Log Convergence All 0 0.5 10 Large mass Box Stack 0.4 Ball Grid θ ( λ i ) = H( λ i ) T H( λ i ) / 2 Wall −10 10 Tower log( θ ( λ i )) 0.3 θ ( λ i ) θ ( λ i ) = H( λ i ) T H( λ i ) / 2 0.2 Large mass −20 10 Box Stack Ball Grid 0.1 Wall Tower −30 0 10 0 100 200 300 0 100 200 300 iterations iterations 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. � K. Erleben c – p. 11/25

Error Correction Penetrations are unavoidable g� Curved surfaces. time t+dt� time t� g� g� Undetected contacts. � K. Erleben c – p. 12/25

Error Correction A first order world simulation: The equations of motion τ = d� L � F = m� and a � (20) dt is replaced by � F = m� v and � τ = I � ω. (21) Used for error correction, yields the scheme, correction () � � (22) given by A = J M − 1 J T , (23a) λ = lcp � A , � d penetration , (23b) s t +1 = � s t + SM − 1 J T � � λ, (23c) � K. Erleben c – p. 13/25

Shock-Propagation The general idea Velocity-based-Shock-Propagation(f,dt)� collision detection at time t� � dynamics(f*dt)� shock-propagration(� Velocity Correction� � dynamics((1-f)*dt)�,� Position Correction� � correction()� )� collision detection at time t + dt� Smoothing� � correction()� t = t + dt� � K. Erleben c – p. 14/25

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. 4� 3� 4� 4� 3� 4� 2� 2� 1� 1� 0� Notice: Stack-height is the path-cost of a breadth-first traversal starting at the fixed bodies. � K. Erleben c – p. 15/25

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 . 4� 3� Edges are given a stack 4� 4� layer number: If stack-heights are Stack Layer 1� 3� equal stack layer 4� 2� number is stack 2� height minus one. Otherwise stack layer 1� 1� number is minimum. 0� � K. Erleben c – p. 16/25

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 � K. Erleben c – p. 17/25

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. � K. Erleben c – p. 18/25

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. � K. Erleben c – p. 19/25

Comparisons Novodex Default � K. Erleben c – p. 20/25

Comparisons Novodex 300 Iterations low separation distance � K. Erleben c – p. 21/25

Comparisons Novodex 30 Iterations, Tweaked Mass, gravity,... � K. Erleben c – p. 22/25

Comparisons Guendelman et. al. (Siggraph 2003) � K. Erleben c – p. 23/25

Comparisons Maya Fatal error... � K. Erleben c – p. 24/25

Comparisons Performance Comparison 0.7 Erleben Guendelman et. al. 0.6 Novodex Default Novodex 30 Novodex 300 0.5 Time (secs) 0.4 0.3 0.2 0.1 0 0 200 400 600 800 1000 Iteration � K. Erleben c – p. 25/25