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 = 1 2 𝑙 𝒚1 − 𝒚2 − 𝑚0 2

𝑌1 𝑌2 𝑦1 𝑦2 𝑚0

12

Representation of a Deformable Body

Ψ 𝐆(𝒚) = 𝜈 𝛝(𝐆) 𝐺

2 + 1

2 𝜇𝑢𝑠2(𝛝(𝐆))

13

Deformation Gradient Strain Tensor Energy Density

Elastic Energy

Rest Pose 𝐹 𝒚 = 0 Deformed Pose 𝐹 𝒚 > 0 𝐹 𝒚 ↓

14

Temporal Discretization

ℎ 𝑦0 𝑦𝑜 𝑦𝑜+1

15

Newton’s 2nd Law of Motion

• 𝑦𝑜+1 = 𝑦𝑜 + ׬

t𝑜 t𝑜+ℎ 𝑤 𝑢 𝑒𝑢

• 𝑤𝑜+1 = 𝑤𝑜 + ׬

t𝑜 t𝑜+ℎ 𝑁−1 𝑔 𝑗𝑜𝑢 𝑦 𝑢

+ 𝑔

𝑓𝑦𝑢 𝑒𝑢

𝑏(𝑢)

16

Time integration: Implicit Euler

• 𝑦𝑜+1 = 𝑦𝑜 + ℎ𝑤𝑜+1
• 𝑤𝑜+1 = 𝑤𝑜 + ℎ𝑁−1 𝑔

𝑗𝑜𝑢 𝑦𝑜+1 + 𝑔 𝑓𝑦𝑢

• 𝑦𝑜+1 = 𝑦𝑜 + ℎ𝑤𝑜 + ℎ2𝑁−1𝑔

𝑓𝑦𝑢 + ℎ2𝑁−1𝑔 𝑗𝑜𝑢(𝑦𝑜+1)

𝑧

17

Variational Implicit Euler

• 𝑦𝑜+1 = 𝑏𝑠𝑕𝑛𝑗𝑜𝑦

1 2ℎ2 𝑦 − 𝑧 𝑁 + 𝐹(𝑦)

• Pick your favorite optimization tool to solve
• Gradient Descent / Newton / Quasi-Newton etc…

inertia elasticity

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 = 𝑦1 − 𝑦2 𝑚0 − 1

𝑌1 𝑌2 𝑦1 𝑦2 𝑚0

Goal: 𝑑 𝒚1, 𝒚2 = 0

23

PBD: Pipeline

Init: 𝑦𝑜+1 = 𝑦𝑜 + ℎ𝑤𝑜 + ℎ2𝑁−1𝑔

𝑓𝑦𝑢

Project: Move 𝑦𝑜+1 to satisfy each constraint

24

PBD: Pipeline (Illustrated version)

𝑦𝑜 𝑦𝑜+1 = 𝑧

Projection

𝑦𝑜+1

25

PBD: Projection

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

26

• Construct 𝜀𝑦𝑘 for the j-th constraint as:
• 𝜀𝑦𝑘 = −𝑁

𝑘 −1∇𝑑 𝑘 𝑈𝜀𝜇𝑘

• Compute the step size 𝜀𝜇𝑘 using 𝑑

𝑘 𝑦 + 𝛼𝑑 𝑘𝜀𝑦𝑘 = 0

• 𝜀𝜇𝑘 =

𝑑𝑘 𝛼𝑑𝑘𝑁𝑘

−1𝛼𝑑𝑘 𝑈

𝜀𝑦𝑘 = −𝑁

𝑘 −1∇𝑑 𝑘 𝑈𝜀𝜇𝑘

PBD: Projection (Cont’d)

27

PBD Projection Example

𝑦1 𝑦2 𝑚0 𝜀𝑦1 𝜀𝑦2

𝑑 = 𝑦1 − 𝑦2 𝑚0 − 1 ∇𝑑 = 1 𝑚0 𝑦1 − 𝑦2 𝑈 𝑦1 − 𝑦2 , − 1 𝑚0 𝑦1 − 𝑦2 𝑈 𝑦1 − 𝑦2 𝐍 = 𝑛1 𝑛2 ⨂𝐉 𝜀𝜇 = 𝑑 ∇𝑑𝐍−1𝛼𝑑 = 𝑚0 𝑛1

−1 + 𝑛2 −1 ( 𝑦1 − 𝑦2 − 𝑚0)

28

PBD Projection Example

𝑦1 𝑦2 𝑚0 𝜀𝑦1 𝜀𝑦2 𝜀𝑦1 = − 𝑚0 𝑛1

−1 + 𝑛2 −1

𝑦1 − 𝑦2 − 𝑚0 𝑛1

−1

𝑚0 𝑦1 − 𝑦2 𝑦1 − 𝑦2 𝜀𝑦2 = − 𝑚0 𝑛1

−1 + 𝑛2 −1

𝑦1 − 𝑦2 − 𝑚0 𝑛2

−1

𝑚0 − 𝑦1 − 𝑦2 𝑦1 − 𝑦2 𝜀𝑦 = −𝐍−1∇𝑑𝑈𝜀𝜇

29

PBD Projection Example

𝑦1 𝑦2 𝑚0 𝜀𝑦1 𝜀𝑦2 𝜀𝑦1 = − 𝑛1

−1

𝑛1

−1 + 𝑛2 −1

𝑦1 − 𝑦2 − 𝑚0 𝑦1 − 𝑦2 𝑦1 − 𝑦2 𝜀𝑦2 = + 𝑛2

−1

𝑛1

−1 + 𝑛2 −1

𝑦1 − 𝑦2 − 𝑚0 𝑦1 − 𝑦2 𝑦1 − 𝑦2

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
• May break the symmetry
• Jacobi

+ Easy to parallelize

• Converges slower
• Less stable

<-- Colored GS [Fratarcangeli et al. 2016] <-- Symmetric GS <-- Chebyshev Acceleration [Wang 2015] <-- Under Relaxation

39

Problems

40

[Liu et al. 2013]

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 𝑦 − 𝑧 𝑁 + 𝐹(𝑦)

• What if 𝐹(𝑦) is (almost) infinitely stiff?

inertia elasticity

42

Constraint-based Variational Implicit Euler

• min

𝑦 1 2ℎ2 𝑦 − 𝑧 𝑁 2 𝑡. 𝑢. 𝑑 𝑦 = 0

𝑌1 𝑌2 𝑦1 𝑦2 𝑚0 𝑑 𝑦1, 𝑦2 = 𝑦1 − 𝑦2 𝑚0 − 1 = 0

43

Constraint-based Variational Implicit Euler

• min

𝑦 1 2ℎ2 𝑦 − 𝑧 𝑁 2 𝑡. 𝑢. 𝑑 𝑦 = 0

• Optimality Condition:
• 𝑁

ℎ2 𝑦 − 𝑧 + 𝛼 𝑦𝑑 𝑦 𝑈𝜇 = 0

• 𝑑 𝑦 = 0

44

Compliant Constraints

• Re-define energy using constraints:
• 𝐹 𝑦 = Σ𝑓𝑗𝜈𝑗𝑑𝑗 𝑦 2 =

1 2 𝑑 𝑦 𝑈𝐋𝑑(𝑦)

𝐋 = 𝜷−1

Stiffness Matrix Compliance Matrix

= 1 2 𝑑 𝑦 𝑈𝜷−𝟐𝑑(𝑦)

45

Variational Implicit Euler with compliant constraints

• min

𝑦 1 2ℎ2 𝑦 − 𝑧 𝑁 2 + 1 2 𝑑 𝑦 𝑈𝜷−𝟐𝑑(𝑦)

• Optimality Condition
• 𝑁

ℎ2 𝑦 − 𝑧 + 𝛼 𝑦𝑑 𝑦 𝑈𝜷−1𝑑 𝑦 = 0

• Optimality Condition with 𝜇 = 𝜷−1𝑑 𝑦
• 𝑁

ℎ2 𝑦 − 𝑧 + 𝛼 𝑦𝑑 𝑦 𝑈𝜇 = 0

• 𝑑 𝑦 − 𝜷𝜇 = 0

46

Numerical Solutions

• 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]

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

47

Newton-Raphson method

• For a given state 𝑦(𝑙) and 𝜇(𝑙)
• Compute Newton-Raphson direction 𝜀𝑦 and 𝜀𝜇 using:

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + 𝛼

𝑦𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙

48

Hard Constraints

• 𝜷 = 𝟏

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0 𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 = 0 min

𝑦

1 2ℎ2 𝑦 − 𝑧 𝑁

2

𝑡. 𝑢. 𝑑 𝑦 = 0

49

Hard Constraints: Geometric Interpretation

𝑑 𝑦 = 0 𝑧 𝑦∗

min

𝑦

1 2ℎ2 𝑦 − 𝑧 𝑁

2

𝑡. 𝑢. 𝑑 𝑦 = 0

50

Step and Project (SAP) [Goldenthal et al. 2007]

𝑑 𝑦 = 0 𝑧 𝑦∗ 𝑑 𝑦 = 0 𝑧 𝑦(2) 𝑦(3)

51

Step and Project (SAP) [Goldenthal et al. 2007]

52

𝑑 𝑦 = 0 𝑧 𝑑 𝑦 = 0 𝑧

𝑦𝑙+1 = min

𝑦

1 2ℎ2 𝑦 − 𝑦(𝑙)

𝑁 2

𝑡. 𝑢. 𝑑 𝑦(𝑙) + 𝛼

𝑦𝑑(𝑦 𝑙 )(𝑦 − 𝑦(𝑙)) = 0

min

𝑦

1 2ℎ2 𝑦 − 𝑧 𝑁

2

𝑡. 𝑢. 𝑑 𝑦 = 0

SAP: Update

53

𝑁 ℎ2 𝑦 − 𝑦(𝑙) + 𝛼

𝑦𝑑 𝑦(𝑙) 𝑈𝜇 = 0

𝑑 𝑦(𝑙) + 𝛼

𝑦𝑑(𝑦 𝑙 )(𝑦 − 𝑦(𝑙)) = 0

Initialize 𝑦(𝑙+1) and 𝜇(𝑙+1) with 𝑦(𝑙) and 0

𝑁/ℎ2 𝛼𝑑𝑈 𝛼𝑑 𝜀𝑦 𝜀𝜇 = − 𝑑(𝑦 𝑙 )

Schur Complement of 0: −𝛼𝑑 (𝑁/ℎ2)−1𝛼𝑑𝑈 [−𝛼𝑑 𝑁/ℎ2)−1𝛼𝑑𝑈 𝜀𝜇 = −𝑑(𝑦 𝑙 )

SAP: Update

54

𝑁 ℎ2 𝑦 − 𝑦(𝑙) + 𝛼

𝑦𝑑 𝑦(𝑙) 𝑈𝜇 = 0

𝑑 𝑦(𝑙) + 𝛼

𝑦𝑑(𝑦 𝑙 )(𝑦 − 𝑦(𝑙)) = 0

Initialize 𝑦(𝑙+1) and 𝜇(𝑙+1) with 𝑦(𝑙) and 0

𝑁/ℎ2 𝛼𝑑𝑈 𝛼𝑑 𝜀𝑦 𝜀𝜇 = − 𝑑(𝑦 𝑙 ) 𝜀𝜇 = 1 ℎ2 𝛼𝑑 𝑁−1𝛼𝑑𝑈 −1𝑑 𝑦 𝑙 𝜀𝑦 = −ℎ2𝑁−1𝛼𝑑𝑈𝜀𝜇

SAP: Solution Drifting

55

𝑑 𝑦 = 0 𝑧 𝑦∗ 𝑑 𝑦 = 0 𝑧 𝑦(2) 𝑦(3)

SAP vs. the Exact Solution

56

SAP Exact Solution

PBD vs. SAP

57

• For each constraint j:
• Compute 𝜀𝜇𝑘:𝜀𝜇𝑘 =

𝑑𝑘 𝛼𝑑𝑘 𝑁𝑘

−1𝛼𝑑𝑘 𝑈

• Compute 𝜀𝑦𝑘: 𝜀𝑦𝑘 = −𝑁

𝑘 −1𝛼𝑑 𝑘 𝑈𝜀𝜇j

• Commit: 𝑦 = 𝑦 + 𝜀𝑦𝑘
• 𝜀𝜇 =

1 ℎ2 𝛼𝑑 𝑁−1𝛼𝑑𝑈 −1𝑑

• 𝜀𝑦 = −ℎ2𝑁−1𝛼𝑑𝑈𝜀𝜇
• 𝑦 = 𝑦 + 𝜀𝑦

PBD (G-S) SAP

PBD v.s. SAP: Geometric Interpretation

𝑑 𝑦 = 0 𝑧 𝑦∗

Ground Truth

𝑑 𝑦 = 0 𝑧 𝑦(2) 𝑦(3)

SAP

𝑑 𝑦 = 0 𝑧

PBD

58

PBD Problems Visualized

• Solution Drifts
• Non-physically-based
• Hard to Control

𝑑 𝑦 = 0 𝑧

PBD

59

PBD Problems Visualized

• Solution Drifts
• Non-physically-based
• Stiffness depends on iteration count
• 𝑧 = 𝑦𝑜 + ℎ𝑤𝑜 + ℎ2𝑁−1𝑔

𝑓𝑦𝑢

• Hard to Control

𝑑 𝑦 = 0 𝑧

PBD

60

PBD Problems Visualized

• Solution Drifts
• Non-physically-based
• Hard to Control
• min

𝑦 1 2ℎ2 𝑦 − 𝑧 𝑁 2

• 𝑡. 𝑢. 𝑑 𝑦 = 0

𝑑 𝑦 = 0 𝑧

PBD

61

Compliant Constraints

• 𝜷 ! = 𝟏, 𝐹 =

1 2 𝒅𝑈𝜷−1𝒅

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

62

Compliant Constraints

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

• 𝜷 ! = 𝟏, 𝐹 =

1 2 𝒅𝑈𝜷−1𝒅

• Mass Spring System

𝐹 𝒚1, 𝒚2 = 1 2 𝑙 𝒚1 − 𝒚2 − 𝑚0 2

𝑌1 𝑌2 𝑦1 𝑦2 𝑚0

𝑑 𝒚1, 𝒚2 =

𝒚1−𝒚2 𝑚0

− 1 𝑙 = 2𝜈𝐵 𝑚0 𝛽 = 2𝜈𝐵𝑚0 −1 = 𝑙𝑚0

2 −1

63

Compliant Constraints

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

• 𝜷 ! = 𝟏, 𝐹 =

1 2 𝒅𝑈𝜷−1𝒅

• FEM (3D case)

𝒅 = 𝜗𝑦 𝜗𝑧 𝜗𝑨 𝜗𝑦𝑧 𝜗𝑦𝑨 𝜗𝑧𝑨 𝑈

𝐋 = 𝑊 𝜇 + 2𝜈 𝜇 𝜇 𝜇 𝜇 + 2𝜈 𝜇 𝜇 𝜇 𝜇 + 2𝜈 𝜈 𝜈 𝜈

𝜷 = 𝐋−1

𝐹 𝐆(𝒚) = 𝑊(𝜈 𝛝 𝐺

2 + 1

2 𝜇𝑢𝑠2 𝛝 ) 𝛝 = 𝜗𝑦 0.5𝜗𝑦𝑧 0.5𝜗𝑦𝑨 0.5𝜗𝑦𝑧 𝜗𝑧 0.5𝜗𝑧𝑨 0.5𝜗𝑦𝑨 0.5𝜗𝑧𝑨 𝜗𝑨

64

𝐹 = 1 2 𝒅T𝐋𝒅

Solve Compliant Constraints with Newton

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + ∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙

SAP Assumption: 𝑧 → 𝑦 𝑙 ∇2𝑑 = 0 SAP Initialization: 𝑦 ← 𝑦 𝑙

𝑧 𝑦(2) 𝑦(3)

65

Solve Compliant Constraints with SAP

𝑁 ℎ2 𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙 𝜀𝜇 = 1 ℎ2 𝛼𝑑 𝑁−1𝛼𝑑𝑈 + ℎ2𝜷

−1 𝑑 𝑦 𝑙

− 𝜷𝜇(𝑙) 𝜀𝑦 = −ℎ2𝑁−1𝛼𝑑𝑈𝜇(𝑙+1) 𝜇(𝑙+1) = 𝜇(𝑙) + 𝜀𝜇 𝑦(𝑙+1) = 𝑦(𝑙+1) + 𝜀𝑦

66

Solve Compliant Constraints with XPBD [Macklin et al. 2016]

𝑁 ℎ2 𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙 𝐺𝑝𝑠 𝑓𝑏𝑑ℎ 𝑑𝑝𝑜𝑡𝑢𝑠𝑏𝑗𝑜𝑢: 𝑘

Compute: 𝜀𝜇𝑘 =

1 ℎ2 𝛼𝑑 𝑘 𝑁 𝑘 −1𝛼𝑑 𝑘 𝑈 + ℎ2𝜷𝑘 −1 𝑑 𝑘 𝑦 𝑙

− 𝜷𝑘𝜇𝑘

(𝑙)

Commit: 𝜇𝑘

(𝑙+1) = 𝜇𝑘 (𝑙) + 𝜀𝜇𝑘

Compute: 𝜀𝑦𝑘 = −ℎ2𝑁

𝑘 −1𝛼𝑑 𝑘 𝑈𝜇𝑘 (𝑙+1)

Commit: 𝑦𝑘

(𝑙+1) = 𝑦𝑘 (𝑙+1) + 𝜀𝑦𝑘

67

PBD v.s. XPBD

• For each constraint j:
• Compute 𝜀𝜇𝑘:𝜀𝜇𝑘 =

1 ℎ2 𝑑𝑘 𝛼𝑑𝑘 𝑁𝑘

−1𝛼𝑑𝑘 𝑈

• Commit: 𝜇𝑘 = 𝜀𝜇𝑘
• Compute 𝜀𝑦𝑘: 𝜀𝑦𝑘 = −ℎ2𝑁

𝑘 −1𝛼𝑑 𝑘 𝑈𝜇𝑘

• Commit: 𝑦 = 𝑦 + 𝜀𝑦𝑘

PBD (G-S) XPBD (G-S)

• For each constraint j:
• Compute 𝜀𝜇𝑘:𝜀𝜇𝑘 =

1 ℎ2 𝑑𝑘−𝛽𝑘𝜇𝑘 𝛼𝑑𝑘 𝑁𝑘

−1𝛼𝑑𝑘 𝑈+ℎ2𝛽𝑘

• Commit: 𝜇𝑘 = 𝜇𝑘 + 𝜀𝜇𝑘
• Compute 𝜀𝑦𝑘: 𝜀𝑦𝑘 = −ℎ2𝑁

𝑘 −1𝛼𝑑 𝑘 𝑈𝜇𝑘

• Commit: 𝑦 = 𝑦 + 𝜀𝑦𝑘

68

XPBD Result

[Macklin et al. 2016]

69

What is left out?

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + ∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

70

Schur Complement of the Upper-left Block

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + ∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙 𝑁 ℎ2 + 𝛼𝑑𝑈𝜷−1𝛼𝑑 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

71

Geometric Stiffness

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + ∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙 𝑁 ℎ2 + 𝛼𝑑𝑈𝜷−1𝛼𝑑 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

72

Geometric Stiffness: Example

73

𝑌1 𝑌2 𝑦1 𝑦2 𝑚0 𝑑 = 𝑦1 − 𝑦2 𝑚0 − 1 𝐹 = 1 2 𝛽−1𝑑2

𝛼

𝑦1,𝑦1 2

𝐹 = 𝛼𝑑𝑈𝛽−1𝛼𝑑 + (𝛽−1𝑑)𝛼

𝑦1,𝑦1 2

𝑑

Geometric Stiffness 𝜇

Stable Constrained Dynamics [Tournier et al. 2015]

• When to use Geometric Stiffness?
• 1. Material highly nonlinear
• 2. Strong stiffness

74

𝑁 ℎ2 + ෍

𝑘=1 𝑛

𝜇𝑘

(𝑙)𝛼2𝑑 𝑘

𝛼𝑑𝑈 𝛼𝑑 −𝜷 𝜀𝑦 𝜀𝜇 = −

𝑁 ℎ2 𝑦(𝑙) − 𝑧 + ∇𝑑𝑈𝜇(𝑙)

𝑑 𝑦 𝑙 − 𝜷𝜇 𝑙

Conclusion

75

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

constraint #iterations/frame per-iteration cost [Goldenthal et al. 2007] hard medium medium [Müller et al. 2007] hard high low [Tournier et al. 2015] compliant

• ne

high [Macklin et al. 2016] compliant high low

Conclusion: PBD

76

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

constraint #iterations/frame per-iteration cost [Goldenthal et al. 2007] hard medium medium [Müller et al. 2007] hard high low [Tournier et al. 2015] compliant

• ne

high [Macklin et al. 2016] compliant high low

Conclusion: PBD

• The Performance of (X)PBD is…
• Similar to one iteration of Jacobi/Gauss-Seidel of SAP

77

Conclusion: SAP

78

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

constraint #iterations/frame per-iteration cost [Goldenthal et al. 2007] hard medium medium [Müller et al. 2007] hard high low [Tournier et al. 2015] compliant

• ne

high [Macklin et al. 2016] compliant high low

Conclusion: Geometric Stiffness

79

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

constraint #iterations/frame per-iteration cost [Goldenthal et al. 2007] hard medium medium [Müller et al. 2007] hard high low [Tournier et al. 2015] compliant

• ne

high [Macklin et al. 2016] compliant high low

Key Takeaway from this Talk

80

• (X)PBD is approximately converging to…

𝑁 ℎ2 𝑦 − 𝑧 + 𝛼

𝑦𝑑 𝑦 𝑈𝜇 = 0

𝑑 𝑦 − 𝜷𝜇 = 0

𝑑 𝑦 = 0 𝑧 𝑦∗

Ground Truth

𝑑 𝑦 = 0 𝑧

PBD

Position Based Dynamics

A fast yet physically plausible method for deformable body simulation

Tiantian Liu GAMES Webinar 03/28/2019