Towards Real-time Simulation of Hyperelastic Materials A - - PowerPoint PPT Presentation

Towards Real-time Simulation

• f Hyperelastic Materials

A Dissertation Presentation

Tiantian Liu April 24th 2018

Deformable Body Simulation

Limited Human Interactivity in Mixed Reality Environments

Limited Materials in Real-time Simulators

Robotics

[Bai et al. 2016]

Rigid Body v.s. Deformable Body

Goal: Fast simulation of general deformable objects

Goal: Fast simulation of general deformable objects

Simple

Related Work: Classic work

Towards Real-time Simulation of Deformable Objects 9

[Baraff and Witkin 1998] [Goldenthal et al. 2007] [Tournier et al. 2015]

Simple

Related Work: Position Based Dynamics

10 Quasi-Newton Methods for Real-time Simulation of Hyperelastic Materials

[Müller et al. 2007] [Macklin et al. 2016]

Simple

Simulation: Prediction of Future

ℎ 𝑦0 𝑦𝑜 𝑦𝑜+1

Temporal Discretization

Newton’s 2nd Law of Motion

𝑦𝑜+1 = 𝑦𝑜 +

𝑢𝑜 𝑢𝑜+ℎ 𝑤 𝑢 𝑒𝑢

𝑤𝑜+1 = 𝑤𝑜 +

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

+ 𝑔

𝑓𝑦𝑢 𝑒𝑢

Temporal Discretization

Implicit Euler Integration

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

𝑗𝑜𝑢 𝑦𝑜+1 + 𝑔 𝑓𝑦𝑢 𝑦𝑜+1 = 𝑦𝑜 + ℎ𝑤𝑜 + ℎ2𝑁−1𝑔

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

𝑧

Temporal Discretization

Variational Implicit Euler

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

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

inertia elasticity 𝑕(𝑦)

Typical workflow of a deformable body simulation

Spatial Discretization Temporal Discretization

min

𝑦

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

Numerical Solution

𝛼𝒉 𝒚

𝒚

Numerical Solution

𝒉 𝒚

𝛼2𝒉 𝒚 ∆𝒚 = −

−1

Numerical Solution: Newton's Method

min

𝑦

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

    Slow

 𝛼2𝐹 depends on 𝑦

 Non-convex

 The Hessian 𝑁/ℎ2 + 𝛼2𝐹 can be indefinite



Non-convex Potential

t

Numerical Solution: Newton's Method

t

Ideal Numerical Problems

Large Convex Quadratic Problem (Ideally with Constant System Matrix) Many Small Non-convex Problems (Ideally Independent)

Mass-spring Systems

Fast Simulation of Mass-Spring Systems Tiantian Liu, Adam W. Bargteil, James F . O'Brien, Ladislav Kavan ACM Transactions on Graphics 32(6) [Proceedings

• f SIGGRAPH Asia], 2013

Mass-spring System: Basis

Hooke’s Law: 𝐹 𝒚𝟐, 𝒚𝟑 = 1 2 𝑙 𝒚𝟐 − 𝒚𝟑 − 𝑚0 2

𝒚1 𝒚2 𝑚0

Hooke’s Law with auxiliary variables

For the j-th spring: 𝐹j 𝒚 =

1 2 𝑙j

𝒚j1 − 𝒚j2 − lj0

2

Introduce auxiliary variable 𝐪j where 𝐪j

= lj0

𝒚𝑘1 𝒚𝑘2 𝒒𝑘

Hooke’s Law with auxiliary variables

𝒒𝑗1 𝒒𝑗2 𝒆𝑗

 min

𝒒𝑘 =𝑚𝑘0 1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘 2

=

1 2 𝑙𝑗

𝒚𝑘1 − 𝒚𝑘2 − 𝑚𝑘0

2

When 𝒒𝑘 = 𝑚𝑘0

𝒚𝑘1−𝒚𝑘2 𝒚𝑘1−𝒚𝑘2

Hooke’s Law with auxiliary variables

𝒚𝑗1 𝒚𝑗2

𝐹 𝒚 =

𝑘

min

𝒒𝑘 =𝑚0𝑘

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘

2

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

1 2 𝑙𝑘 𝒚𝑘1

𝑈 − 𝒚𝑘2 𝑈 − 𝒒𝑘 𝑈 2

Hooke’s Law with auxiliary variables

𝒚𝑘1

𝑈 − 𝒚𝑘2 𝑈 : Discrete Shape Descriptor

𝒚𝑘1

𝑈 − 𝒚𝑘2 𝑈 = 𝑯𝒌 𝑼𝒚

⋮ 1 ⋮ −1 ⋮ 𝑯𝒌 ∈ ℝ𝒐×𝟐 𝒌𝟐 𝒌𝟑 𝒚1

𝑈

𝒚2

𝑈

⋮ ⋮ ⋮ ⋮ ⋮ 𝒚𝑜

𝑈

𝒚 ∈ ℝ𝒐×𝟒

Hooke’s Law with auxiliary variables

𝒒𝑘

𝑈: Projection

𝒒𝑘

𝑈 = 𝑻𝒌 𝑼𝒒

⋮ 1 ⋮ ⋮ ⋮ 𝑻𝒌 ∈ ℝ𝒏×𝟐 𝒌 𝒒1

𝑈

𝒒2

𝑈

⋮ ⋮ ⋮ ⋮ ⋮ 𝒒𝑜

𝑈

𝒒 ∈ ℝ𝒏×𝟒

Hooke’s Law with auxiliary variables

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

1 2 𝑙𝑘 𝒚𝑘1

𝑈 − 𝒚𝑘2 𝑈 − 𝒒𝑘 𝑈 2

𝐹 𝒚 = min

𝒒∈ℳ

1 2 𝑢𝑠 𝒚𝑈

𝑘

𝑙𝑘𝑯𝑘𝑯𝑘

𝑈

𝒚 − 𝑢𝑠 𝒚𝑈

𝑘

𝑙𝑘𝑯𝑘𝑻𝑘

𝑈

𝒒 + 𝑫

𝑴 𝑲 𝐹 𝒚 = min

𝒒∈ℳ

1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Variational Time Integration with Auxiliary Variable

min

𝑦

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 𝐹(𝑦) 𝒒 ∈ 𝓝 𝒒𝒌 = 𝑚𝑘0

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Optimization

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

 𝑵, 𝑴, 𝑲, 𝒅 does not depend on 𝒚 or 𝒒  If we fix 𝒚 -> easy to solve for 𝒒  If we fix 𝒒 -> easy to solve for 𝒚  Invites alternate solver (local/global)

Local Step

𝒚𝑗1 𝒚𝑗2 𝒒𝑗

For each spring, project to unit length using the

current 𝒚 to find 𝒒𝒋

Trivially Parallelizable

Global Step

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

System matrix (𝐍/h2 + 𝑴) is: Independent of 𝒚 and 𝒒 (Constant) Positive Definite Thus can be pre-factorized (using e.g. Cholesky)

𝐺𝑗𝑦 𝒒: 𝒚∗ = 𝑵 ℎ2 + 𝑴

−1

𝑵 ℎ2 𝒛 + 𝑲𝒒

Alternating Solver

Large Convex Quadratic Problem (with a Constant System Matrix) Many Small Non-convex Problems

Performance

Our Method Newton’s Method

Performance

Our Method Newton’s Method

Results: Mass-spring Systems

Results: Mass-spring Systems

Results: Mass-spring Systems

Remark: Fast Mass-spring Systems

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝑚𝑘0

2 =

min

𝒒𝑘 =𝑚𝑘0

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘

2

min

𝑦

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 𝐹(𝑦)

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Remark: Fast Mass-spring Systems

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝑚𝑘0

2 =

min

𝒒𝑘 =𝑚𝑘0

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘

2

min

𝑦

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 𝐹(𝑦)

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Formulate IE as an Optimization Problem Formulate Hooke’s Law with an Auxiliary Variable Local/global Solve

Remark: Fast Mass-spring Systems

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝑚𝑘0

2 =

min

𝒒𝑘 =𝑚𝑘0

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘

2

min

𝑦

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 𝐹(𝑦)

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Remark: Fast Mass-spring Systems

min

𝒒𝑘 =𝑚𝑘0

1 2 𝑙𝑘 𝒚𝑘1 − 𝒚𝑘2 − 𝒒𝑘

2

Projective Dynamics

Projective Dynamics: Fusing Constraint Projections for Fast Simulation Sofien Bouaziz, Sebastian Martin, Tiantian Liu, Ladislav Kavan, Mark Pauly ACM Transactions on Graphics 33(4) [Proceedings

• f SIGGRAPH], 2014

Key Idea of Fast Mass-Spring Systems

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

𝑥

𝑘 𝑯𝒌 𝑼𝒚 − 𝒒𝑘 2

𝐸𝑗𝑡𝑑𝑠𝑓𝑢𝑓 𝑇ℎ𝑏𝑞𝑓 𝐸𝑓𝑡𝑑𝑠𝑗𝑞𝑢𝑝𝑠 − 𝑄𝑠𝑝𝑘𝑓𝑑𝑢𝑗𝑝𝑜 2

Other Discrete Shape Descriptors

Rest pose 𝒀 Current pose 𝒚 𝑯𝑼𝒚 𝐹 𝒚 = min

𝒒∈ℳ 𝑘

𝑥

𝑘 𝑯𝒌 𝑼𝒚 − 𝒒𝑘 2

Other Manifolds

𝑯𝒌

𝑼𝒚

𝒒𝒌

𝑯𝒌

𝑼𝒚

𝒒𝑘

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

𝑥

𝑘 𝑯𝒌 𝑼𝒚 − 𝒒𝑘 2

Intuitive Projection Manifold: SO(3)

SO(3) … Best Fit Rotation Matrix “As Rigid As Possible” [Chao et al. 2010]

Intuitive Projection Manifold: SL(3)

SL(3) … Group of Matrices with det = 1 Volume Preservation

More Discrete Shape Descriptors:

Laplace-Beltrami operator

Results: Projective Dynamics

Remark: Projective Dynamics

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

𝑥

𝑘 𝑯𝒌 𝑼𝒚 − 𝒒𝑘 2

 Like before, 𝑵, 𝑴, 𝑲, 𝒅 does not depend on 𝒚 and 𝒒  If we fix 𝒚 -> easy to solve for 𝒒: Projection  If we fix 𝒒 -> easy to solve for 𝒚: 𝒚∗ =

𝑵 ℎ2 + 𝑴 −1 𝑵 ℎ2 𝒛 + 𝑲𝒒

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷

Limitation: Projective Dynamics

𝐸𝑗𝑡𝑑𝑠𝑓𝑢𝑓 𝑇ℎ𝑏𝑞𝑓 𝐸𝑓𝑡𝑑𝑠𝑗𝑞𝑢𝑝𝑠 − 𝑄𝑠𝑝𝑘𝑓𝑑𝑢𝑗𝑝𝑜 2

𝐹 𝒚 = min

𝒒∈ℳ 𝑘

𝑥

𝑘 𝑯𝒌 𝑼𝒚 − 𝒒𝑘 2 Special Requirement for the Energy Representation

More Materials?

Soft ARAP Stiff ARAP

Spline-Based Materials [Xu et al. 2015]

Soft ARAP Stiff ARAP

Polynomial Material [Xu et al. 2015]

Quasi-Newton Methods for Real-Time Simulation of Hyperelastic Materials

Quasi-Newton Methods for Real-Time Simulation

• f Hyperelastic Materials

Tiantian Liu, Sofien Bouaziz, Ladislav Kavan ACM Transactions on Graphics 36(3) [Presented at SIGGRAPH],2017.

Reformulation of Projective Dynamics

min

𝑦∈ℝ𝑜×3,𝒒∈𝓝

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒 + 𝐷 min

𝑦∈ℝ𝑜×3

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒(𝒚) + 1 2 𝑢𝑠 𝒒(𝒚)𝑼𝑻𝒒(𝒚)

𝒉(𝒚)

Reformulation of Projective Dynamics

𝛼𝒉 𝒚 = 𝑵 ℎ2 𝒚 − 𝒛 + 𝑴𝒚 − 𝑲𝒒 𝒚 + 𝜖𝒒 𝒚 𝜖𝒚 : (𝑻𝒒 𝒚 − 𝑲𝑼𝒚)

min

𝑦∈ℝ𝑜×3

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒(𝒚) + 1 2 𝑢𝑠 𝒒(𝒚)𝑼𝑻𝒒(𝒚)

𝒉(𝒚)

Projection Differential

𝜀 𝑯𝑼𝒚 − 𝒒 𝒚

2 = 𝑯𝑼𝒚 − 𝒒 𝒚 𝑈

𝑯𝑼𝜀𝒚 −𝜀𝒒 𝒚 𝑈 𝑯𝑼𝒚 − 𝒒 𝒚 𝑯𝑼𝒚 𝒒(𝒚) 𝜀𝒒 𝒚

Reformulation of Projective Dynamics

𝛼𝒉 𝒚 = 𝑵 ℎ2 𝒚 − 𝒛 + 𝑴𝒚 − 𝑲𝒒 𝒚 + 𝜖𝒒 𝒚 𝜖𝒚 : (𝑻𝒒 𝒚 − 𝑲𝑼𝒚)

(𝑵 ℎ2 + 𝐌)−1𝛼𝒉 𝒚 = 𝒚 −

𝑵 ℎ2 + 𝑴

−1 𝑵

ℎ2 𝒛 + 𝑲𝒒 𝒚∗ 𝒚∗ = 𝒚 − (𝑵/ℎ2 + 𝑴)−1𝛼𝒉 𝒚

min

𝑦∈ℝ𝑜×3

1 2ℎ2 𝑢𝑠 𝑦 − 𝑧 𝑈𝑵 𝑦 − 𝑧 + 1 2 𝑢𝑠 𝒚𝑈𝑴𝒚 − 𝑢𝑠 𝒚𝑈𝑲𝒒(𝒚) + 1 2 𝑢𝑠 𝒒(𝒚)𝑼𝑻𝒒(𝒚)

𝒉(𝒚)

𝒚∗ = 𝒚 − (𝑵/ℎ2 + 𝑴)−1𝛼𝒉 𝒚

𝒚∗ = 𝒚 − 𝜷 𝛼2𝑕(𝒚)

−1𝛼𝑕 𝒚

Compare to one Newton step:

Reformulation of Projective Dynamics

 𝛽: Step size, usually decided by linesearch, typical value is 1.  𝛼2𝑕 𝒚 : Hessian Matrix, 𝑵/ℎ2 + 𝛼2𝐹(𝒚)

Projective Dynamics: A Quasi Newton method applied

• n a special type of energy

Quasi-Newton Formulation

𝒚∗ = 𝒚 − 𝛽(𝑵/ℎ2 + 𝑴)−1𝛼𝒉 𝒚 𝛽 = 1

Supporting More General Materials

𝒚∗ = 𝒚 − 𝛽(𝑵/ℎ2 + 𝑴)−1𝛼𝒉 𝒚 This quasi-Newton formulation can be used for any hyperelastic material, but:

• We need to do line-search
• 𝛽 = 1 only works for Projective Dynamics
• We need to define the proper weights 𝑥𝑗
• 𝑵/ℎ2 + 𝑘 𝑥

𝑘𝑯𝒌𝑯𝒌 𝑼

𝑥

𝑘

Strain-Stress Curve for PD

• 𝑵/ℎ2 + 𝑘 𝑥

𝑘𝑯𝒌𝑯𝒌 𝑼

Strain Stress

Supporting More General Materials

Strain Stress

𝑥

𝑘

• 𝑵/ℎ2 + 𝑘 𝑥

𝑘𝑯𝒌𝑯𝒌 𝑼

Supporting More General Materials

Quasi-Newton Algorithm

Compute Gradient

Quasi-Newton Algorithm

Evaluate Descent Direction

Quasi-Newton Algorithm

Line Search

Quasi-Newton Algorithm

We can do more

L-BFGS Acceleration

Projective Dynamics Quasi-Newton Method Exact Solution

L-BFGS Acceleration

Quasi-Newton Method Projective Dynamics

𝛼𝒉 𝒚

𝒚

Core of Quasi-Newton Methods

𝒉 𝒚

∆𝒚 = −

−1

𝑵 𝒊𝟑 + 𝑴 𝑩

L-BFGS with rest-pose Hessian

L-BFGS with rest-pose Hessian

L-BFGS with Scaled Identity

L-BFGS with updating Hessian

Performance of L-BFGS family

Results: Accuracy

Results: Robustness

Results: Collision

Results: Anisotropy

Results: Spline-Based Materials

Remark

 Our method is:

 General: supports a variety types of hyperelastic materials  Fast: >10x faster compared to Newton’s method to achieve similar accuracy level  Simple: avoids Hessian computation, avoids definiteness fix

Simple

Towards Real-time Simulation of Deformable Objects:

Generalization of Spatial Discretization Models

Fast Mass Spring System Projective Dynamics

Towards Real-time Simulation of Deformable Objects:

Generalization of Material Models + Acceleration

Projective Dynamics Quasi-Newton Methods

Towards Real-time Simulation of Deformable Objects:

What’s Next?

Quasi-Newton Methods ?

Core of Our Methods

∆𝒚 = − 𝑵 𝒊𝟑 + 𝑴

−1

𝛼g(𝐲) ×

𝑵 𝒊𝟑 + 𝑴 =

Core of Our Methods

∆𝒚 = − 𝑵 𝒊𝟑 + 𝑴

−1

𝛼g(𝐲) × =

Time Varying Events

 Collisions  Tearing or Cutting

×

Collisions

Collisions

Collision: Soft Constraint 𝑜 𝑦𝑡 𝑦

𝐹𝑑𝑝𝑚 = 𝑙𝑑𝑝𝑚 2 𝑦 − 𝑦𝑡 𝑈𝑜

2

, 𝑗𝑔 𝑦 − 𝑦𝑡 𝑈𝑜 < 0 , 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓

Collision: Soft Constraint

𝐹𝑑𝑝𝑚 = 𝑙𝑑𝑝𝑚 2 𝑦 − 𝑦𝑡 𝑈𝑜

2

, 𝑗𝑔 𝑦 − 𝑦𝑡 𝑈𝑜 < 0 , 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 𝛼𝐹𝑑𝑝𝑚 = 𝑙𝑑𝑝𝑚 𝑦 − 𝑦𝑡 𝑈𝑜 𝑜 , 𝑗𝑔 𝑦 − 𝑦𝑡 𝑈𝑜 < 0 , 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓 𝛼2𝐹𝑑𝑝𝑚 = 𝑙𝑑𝑝𝑚𝑜𝑜𝑈 , 𝑗𝑔 𝑦 − 𝑦𝑡 𝑈𝑜 < 0 , 𝑝𝑢ℎ𝑓𝑠𝑥𝑗𝑡𝑓

Quasi-Newton Algorithm with Collisions

𝛼𝐹𝑑𝑝𝑚 𝐹𝑑𝑝𝑚

Tearing

Tearing

Tearing

Tearing

Tearing

Quasi-Newton Algorithm with Tearing

Original L

Parallelization

 Local Step / Gradient Evaluation Step ?

 Yes

 Global Step / Decent Direction Evaluation ?

 Dependents on the Linear Solver

Choice of Linear Solver: Direct Solver

× =

Pros: Accurate Fast in CPU if Prefactorized Cons: Hard to Parallelize Memory Consuming

Choice of Linear Solver: Iterative Solvers

[Wang 2015] [Fratarcangeli et al. 2016]

Choice of Linear Solver: CG

Simulating Stiff/Rigid Materials

Simulating Stiff/Rigid Materials

[Image courtesy of FistfulOfTalent.com]

Increasing Stiffness Directly?

Using Hard Constraints

[Tournier et al. 2015]

𝑑 𝐺 = 𝐺 − 𝑆 = 0

Using Hard Constraints (Cont’d)

𝑦 = 𝑆𝑌 + 𝑢

Damping

 Current damping model: post-processing models – Ether drag, PBD damping

[Li et al. 2018]

Other Time Integrators

 More vivid motion?

 Other Integrators

 Implicit Midpoint  Newmark-Beta  BDF2  [Bathe 2007] Integrator

 Energy Momentum Methods

 [Dimitar et al. 2018]

What’s Next?

 Bring Machine Learning to Physics?

[Video courtesy of Junior Rojas]

Forward Physics Inverse Physics

Our Method

A Bigger Picture

A Bigger Picture

Phys-based Simulation

Thank You