Physics-Based Simulation

Continuum Mechanics-Based Constraints

While simple geometric constraints such as distance and volume are effective for many materials, a more physically rigorous and expressive simulation can be achieved by deriving constraints directly from the principles of continuum mechanics. This approach connects Position-Based Dynamics to the rich theoretical foundation of established methods like the Finite Element Method (FEM), allowing for the simulation of complex material behaviors including anisotropy, elastoplasticity, and nonlinear elasticity. The core idea is to define constraints based on the deformation gradient, a tensor that measures the local stretching and shearing of the material.

This section will first review the necessary concepts of discretized deformation and strain. It will then detail two primary methods for formulating continuum-based constraints: a holistic approach based on strain energy and a more granular method that constrains individual components of the strain tensor directly.

The Strain Energy Constraint

The most direct way to incorporate a standard hyperelastic material model into PBD is to use its strain energy potential as a constraint function.

In a discrete setting, such as a tetrahedral mesh, the deformation gradient is assumed to be constant within each element. For a tetrahedron with material-space vertices $({\mathbf{P}}_1,{\mathbf{P}}_2,{\mathbf{P}}_3,{\mathbf{P}}_4)$ and world-space vertices $({\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3,{\mathbf{p}}_4)$, we can define a material-space shape matrix ${\mathbf{D}}_m$ and a world-space shape matrix ${\mathbf{D}}_s$ from the edge vectors: $$\begin{array}{rl}{\mathbf{D}}_m& =[{\mathbf{P}}_1-{\mathbf{P}}_4,{\mathbf{P}}_2-{\mathbf{P}}_4,{\mathbf{P}}_3-{\mathbf{P}}_4]\\ {\mathbf{D}}_s& =[{\mathbf{p}}_1-{\mathbf{p}}_4,{\mathbf{p}}_2-{\mathbf{p}}_4,{\mathbf{p}}_3-{\mathbf{p}}_4]\end{array}$$ The deformation gradient for the tetrahedron is then given by $\mathbf{F}={\mathbf{D}}_s{\mathbf{D}}_m^{-1}$. Note that ${\mathbf{D}}_m^{-1}$ is constant and can be precomputed.

From $\mathbf{F}$, we can derive various strain measures. One of the most common is the Green-St. Venant strain tensor, $\mathbf{\epsilon }$, which is independent of rigid-body rotations: $$\mathbf{\epsilon }=\frac{1}{2}({\mathbf{F}}^T\mathbf{F}-\mathbf{I})\text{\hspace{1em}(33.5.1)}$$ The related right Cauchy-Green deformation tensor, $\mathbf{S}={\mathbf{F}}^T\mathbf{F}=2\mathbf{\epsilon }+\mathbf{I}$, will be particularly useful for formulating direct strain constraints. If there is no deformation, $\mathbf{F}=\mathbf{I}$, $\mathbf{S}=\mathbf{I}$, and $\mathbf{\epsilon }=0$.

By Hooke’s generalized law, we can relate stress to strain:

$$\mathbf{S}=\mathbf{C}\mathbf{\epsilon }$$ where $\mathbf{C}$ is the elasticity tensor of the material. Recall that the energy of a deformed solid is defined by integrating the scalar strain energy density field over the whole body $\Omega$:

$$E_s={\int }_{\Omega }{\Psi }_{s}d\mathbf{P}$$ where the strain energy density can be defined in terms of stress and strain:
$${\Psi }_s=\frac{1}{2}\epsilon :\mathbf{S}=\frac{1}{2}tr({\epsilon }^{\top }\mathbf{S})$$

So the total strain energy for a single tetrahedral element is its rest volume $V_0$ multiplied by the energy density. For an undeformed state, this energy is zero. This naturally leads to the following constraint formulation: $$C({\mathbf{p}}_1,...,{\mathbf{p}}_4)=E_s=V_0{\Psi }_s(\mathbf{F}({\mathbf{p}}_1,...,{\mathbf{p}}_4))=0\text{\hspace{1em}(33.5.2)}$$

The PBD solver requires the gradient of the constraint with respect to the positions of the vertices. This gradient has a deep physical meaning, as it is directly related to the first Piola-Kirchhoff stress tensor, $\mathbf{P}$, which relates forces in the current configuration to areas in the reference configuration. The stress tensor is a function of the deformation gradient, $\mathbf{P}(\mathbf{F})$. Using the chain rule, the gradients of the energy for the first three vertices can be computed in a compact matrix form: $$\left[\frac{\partial E_s}{\partial {\mathbf{p}}_1},\frac{\partial E_s}{\partial {\mathbf{p}}_2},\frac{\partial E_s}{\partial {\mathbf{p}}_3}\right]=V_0\mathbf{P}(\mathbf{F}){\mathbf{D}}_m^{-T}\text{\hspace{1em}(33.5.3)}$$ The gradient for the fourth vertex is determined by the condition of translational invariance (i.e., the net force on the element must be zero): $$\frac{\partial E_s}{\partial {\mathbf{p}}_4}=-\sum _{i=1}^3\frac{\partial E_s}{\partial {\mathbf{p}}_i}\text{\hspace{1em}(33.5.4)}$$ This formulation extends naturally to surface meshes composed of triangles. For a triangle with rest area $A_0$, the constraint is $C(\mathbf{p})=A_0{\Psi }_s({\mathbf{F}}_{\text{tri}})$. The gradients are computed analogously: $$\left[\frac{\partial E_s}{\partial {\mathbf{p}}_1},\frac{\partial E_s}{\partial {\mathbf{p}}_2}\right]=A_0\mathbf{P}({\mathbf{F}}_{\text{tri}}){\mathbf{D}}_m^{-T}\text{and}\frac{\partial E_s}{\partial {\mathbf{p}}_3}=-\sum _{i=1}^2\frac{\partial E_s}{\partial {\mathbf{p}}_i}\text{\hspace{1em}(33.5.5)}$$

Remark 33.5.1 (Handling Element Inversion). Standard constitutive models are often not defined for degenerate or inverted elements (where $\det (\mathbf{F})\le 0$), which can cause instabilities in a simulation. This is a critical practical issue. This problem can be robustly addressed by augmenting the constitutive model with techniques designed for inversion handling.

The elegance of this strain energy constraint lies in its generality. It is not limited to a single material type, like the simple St. Venant-Kirchhoff model. More complex and physically accurate models, such as the Neo-Hookean model for rubber-like materials, can be readily supported by simply substituting the appropriate energy density function ${\Psi }_s$. This allows PBD to simulate complex physical effects like lateral contraction (Poisson's effect), anisotropy, and elastoplasticity with high fidelity and efficiency.

Direct Strain-Based Constraints

For applications requiring more granular or intuitive control over material properties, particularly for simulating anisotropy, it is advantageous to define separate constraints for each mode of deformation. This approach operates on the components of the right Cauchy-Green deformation tensor, $\mathbf{S}={\mathbf{F}}^T\mathbf{F}$.

Stretch Constraints

The diagonal entries of $\mathbf{S}$, $S_{ii}$, represent the squared stretch along the material's reference axes. In an undeformed state, $S_{ii}=1$. To ensure the constraint projection is linear and converges in a single step, the stretch constraint is formulated based on the stretch itself, rather than its square: $$C_{\text{stretch},i}(\mathbf{x})=\sqrt{S_{ii}}-1=0\text{\hspace{1em}(33.5.6)}$$ By defining one such constraint for each principal axis, one can assign independent stiffness parameters, which is the key to modeling anisotropic stretch resistance.

Shear Constraints

The off-diagonal entries of $\mathbf{S}$, $S_{ij}={\mathbf{f}}_i\cdot {\mathbf{f}}_j$ (where ${\mathbf{f}}_k$ are the column vectors of $\mathbf{F}$), measure the shear between material axes. A naive constraint $S_{ij}=0$ problematically couples shearing and stretching. To isolate the shearing deformation, the constraint must be normalized, thus depending only on the angle between the deformed axes: $$C_{\text{shear},ij}(\mathbf{x})=\frac{{\mathbf{f}}_i\cdot {\mathbf{f}}_j}{|{\mathbf{f}}_i||{\mathbf{f}}_j|}=0\text{\hspace{1em}(33.5.7)}$$ This formulation ensures that the shear constraints enforce orthogonality between the material axes without interfering with the stretch constraints that control their length. This decoupling is critical for intuitive and independent tuning of material behavior.