Position Based Dynamics A fast yet physically plausible method for - PowerPoint PPT Presentation

Position Based Dynamics A fast yet physically plausible method for deformable body simulation Tiantian Liu GAMES Webinar 03/28/2019

Position Based Dynamics A fast yet physically plausible method for deformable body simulation Tiantian Liu GAMES Webinar 03/28/2019

Intended Takeaway from this Talk… • For Rookies… • Basic idea of a deformable body simulation pipeline • What is Position Based Dynamics (PBD) • How to implement the basic building blocks of PBD • For Veterans… • A physically correct understanding of PBD • Insights and potential improvements of PBD 3

Major References of this Talk [Müller et al. 2007] [Macklin et al. 2016] [Goldenthal et al. 2007] [Tournier et al. 2015] 4

Computer Graphics (Past and Now) 1986 2017 5

Deformable Body Simulation Game Movie/Animation AR/VR 6

Deformable Body Simulation Game Movie/Animation AR/VR 7

Deformable Body Simulation Game Movie/Animation AR/VR 8

Rigid Body v.s. Deformable Body 9

Deformable Body Simulation [Liu et al. 2017] 10

Spatial Discretization 11

Representation of a Deformable Body 𝑌 1 𝑌 2 𝑦 2 𝑦 1 𝑚 0 𝐹 𝒚 1 , 𝒚 2 = 1 𝒚 1 − 𝒚 2 − 𝑚 0 2 2 𝑙 12

Representation of a Deformable Body 2 + 1 2 𝜇𝑢𝑠 2 (𝛝(𝐆)) Ψ 𝐆(𝒚) = 𝜈 𝛝(𝐆) 𝐺 Energy Density Deformation Strain Tensor Gradient 13

Elastic Energy 𝐹 𝒚 ↓ Deformed Pose Rest Pose 𝐹 𝒚 > 0 𝐹 𝒚 = 0 14

Temporal Discretization 𝑦 𝑜 𝑦 𝑜+1 𝑦 0 ℎ 15

Newton’s 2 nd Law of Motion t 𝑜 +ℎ 𝑤 𝑢 𝑒𝑢 • 𝑦 𝑜+1 = 𝑦 𝑜 + ׬ t 𝑜 t 𝑜 +ℎ 𝑁 −1 𝑔 • 𝑤 𝑜+1 = 𝑤 𝑜 + ׬ 𝑗𝑜𝑢 𝑦 𝑢 + 𝑔 𝑓𝑦𝑢 𝑒𝑢 t 𝑜 𝑏(𝑢) 16

Time integration: Implicit Euler • 𝑦 𝑜+1 = 𝑦 𝑜 + ℎ𝑤 𝑜+1 • 𝑤 𝑜+1 = 𝑤 𝑜 + ℎ𝑁 −1 𝑔 𝑗𝑜𝑢 𝑦 𝑜+1 + 𝑔 𝑓𝑦𝑢 • 𝑦 𝑜+1 = 𝑦 𝑜 + ℎ𝑤 𝑜 + ℎ 2 𝑁 −1 𝑔 𝑓𝑦𝑢 + ℎ 2 𝑁 −1 𝑔 𝑗𝑜𝑢 (𝑦 𝑜+1 ) 𝑧 17

Variational Implicit Euler 1 • 𝑦 𝑜+1 = 𝑏𝑠𝑕𝑛𝑗𝑜 𝑦 2ℎ 2 𝑦 − 𝑧 𝑁 + 𝐹(𝑦) inertia elasticity • Pick your favorite optimization tool to solve • Gradient Descent / Newton / Quasi-Newton etc … 18

Deformable Body Simulation [Liu et al. 2017] 19

State-of-the-Art Real-time Simulators NVIDIA FleX [Video courtesy of NVIDIA] 20

State-of-the-Art Real-time Simulators Maya nDynamics [Video courtesy of Pixclue Studios] 21

Position Based Dynamics 22

Position Based Dynamics 𝑌 1 𝑌 2 𝑦 2 𝑦 1 𝑚 0 𝑦 1 − 𝑦 2 𝑑 𝒚 1 , 𝒚 2 = − 1 𝑚 0 Goal: 𝑑 𝒚 1 , 𝒚 2 = 0 23

PBD: Pipeline Init: 𝑦 𝑜+1 = 𝑦 𝑜 + ℎ𝑤 𝑜 + ℎ 2 𝑁 −1 𝑔 𝑓𝑦𝑢 Project: Move 𝑦 𝑜+1 to satisfy each constraint 24

PBD: Pipeline (Illustrated version) Projection 𝑦 𝑜 𝑦 𝑜+1 = 𝑧 𝑦 𝑜+1 25

PBD: Projection • Find a projection direction 𝜀𝑦 to: • Satisfy 𝑑 𝑦 + 𝜀𝑦 ≈ 𝑑(𝑦) + 𝛼𝑑 𝑦 𝜀𝑦 = 0 • Conserve linear momentum: Σ𝑛 𝑗 𝜀𝑦 𝑗 = 0 • Conserve angular momentum: Σ𝑛 𝑗 𝑦 𝑗 × 𝜀𝑦 𝑗 = 0 26

PBD: Projection (Cont’d) • Construct 𝜀𝑦 𝑘 for the j-th constraint as: −1 ∇𝑑 −1 ∇𝑑 𝑈 𝜀𝜇 𝑘 𝑈 𝜀𝜇 𝑘 • 𝜀𝑦 𝑘 = −𝑁 𝜀𝑦 𝑘 = −𝑁 𝑘 𝑘 𝑘 𝑘 • Compute the step size 𝜀𝜇 𝑘 using 𝑑 𝑘 𝑦 + 𝛼𝑑 𝑘 𝜀𝑦 𝑘 = 0 𝑑 𝑘 • 𝜀𝜇 𝑘 = −1 𝛼𝑑 𝑘 𝑈 𝛼𝑑 𝑘 𝑁 𝑘 27

PBD Projection Example 𝜀𝑦 1 𝜀𝑦 2 𝑦 1 𝑦 2 𝑚 0 𝑦 1 − 𝑦 2 𝑈 𝑦 1 − 𝑦 2 𝑈 ∇𝑑 = 1 , − 1 𝑦 1 − 𝑦 2 𝐍 = 𝑛 1 0 𝑑 = − 1 𝑛 2 ⨂𝐉 0 𝑚 0 𝑦 1 − 𝑦 2 𝑚 0 𝑦 1 − 𝑦 2 𝑚 0 𝑑 𝑚 0 𝜀𝜇 = ∇𝑑𝐍 −1 𝛼𝑑 = −1 ( 𝑦 1 − 𝑦 2 − 𝑚 0 ) −1 + 𝑛 2 𝑛 1 28

PBD Projection Example 𝜀𝑦 1 𝜀𝑦 2 𝑦 1 𝑦 2 𝑚 0 𝜀𝑦 = −𝐍 −1 ∇𝑑 𝑈 𝜀𝜇 −1 𝑚 0 𝑛 1 𝑦 1 − 𝑦 2 𝜀𝑦 1 = − 𝑦 1 − 𝑦 2 − 𝑚 0 −1 + 𝑛 2 −1 𝑚 0 𝑦 1 − 𝑦 2 𝑛 1 −1 𝑚 0 𝑛 2 − 𝑦 1 − 𝑦 2 𝜀𝑦 2 = − 𝑦 1 − 𝑦 2 − 𝑚 0 −1 + 𝑛 2 −1 𝑚 0 𝑦 1 − 𝑦 2 𝑛 1 29

PBD Projection Example 𝜀𝑦 1 𝜀𝑦 2 𝑦 1 𝑦 2 𝑚 0 −1 𝑛 1 𝑦 1 − 𝑦 2 𝜀𝑦 1 = − 𝑦 1 − 𝑦 2 − 𝑚 0 −1 + 𝑛 2 −1 𝑦 1 − 𝑦 2 𝑛 1 −1 𝑛 2 𝑦 1 − 𝑦 2 𝜀𝑦 2 = + 𝑦 1 − 𝑦 2 − 𝑚 0 −1 + 𝑛 2 −1 𝑦 1 − 𝑦 2 𝑛 1 30

PBD: Pipeline Project: Move 𝑦 𝑜+1 to satisfy each constraint 31

Iteration Strategy: Gauss-Seidel v.s. Jacobi 32

Gauss-Seidel 33

Gauss-Seidel 34

Gauss-Seidel 35

Gauss-Seidel 36

Gauss-Seidel 37

Jacobi 38

Gauss-Seidel v.s. Jacobi • Gauss-Seidel + Converges faster - Hard to parallelize <-- Colored GS [Fratarcangeli et al. 2016] - May break the symmetry <-- Symmetric GS • Jacobi + Easy to parallelize - Converges slower <-- Chebyshev Acceleration [Wang 2015] - Less stable <-- Under Relaxation 39

Problems [Liu et al. 2013] 40

Conclusion • Position Based Dynamics is: • Fast • Simple to implement • Stable (for most cases) • It was also considered: • Non-physically-based (stiffness related to iteration count) • Hard to control 41

Variational Implicit Euler 1 • 𝑦 𝑜+1 = 𝑏𝑠𝑕𝑛𝑗𝑜 𝑦 2ℎ 2 𝑦 − 𝑧 𝑁 + 𝐹(𝑦) inertia elasticity • What if 𝐹(𝑦) is (almost) infinitely stiff? 42

Constraint-based Variational Implicit Euler 2 𝑡. 𝑢. 𝑑 𝑦 = 0 1 • min 2ℎ 2 𝑦 − 𝑧 𝑁 𝑦 𝑌 1 𝑌 2 𝑦 2 𝑦 1 𝑚 0 𝑦 1 − 𝑦 2 𝑑 𝑦 1 , 𝑦 2 = − 1 = 0 𝑚 0 43

Constraint-based Variational Implicit Euler 2 𝑡. 𝑢. 𝑑 𝑦 = 0 1 • min 2ℎ 2 𝑦 − 𝑧 𝑁 𝑦 • Optimality Condition: 𝑁 • 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑧 + 𝛼 • 𝑑 𝑦 = 0 44

Compliant Constraints • Re-define energy using constraints: = 1 • 𝐹 𝑦 = Σ 𝑓 𝑗 𝜈 𝑗 𝑑 𝑗 𝑦 2 = 1 2 𝑑 𝑦 𝑈 𝐋𝑑(𝑦) 2 𝑑 𝑦 𝑈 𝜷 −𝟐 𝑑(𝑦) 𝐋 = 𝜷 −1 Stiffness Matrix Compliance Matrix 45

Variational Implicit Euler with compliant constraints 1 2 + 1 2 𝑑 𝑦 𝑈 𝜷 −𝟐 𝑑(𝑦) • min 2ℎ 2 𝑦 − 𝑧 𝑁 𝑦 • Optimality Condition 𝑁 𝑦 𝑑 𝑦 𝑈 𝜷 −1 𝑑 𝑦 = 0 • ℎ 2 𝑦 − 𝑧 + 𝛼 • Optimality Condition with 𝜇 = 𝜷 −1 𝑑 𝑦 𝑁 • 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑧 + 𝛼 • 𝑑 𝑦 − 𝜷𝜇 = 0 46

Numerical Solutions 𝑁 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑧 + 𝛼 𝑑 𝑦 − 𝜷𝜇 = 0 • To achieve the exact solution • Newton-Raphson • To achieve an approximated solution • Step-and-Project [Goldenthal et al. 2007] • Semi-Implicit Euler [Tournier et al. 2015] • Position Based Dynamics [Müller et al. 2007, Macklin et al. 2016] 47

Newton-Raphson method 𝑁 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑧 + 𝛼 𝑑 𝑦 − 𝜷𝜇 = 0 • For a given state 𝑦 (𝑙) and 𝜇 (𝑙) • Compute Newton-Raphson direction 𝜀𝑦 and 𝜀𝜇 using: 𝑛 𝑁 𝑁 ℎ 2 𝑦 (𝑙) − 𝑧 + 𝛼 (𝑙) 𝛼 2 𝑑 𝑦 𝑑 𝑈 𝜇 (𝑙) 𝛼𝑑 𝑈 ℎ 2 + ෍ 𝜇 𝑘 𝜀𝑦 𝑘 𝜀𝜇 = − 𝑘=1 𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙 𝛼𝑑 −𝜷 48

Hard Constraints 𝑁 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑧 + 𝛼 𝑑 𝑦 − 𝜷𝜇 = 0 • 𝜷 = 𝟏 1 𝑁 2 𝑦 𝑑 𝑦 𝑈 𝜇 = 0 min 2ℎ 2 𝑦 − 𝑧 𝑁 ℎ 2 𝑦 − 𝑧 + 𝛼 𝑦 𝑡. 𝑢. 𝑑 𝑦 = 0 𝑑 𝑦 = 0 49

Hard Constraints: Geometric Interpretation 𝑧 𝑦 ∗ 1 2 min 2ℎ 2 𝑦 − 𝑧 𝑁 𝑦 𝑡. 𝑢. 𝑑 𝑦 = 0 𝑑 𝑦 = 0 50

Step and Project (SAP) [Goldenthal et al. 2007] 𝑧 𝑧 𝑦 (2) 𝑦 (3) 𝑦 ∗ 𝑑 𝑦 = 0 𝑑 𝑦 = 0 51

Step and Project (SAP) [Goldenthal et al. 2007] 𝑧 𝑧 𝑑 𝑦 = 0 𝑑 𝑦 = 0 1 1 2 2 2ℎ 2 𝑦 − 𝑦 (𝑙) min 2ℎ 2 𝑦 − 𝑧 𝑁 𝑦 𝑙+1 = min 𝑁 𝑦 𝑦 𝑡. 𝑢. 𝑑 𝑦 (𝑙) + 𝛼 𝑦 𝑑(𝑦 𝑙 )(𝑦 − 𝑦 (𝑙) ) = 0 𝑡. 𝑢. 𝑑 𝑦 = 0 52

SAP: Update 𝑁 𝑦 𝑑 𝑦 (𝑙) 𝑈 𝜇 = 0 ℎ 2 𝑦 − 𝑦 (𝑙) + 𝛼 Initialize 𝑦 (𝑙+1) and 𝜇 (𝑙+1) with 𝑦 (𝑙) and 0 𝑑 𝑦 (𝑙) + 𝛼 𝑦 𝑑(𝑦 𝑙 )(𝑦 − 𝑦 (𝑙) ) = 0 0 𝑁/ℎ 2 𝛼𝑑 𝑈 𝜀𝑦 𝜀𝜇 = − 𝑑(𝑦 𝑙 ) 𝛼𝑑 0 Schur Complement of 0 : −𝛼𝑑 (𝑁/ℎ 2 ) −1 𝛼𝑑 𝑈 [−𝛼𝑑 𝑁/ℎ 2 ) −1 𝛼𝑑 𝑈 𝜀𝜇 = −𝑑(𝑦 𝑙 ) 53