Paper Reading Note
Implicit Euler Solver
$$ \begin{gather*} \mathbf{q}_{n+1}=\mathbf{q}_n+h\mathbf{v}_{n+1} \\ \mathbf{v}_{n+1}=\mathbf{v}_n+h\mathbf{M}^{-1}(\mathbf{f}_{\text{int}}(\mathbf{q}_{n+1})+\mathbf{f}_{\text{ext}}) \end{gather*} $$$$ \require{mathtools} \DeclarePairedDelimiters\norm{\lVert}{\rVert} \min_{\mathbf{q}_{n+1}} \dfrac{1}{2h^2} \norm*{\mathbf{M}^{\frac{1}{2}} (\mathbf{q}_{n+1} - \mathbf{s}_n)}_F^2 + \sum_i W_i(\mathbf{q}_{n+1}) $$where $\mathbf{s}_n = \mathbf{q}_n + h\mathbf{v}_n + h^2\mathbf{M}^{-1}\mathbf{f}_{\text{ext}}$ and $W_i(\mathbf{q})$ is a scalar potential energy function.
$$ \begin{equation}\label{potential} \dfrac{1}{2h^2} \norm*{\mathbf{M}^{\frac{1}{2}} (\mathbf{q}_{n+1} - \mathbf{s}_n)}_F^2 \end{equation} $$which states that the solution should follow its momentum (plus external forces), and the elastic potential, that requires the solution to minimize the elastic deformation. The corresponding weighting terms, i.e. the mass distribution in $M$, the time step $h$ and the material stiffness of $W$, determine which potential has more importance in this balance.
Nonlinear Elasticity
In nonlinear continuum mechanics elastic potentials are $W(\mathbf{q}) = \Psi (\mathbf{E}(\mathbf{q}))$, where $\mathbf{E}(\mathbf{q})$ is the strain. We can observe that $\mathbf{E}(\mathbf{q}) = 0$ defines a constraint manifold of all possible undeformed configurations, while $\Psi (\mathbf{E}(\mathbf{q}))$ measures how far the deformed configuration is from this manifold (like level set). Therefore, these two concepts can be decoupled; the distance metric does not have to be a complicated nonlinear function because the nonlinearities are already captured by the constraint manifold.
$$ W(\mathbf{q}, \mathbf{p}) = d(\mathbf{q}, \mathbf{p}) + \delta_\mathbf{E}(\mathbf{p}) $$Here, $\delta_\mathbf{E}(\mathbf{p})$ is an indicator function that evaluates to zero if $\mathbf{E}(\mathbf{p}) = 0$ and to $+\infty$ otherwise, and formalizes the requirement that $\mathbf{p}$ should lie on the constraint manifold. The function $d(\mathbf{q}, \mathbf{p})$ then measures a distance between ${\bf q}$ and ${\bf p}$. Minimizing it over $\mathbf{p}$ corresponds to a projection of $\mathbf{q}$ onto the constraint manifold. An analogous elastic potential can be defined as $\tilde{W}(\mathbf{q}) = \min_\mathbf{p} W(\mathbf{q}, \mathbf{p})$
Quadratic distance measures:
$$ W(\mathbf{q}, \mathbf{p}) = \frac{w}{2}\norm*{\mathbf{Aq} - \mathbf{Bp}}_F^2 + \delta_\mathbf{C}(\mathbf{p}) $$where $\mathbf{A}$ and $\mathbf{B}$ are constant matrices and $w$ is a nonnegative weight. We are not restricted to Green’s strains but can use any constraint definition $\mathbf{C}(\mathbf{q}) = 0$ for the set of desired configurations.
$$ \frac{1}{2h^2} \norm*{\mathbf{M}^{\frac{1}{2}} (\mathbf{q} - \mathbf{s}_n)}_F^2 + \sum_i \frac{w_i}{2} \norm*{\mathbf{A}_i \mathbf{S}_i \mathbf{q} - \mathbf{B}_i \mathbf{p}_i }_F^2 + \delta_{\mathbf{C}_i}(\mathbf{p}_i) $$where $\mathbf{S}_i$ is a constant selection matrix that selects the vertices involved in the $i$-th constraint.
$$ \min_{\mathbf{p}_i} \dfrac{w_i}{2} \norm*{\mathbf{A}_i \mathbf{S}_i \mathbf{q} - \mathbf{B}_i \mathbf{p}_i }_F^2 + \delta_{\mathbf{C}_i}(\mathbf{p}_i) $$$$ \left(\frac{\mathbf{M}}{h^2} + \sum_i w_i \mathbf{S}_i^T \mathbf{A}_i^T \mathbf{A}_i \mathbf{S}_i \right) \mathbf{q} = \frac{\mathbf{M}}{h^2} \mathbf{s}_n + \sum_i w_i \mathbf{S}_i^T \mathbf{A}_i^T \mathbf{B}_i \mathbf{p}_i $$The system matrix (left hand side) is constant as long as the constraints are not changing and therefore can be prefactored at initializations. The right hand side requires recomputation in each iteration after the projection variables have been updated in the local step.
Can choose $\mathbf{A}_i = \mathbf{B}_i$ as differential coordinate matrices to improve convergence since the internal physical constraints are translation invariant. But needed to be derived from continuum formulations to avoid tessellation and resolution dependent. Note that the choice of $\mathbf{A}_i$ and $\mathbf{B}_i$ only impacts the numerical solution procedure and does not affect the conservation of momentum.
PBD View
$$ \dfrac{1}{2} \sum_i \norm*{\mathbf{M}^{\frac{1}{2}} (\mathbf{S}_i \mathbf{q} - \mathbf{p}_i)}_F^2 + \delta_{C_i}(\mathbf{p}_i) $$$$ \dfrac{1}{2} \norm*{\mathbf{M}^{\frac{1}{2}} \Delta \mathbf{q}}_F^2 + \lambda \left(C(\mathbf{q}) + \operatorname{tr}(\nabla C(\mathbf{q})^T \Delta \mathbf{q})\right) $$$$ \Delta \mathbf{q} = - \mathbf{M}^{-1} \nabla C(\mathbf{q}) \frac{C(\mathbf{q})}{\norm*{\mathbf{M}^{-\frac{1}{2}} \nabla C(\mathbf{q})}_F^2} $$corresponding exactly to the mass-weighted update rule of PBD.
Jacobi Solver:
$$ \dfrac{1}{2h^2} \norm*{\mathbf{M}^{\frac{1}{2}} (\mathbf{q}_{n+1} - \mathbf{s}_n)}_F^2 + \sum_i \dfrac{w_i}{2} \norm*{\mathbf{M}_i^{\frac{1}{2}} (\mathbf{S}_i \mathbf{q} - \mathbf{p}_i)}_F^2 + \delta_{C_i}(\mathbf{p}_i) $$Its connection to projective implicit Euler is by choosing $\mathbf{A}_i = \mathbf{B}_i = \mathbf{M}_i^{\frac{1}{2}}$
Discussion
The rest of the paper implements several constraints from continuum mechanics’ view. Discrete constraints like positional constraints and collisions are also discussed.
This method has a great generality and numerical stability. The only requirement is that the mesh elements of the input model are well behaved in order to compute the discretization of the gradient and Laplace-Beltrami operators of the original manifold.
Comparison with Newton: Newton’s method exhibits quadratic convergence while local/global solvers have linear convergence.
Comparison with PBD: PBD does not include the momentum constraints making the material stiffness dependent of the number of iterations.
Limitation:
- Implicit damping
- Trading hard constraints for simplicity and efficiency. Although they can be approximated by increasing the weight of the constraints, this can degrade the conditioning of the linear system and can result in locking artifacts.