Projective Dynamics: Fusing Constraint Projections for Fast Simulation

liblaf5/18/23About 5 min

Formula

G = \frac{E}{2 (1 + ν)}

$G$ — shear modulus^[1]
$E$ — Young’s modulus
$ν$ — Poisson’s ratio

Energy

W (\vb q, \vb T) = \frac{w}{2} A \norm {\vb X_{f} \vb X_{g}^{- 1} - \vb T}_{F}^{2}

$W$ — potential energy
$w$ (aka. $G$ ) — shear modulus
(constant) $A$ — triangle area / tetrahedron volume
$\vb X_{f} = \bmqty \vb q_{0, 1, 2} - \vb q_{3}$ — deformed shape
(constant) $\vb X_{g} = \bmqty {\tilde{\vb q}}_{0, 1, 2} - {\tilde{\vb q}}_{3}$ — undeformed shape
(auxiliary) $\vb T$ — rotation matrix

Force

\vb f_{i} = \pdv W (\vb q, \vb T) \vb q_{i} = w A \pqty \vb X_{f} \vb X_{g}^{- 1} - \vb T ⊙ \pqty \pdv \vb X_{f} \vb q_{i} \vb X_{g}^{- 1}

$\vb f$ — force
(variable) $\vb q$ — deformed position

Hessian Matrix

\vb H_{i j} = \pdv f_{i} \vb q_{j} = w A \pqty \pdv \vb X_{f} \vb q_{j} \vb X_{g}^{- 1} ⊙ \pqty \pdv \vb X_{f} \vb q_{i} \vb X_{g}^{- 1}

$\vb H$ — Hessian matrix

Tips

$\vb H$ is constant. The non-zero values of $\vb H_{i j}$ only appear in cases when $\vb q_{i}$ and $\vb q_{j}$ are both located within the same tetrahedron.

Detailed Formula

Force

f_{a i} = \pdv W q_{a i} = w A \pqty \vb X_{f} \vb X_{g}^{- 1} - \vb T ⊙ \pqty \pdv \vb X_{f} q_{a i} \vb X_{g}^{- 1}

$i = 0, 1, 2, 3$ — id of vertex
$a = x, y, z$ — axis

\vb X_{f} = [\begin{matrix} q_{x 0} - q_{x 0} & q_{x 0} - q_{x 0} & q_{x 0} - q_{x 0} \\ q_{y 0} - q_{y 0} & q_{y 0} - q_{y 0} & q_{y 0} - q_{y 0} \\ q_{z 0} - q_{z 0} & q_{z 0} - q_{z 0} & q_{z 0} - q_{z 0} \end{matrix}]

\pqty {\pdv \vb X_{f} q_{a i}}_{l m} = 1 ⟺ i = 0, 1, 2 \qc a = x, y, z \qc l = a \qc m = i \pqty {\pdv \vb X_{f} q_{a i}}_{l m} = - 1 ⟺ i = 3 \qc a = x, y, z \qc l = a \qc m = 0, 1, 2

Therefore, when $i = 0, 1, 2$

\begin{aligned} f_{a i} & = \sum_{l n} w A \pqty {\vb X_{f} \vb X_{g}^{- 1} - \vb T}_{l n} \pqty {\pdv \vb X_{f} q_{0 x} \vb X_{g}^{- 1}}_{l n} \\ = \sum_{l m n} w A \pqty {\vb X_{f} \vb X_{g}^{- 1} - \vb T}_{l n} \pqty {\pdv \vb X_{f} q_{0 x}}_{l m} \pqty {\vb X_{g}^{- 1}}_{m n} \\ = \sum_{n} w A \pqty {\vb X_{f} \vb X_{g}^{- 1} - \vb T}_{a n} \pqty {\vb X_{g}^{- 1}}_{i n} \\ = \sum_{n} w A \pqty {\vb X_{f} \vb X_{g}^{- 1} - \vb T}_{a n} \pqty {\vb X_{g}^{- T}}_{n i} \end{aligned}

\bmqty{\vb{f}_0 & \vb{f}_1 & \vb{f}_2} = f_{a i} = w A \pqty{\vb{X}_f \vb{X}_g^{-1} - \vb{T}} \vb{X}_g^{-T}

Hessian Matrix

When $i, j = 0, 1, 2$

\begin{aligned} H_{a i b j} & = \sum_{l n} w A \pqty {\pdv \vb X_{f} q_{a i} \vb X_{g}^{- 1}}_{l n} \pqty {\pdv \vb X_{f} q_{b j} \vb X_{g}^{- 1}}_{l n} \\ = \sum_{l m n k} w A \pqty {\pdv \vb X_{f} q_{a i}}_{l m} \pqty {\vb X_{g}^{- 1}}_{m n} \pqty {\pdv \vb X_{f} q_{b j}}_{l k} \pqty {\vb X_{g}^{- 1}}_{k n} \\ = \sum_{n} w A \pqty {\pdv \vb X_{f} q_{a i}}_{a i} \pqty {\vb X_{g}^{- 1}}_{i n} \pqty {\pdv \vb X_{f} q_{b j}}_{b j} \pqty {\vb X_{g}^{- 1}}_{j n} \pqty l = a = b \\ = \sum_{n} w A \pqty {\vb X_{g}^{- 1}}_{i n} \pqty {\vb X_{g}^{- T}}_{n j} \end{aligned}

\vb H = H_{i j} = w A \vb X_{g}^{- 1} \vb X_{g}^{- T}

When $i = 0, 1, 2 \qc j = 3$

\begin{aligned} H_{a i b j} & = \sum_{l n} w A \pqty {\pdv \vb X_{f} q_{a i} \vb X_{g}^{- 1}}_{l n} \pqty {\pdv \vb X_{f} q_{b j} \vb X_{g}^{- 1}}_{l n} \\ = \sum_{l m n k} w A \pqty {\pdv \vb X_{f} q_{a i}}_{l m} \pqty {\vb X_{g}^{- 1}}_{m n} \pqty {\pdv \vb X_{f} q_{b j}}_{l k} \pqty {\vb X_{g}^{- 1}}_{k n} \\ = \sum_{n k} w A \pqty {\pdv \vb X_{f} q_{a i}}_{a i} \pqty {\vb X_{g}^{- 1}}_{i n} \pqty {\pdv \vb X_{f} q_{b j}}_{b k} \pqty {\vb X_{g}^{- 1}}_{k n} \pqty l = a = b \\ = \sum_{n k} - w A \pqty {\vb X_{g}^{- 1}}_{i n} \pqty {\vb X_{g}^{- T}}_{n k} \end{aligned}

H_{i 3} = \sum_{k} - w A \pqty {\vb X_{g}^{- 1} \vb X_{g}^{- T}}_{i k}

When $i = j = 3$

\begin{aligned} H_{a i b j} & = \sum_{l n} w A \pqty {\pdv \vb X_{f} q_{a i} \vb X_{g}^{- 1}}_{l n} \pqty {\pdv \vb X_{f} q_{b j} \vb X_{g}^{- 1}}_{l n} \\ = \sum_{l m n k} w A \pqty {\pdv \vb X_{f} q_{a i}}_{l m} \pqty {\vb X_{g}^{- 1}}_{m n} \pqty {\pdv \vb X_{f} q_{b j}}_{l k} \pqty {\vb X_{g}^{- 1}}_{k n} \\ = \sum_{m n k} w A \pqty {\pdv \vb X_{f} q_{a i}}_{a m} \pqty {\vb X_{g}^{- 1}}_{m n} \pqty {\pdv \vb X_{f} q_{b j}}_{b k} \pqty {\vb X_{g}^{- 1}}_{k n} \pqty l = a = b \\ = \sum_{m n k} w A \pqty {\vb X_{g}^{- 1}}_{m n} \pqty {\vb X_{g}^{- T}}_{n k} \end{aligned}

H_{33} = \sum_{m k} - w A \pqty {\vb X_{g}^{- 1} \vb X_{g}^{- T}}_{m k}

Solver

\vb s_{n} = \vb q_{n} + h \vb v_{n} + h^{2} \vb M^{- 1} \vb f

$h$ (aka. $Δ t$ ) — time step
$\vb M$ — mass matrix

\begin{aligned} For k = 0 \dots K \\ Solve \pqty \frac{1}{h^{2}} \vb M + \vb H Δ \vb q = - \frac{1}{h^{2}} \vb M \pqty \vb q^{(k)} - \vb s_{n} + \vb f \pqty \vb q^{(k)} \\ \vb q^{(k + 1)} \leftarrow \vb q^{(k)} + Δ \vb q \\ \vb q_{n + 1} = \vb q^{(k + 1)} \\ \vb v_{n + 1} = (\vb q_{n + 1} - \vb q_{n}) / h \end{aligned}

Conjugate gradient method

used to solve $\vb A \vb x = \vb b$
$\vb A = \frac{1}{h^{2}} \vb M + \vb H$
$\vb b = - \frac{1}{h^{2}} \vb M \pqty \vb q^{(k)} - \vb s_{n} + \vb f \pqty \vb q^{(k)}$

\begin{aligned} \vb r_{0} = \vb b - \vb A \vb x_{0} \\ if \vb r_{0} is sufficiently small, then return \vb x_{0} as the result \\ \vb p_{0} = \vb r_{0} \\ repeat \\ α_{k} = \frac{\vb r_{k}^{T} \vb r_{k}}{\vb p_{k}^{T} \vb A \vb p_{k}} = \frac{\vb r_{k} \vdot \vb r_{k}}{\vb p_{k} \vdot \vb A \vb p_{k}} \\ \vb x_{k + 1} = \vb x_{k} + α_{k} \vb p_{k} \\ \vb r_{k + 1} = \vb r_{k} - α_{k} \vb A \vb p_{k} \\ if \vb r_{k + 1} is sufficiently small, then exit loop \\ β_{k} = \frac{\vb r_{k + 1}^{T} \vb r_{k + 1}}{\vb r_{k}^{T} \vb r_{k}} = \frac{\vb r_{k + 1} \vdot \vb r_{k + 1}}{\vb r_{k} \vdot \vb r_{k}} \\ \vb p_{k + 1} = \vb r_{k + 1} + β_{k} \vb p_{k} \\ k = k + 1 \\ end repeat \\ return \vb x_{k + 1} as the result \end{aligned}

Shear modulus - Wikipedia ↩︎