Neo-Hookean elastic material with XPBD

We inherit from multi-affine models.

Extended Position-Based Dynamics (XPBD) is an optimization method designed to iteratively enforce constraints by applying geometric updates directly on point positions. It is particularly suited for simulations requiring stable and realistic physical deformations, and is very robust to long time steps and low iteration counts.

Mathematical Background

XPBD extends classical Position-Based Dynamics (PBD) by incorporating a compliance term to constraints, introducing elasticity into constraints. It has a direct correspondence to well-known energy potential and allows to solve constraints in a time step and iteration count independent manner.

Classical Position-Based Dynamics (PBD)

Mathematically, classical PBD solves constrained optimization problems through iterative projections of the current solution onto constraint manifolds defined by energy potentials:

\[\mathbf{x}^{t+1} \gets \mathbf{x}^{t} + \Delta \mathbf{x}, \quad\text{s.t.}\quad \mathcal{C}(\mathbf{x}^{t+1}) = 0\]

where:

• \(\mathbf{x}^{t}\) is the position of the particles at iteration \(t\).

• \(\mathcal{C}\) represents the vector of constraint functions.

• \(\Delta \mathbf{x}\) is the position correction, ensuring constraint satisfaction.

In classical PBD, each iteration applies a constraint correction step computed as:

\[\Delta \mathbf{x} = k_j s_j \mathbf{M}^{-1}\nabla \mathcal{C}_j(\mathbf{x}^{t})\]

with the scaling factor \(s_j\) calculated through a single Newton iteration step:

\[s_j = \frac{-\mathcal{C}_j(\mathbf{x}^{t})} {\nabla \mathcal{C}_j(\mathbf{x}^{t}) \mathbf{M}^{-1}\nabla \mathcal{C}_j(\mathbf{x}^{t})^T}\]

where:

• \(j\) indexes the constraints.

• \(k_j \in [0, 1]\) is the constraint stiffness.

• \(\mathbf{M}\) is the diagonal mass matrix.

• \(\nabla \mathcal{C}_j(\mathbf{x}^{t})\) denotes the gradient of the constraint \(\mathcal{C}_j\).

XPBD

XPBD extends the classical Position-Based Dynamics by incorporating a compliance parameter \(\alpha_j \geq 0\) for each constraint, introducing elasticity into constraints. Constraints can thus behave elastically, allowing minor violations to improve numerical stability and physical realism:

\[\mathcal{C_j}(\mathbf{x^{t+1}}) = -\tilde\alpha_j \lambda_j^{t+1}\]

where \(\lambda^{t+1}\) is a Lagrange multiplier at iteration \(t+1\). Note that we use here \(\tilde\alpha_j = \alpha_j / \Delta t^2\) to ensure that the compliance parameter is independent of the time step.

The compliance parameter is also directly correlated to the elastic energy potential derived from a constraint \(\mathcal{C}\) by:

\begin{align} U(\mathbf{x}) = \frac{1}{2}\mathcal{C}(\mathbf{x})^T \alpha^{-1}\mathcal{C}(\mathbf{x}) \end{align}

where \(\alpha\) is the block diagonal compliance matrix representing inverse constraint stiffness.

XPBD iteratively solves the constraint system using incremental updates of the Lagrange multipliers \(\lambda_j\). At each iteration \(i\), the incremental Lagrange multiplier \(\Delta \lambda_j\) for a constraint \(j\) is computed by solving the regularized linear equation derived from a Gauss-Seidel style update step:

\[\Delta \lambda_j = \frac{-\mathcal{C}_j(\mathbf{x}^{t}) - \tilde\alpha_j \lambda^t_{j}} {\nabla \mathcal{C}_j(\mathbf{x}^{t}) \mathbf{M}^{-1}\nabla \mathcal{C}_j(\mathbf{x}^{t})^T + \tilde\alpha_j}\]

• \(\lambda^t_{j}\) represents the accumulated Lagrange multiplier at iteration \(t\) for constraint \(j\).

• When \(\alpha_j = 0\), XPBD reduces exactly to classical PBD.

• For \(\alpha_j > 0\), an additional regularization term is included in both the numerator and denominator, limiting the magnitude of the constraint force.

The corresponding position update is computed as:

\[\Delta \mathbf{x} = \Delta \lambda_j \mathbf{M}^{-1}\nabla \mathcal{C}_j(\mathbf{x}^{t}).\]