Position-Based Dynamics Analysis and Implementation Miles Macklin - PowerPoint PPT Presentation

Position-Based Dynamics Analysis and Implementation Miles Macklin

Analysis

Position-Based Dynamics • Very stable • Highly damped • Example

Continuous Equations of Motion • Newton’s second law • Will consider forces which we can derive from an energy potential E(x) M ¨ x = f ( x ) • Our path: start with implicit Euler and transform it into PBD • Why implicit Euler? Also highly stable, damped.


 
 
 
 Implicit Euler Integration • Implicit Euler: 
 v n +1 = v n + ∆ t M − 1 f ( x n +1 ) x n +1 = x n + ∆ t v n +1 • Equivalent to: 
 ✓ x n +1 − 2 x n + x n − 1 ◆ = f ( x n +1 ) M ∆ t 2 • Forces evaluated at end of the time-step • Implicit, position-level, time-discretization of Newton’s equations

Variational Implicit Euler • Discrete equations of motion M ( x n +1 − 2 x n + x n − 1 ) = ∆ t 2 f ( x n +1 ) • Are the first order optimality conditions for a non-linear minimization argmin 1 x ) T M ( x n +1 − ˜ 2( x n +1 − ˜ x ) − ∆ t 2 E ( x n +1 ) • [Goldenthal et al. 2007] 
 [Liu et al. 2013] x = 2 x n − x n − 1 + M − 1 f ext ˜ = x n + ∆ t v n + M − 1 f ext

Variational Implicit Euler • In the limit of infinite 1 x ) T M ( x n +1 − ˜ 2( x n +1 − ˜ x ) − ∆ t 2 E ( x n +1 ) argmin stiffness we obtain a constrained minimization 
 E → ∞ 1 x ) T M ( x n +1 − ˜ 2( x n +1 − ˜ x ) argmin subject to C ( x n +1 ) = 0

Geometric Interpretation x ˜ x n 1 x ) T M ( x n +1 − ˜ 2( x n +1 − ˜ x ) argmin subject to C ( x n +1 ) = 0 x ∗ x 1 • Variational form gives a “step and project” interpretation for C ( x ) = 0 x 2 implicit Euler • PBD performs approximate projection

Solving Discrete constrained equations of motion • Implicit time discretization produces a non-linear M ( x n +1 � ˜ x ) � ∆ t 2 r C ( x n +1 ) T λ = 0 system of equations C ( x n +1 ) = 0 • How do we solve such a system? Non-Linear System • Newton’s method g ( x i , λ i ) = 0 h ( x i , λ i ) = 0

Approximate Newton Step Full Newton System First approximation: 
  r C T �  �  � g ( x i , λ i ) ∆ x K • M = K + O(dt^2) = � h ( x i , λ i ) ∆ λ r C 0 • Common Quasi-Newton simplification 
 Approximate System  r C T �  �  � ∆ x M 0 Second approximation: 
 = � h ( x i , λ i ) ∆ λ r C 0 • Assume g = 0 • True for first iteration (Schur Complement) PBD System • Typically remains small r C ( x i ) M − 1 r C ( x i ) T ⇤ ⇥ ∆ λ = � C ( x i )

Variational Interpretation of Approximate Projection Implicit Euler x ˜ x n 1 x ) T M ( x − ˜ argmin 2( x − ˜ x ) x ∗ subject to C ( x ) = 0 x 1 C ( x ) = 0 PBD (each iteration) x 2 1 2( x − x i ) T M ( x − x i ) argmin subject to C ( x ) = 0

Problems • To arrive at PBD we had to assume infinitely stiff energy potentials • This means PBD converges to an infinitely stiff solution regardless of stiffness coefficient • Stiffness dependent on iteration count and time-step • No concept of total constraint force • Fully implicit -> severe energy dissipation

Iteration Count Dependent Stiffness 20 ITERATIONS 160 ITERATIONS

PBD Extensions • Projective Dynamics [Bouaziz et al. 2014] • XPBD [Macklin et al. 2016] • Second order PBD

XPBD • Instead of assuming infinite stiffness, Potential allow constraints to be compliant E = 1 2 C T ( x n +1 ) α − 1 C ( x n +1 ) • Leads to a modified / regularized non-linear system • Direct correspondence to engineering Compliance stiffness (Young’s modulus) α = k − 1 • Compliance is simply inverse stiffness • [Servin et al. 2006]

XPBD Newton Step Modified Newton System • Take Schur complement of approximate system  r C T �  �  � ∆ x M 0 = � with respect to M ˜ h ( x i , λ i ) ∆ λ r C α • Obtain PBD or Fast Projection form Schur complement r C ( x i ) M − 1 r C ( x i ) T + ˜ ⇥ ⇤ ∆ λ = � C ( x i ) � ˜ α αλ i • [Goldenthal et al 2007]

XPBD Gauss-Seidel Update PBD • View PBD “scaling fator” s as incremental Lagrange � C j ( x i ) s j = multiplier r C j M − 1 r C T j • Additional compliance terms XPBD • Must store Lagrange multiplier for each constraint � C j ( x i ) � ˜ α j λ ij ∆ λ j = r C j M − 1 r C T j + ˜ α j • PBD solves the infinite stiffness case

XPBD Algorithm • Only two differences from PBD: ‣ Lagrange multiplier calculation (include compliance terms) ‣ Lagrange multiplier update (store instead of discard)

RESULTS • Contact / friction

XPBD - FEM Elastic Energy Potential E tri = V 1 2 ✏ T K ✏ • Generalizes to arbitrary constitutive models Constraint Vector • Treat strain as vector of   ✏ x constraints C tri ( x ) = ✏ tri = ✏ y   ✏ xy • Compliance matrix is inverse Compliance Matrix stiffness − 1   λ + 2 µ λ 0 α tri = K − 1 = λ λ + 2 µ 0   0 0 2 µ

RESULTS • Contact / friction

Results - XPBD vs Implicit Euler • Compare solver output to a non-linear Newton method • Close agreement for primal and dual variables


 
 
 Second Order Implicit Euler • First order backward Euler (BDF1): 
 v n +1 = v n + ∆ t M − 1 f ( x n +1 ) x n +1 = x n + ∆ t v n +1 • Second order backward Euler (BDF2) v n +1 = 4 3 v n − 1 3 v n − 1 + 2 3 ∆ t M − 1 f ( x n +1 ) x n +1 = 4 3 x n − 1 3 x n − 1 + 2 3 ∆ t v n +1

Second Order PBD • Second order prediction: • First order prediction: x = x n + ∆ t v n + ∆ t 2 M − 1 f ext x =4 3 x n − 1 3 x n − 1 + 8 ˜ 9 ∆ t v n ˜ − 2 9 ∆ t v n − 1 + 4 9 ∆ t 2 M − 1 f ext • First order velocity update: v n +1 = 1 • Second order velocity update: x n +1 − x n ⇤ ⇥ ∆ t  3 � v n +1 = 1 2 x n +1 − 2 x n + 1 2 x n − 1 . ∆ t • See [English 08]

Second Order PBD First Order Second Order

Second Order PBD First Order Second Order

Second Order PBD First Order Second Order

Second Order PBD • Significantly less damping • Positions stay closer to constraint manifold • Requires fewer constraint iterations! • Non-smooth events (contact) need special handling

Implementation

Parallel PBD • Gauss-Seidel inherently serial • Parallel options: ‣ Graph coloring methods ‣ Jacobi methods ‣ Hybrid methods

Graph Coloring Methods • Break constraint graph into independent sets • Solve the constraints in a set in parallel • “Batched” Gauss-Seidel • Requires synchronization between each set • Size of sets decreases -> poor utilisation 3 Color Graph

Jacobi Methods • Process each constraint or particle in parallel • Sum up contributions on each particle Particle-centric approach Constraint-centric approach (gather) (scatter) foreach particle (in parallel ) foreach constraint (in parallel ) { { foreach constraint calculate constraint error { foreach particle calculate constraint error { update delta update delta ( atomically ) } } } }


 Jacobi Methods • Problem: system matrix can be indefinite, Jacobi will not converge, e.g.: for redundant constraints (cf. figure) • Regularized Jacobi iteration via averaging [Bridson et al. 02] • Sum all constraint deltas together and divide by constraint count for that particle 
 x i x i + 1 X λ j r C j n i n i • Successive-over relaxation by user parameter omega [0,2]: x i x i + ω X λ j r C j n i n i

Parallel Methods Comparison Method Advantages Disadvantages Good Convergence Graph Coloring Batched Gauss-Seidel Very Robust Synchronization Slow Convergence Jacobi Trivial Parallelism Less Robust

Hybrid Parallel Methods • Best of both worlds • Perform graph-coloring • Upper limit on number of colors • Process everything else with Jacobi • [Fratarcangeli & Pellacini 2015]

Solver Framework

Unified Solver Everything is a set of particles connected by constraints • Simplifies collision detection • Two-way interaction of all object types: 
 ‣ Cloth ‣ Deformables ‣ Fluids ‣ Rigid Bodies • Fits well on the GPU

Examples • Show some neat examples of what we can do with Flex that would not be possible in other frameworks • Smoke / Cloth • Water / Buoyancy

Particles • Velocity stored explicitly struct Particle { • Phase-ID used to control collision float pos[3]; filtering float vel[3]; float invMass; • Global radius int phase; • SOA layout };

Constraints • Constraint types: ‣ Distance (clothing) ‣ Shape (rigids, plastics) ‣ Density (fluids) ‣ Volume (inflatables) ‣ Contact (non-penetration) • Combine constraints ‣ Melting, phase-changes ‣ Stiff cloth, bent metal

Contact and Friction

Collision Detection Between Particles • All dynamics represented as particles • Kinematic objects represented as meshes C contact = | x i − x j | − 2 r ≥ 0 • Two types of collision detection: ‣ Particle-Particle ‣ Particle-Mesh C contact = n · x − r ≥ 0