Physics-Based Simulation

Cloth: Stretching and Bending

In the context of Position-Based Dynamics (PBD), the complex mechanical behaviors of cloth, such as its resistance to stretching and bending, are modeled through a system of geometric constraints. Instead of accumulating forces, the PBD framework directly manipulates the positions of the mesh vertices to satisfy these constraints in an iterative manner. This section will detail the formulation of fundamental constraints for cloth simulation.

Stretching Resistance via Distance Constraints

The primary characteristic of most textiles is their high resistance to stretching. In PBD, this property is enforced by constraining the distance between connected particles to remain close to its initial, or rest, distance. This is one of the simplest yet most crucial constraints in the PBD ecosystem (we have already seen this constraint in Example 32.4.1)

Example 33.1.1 (Stretching Constraint). Consider two particles, $i=1,2$, with positions ${\mathbf{p}}_1$ and ${\mathbf{p}}_2$, masses $m_1$ and $m_2$, and a rest distance $d$ between them. The stretching constraint function $C$ is defined as the difference between the current distance and the rest distance: $$C({\mathbf{p}}_1,{\mathbf{p}}_2)=\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert -d\text{\hspace{1em}(33.1.1)}$$ The goal is to find corrections $\Delta {\mathbf{p}}_1$ and $\Delta {\mathbf{p}}_2$ such that $C({\mathbf{p}}_1+\Delta {\mathbf{p}}_1,{\mathbf{p}}_2+\Delta {\mathbf{p}}_2)=0$. The gradients of the constraint function with respect to the particle positions are: $${\nabla }_{{\mathbf{p}}_1}C=\frac{{\mathbf{p}}_1-{\mathbf{p}}_2}{\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert }=\mathbf{n}\text{and}{\nabla }_{{\mathbf{p}}_2}C=-\frac{{\mathbf{p}}_1-{\mathbf{p}}_2}{\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert }=-\mathbf{n}$$ where $\mathbf{n}$ is the unit vector along the axis connecting the two particles. Following the general PBD projection formula (32.4.4), the scalar Lagrange multiplier $\lambda$ is computed as: $$\lambda =\frac{C({\mathbf{p}}_1,{\mathbf{p}}_2)}{w_1\Vert {\nabla }_{{\mathbf{p}}_1}C{\Vert }^2+w_2\Vert {\nabla }_{{\mathbf{p}}_2}C{\Vert }^2}=\frac{\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert -d}{w_1+w_2}$$ where $w_i=1/m_i$ is the inverse mass of particle $i$. The position corrections are then found by moving the particles along their respective gradient directions, scaled by their inverse mass and $\lambda$: $$\begin{gathered} \Delta {\mathbf{p}}_1=-\frac{w_1}{w_1+w_2}(\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert -d)\mathbf{n} \\ \Delta {\mathbf{p}}_2=+\frac{w_2}{w_1+w_2}(\Vert {\mathbf{p}}_1-{\mathbf{p}}_2\Vert -d)\mathbf{n} \end{gathered}$$ These corrections, when applied, will move the particles to exactly satisfy the rest length. Note that the total correction is distributed between the particles based on their inverse mass, ensuring that lighter particles move more than heavier ones and that linear momentum is conserved ($\sum m_i\Delta {\mathbf{p}}_i=0$).

To achieve different levels of elasticity, a stiffness parameter $k\in [0,1]$ can be introduced by scaling the corrections $\Delta {\mathbf{p}}_i$ by $k$. This allows for materials with varying elasticity, from perfectly rigid ($k=1$) to completely unstiff ($k=0$).

Dihedral Angle Bending Constraints

While stretching constraints maintain the structural integrity of the cloth mesh, they do not prevent it from folding unnaturally. Bending resistance, which dictates how the cloth wrinkles and drapes, is modeled by constraining the angle between adjacent triangles.

The constraint is defined for a pair of triangles $({\mathbf{p}}_1,{\mathbf{p}}_3,{\mathbf{p}}_2)$ and $({\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_4)$ sharing a common edge $({\mathbf{p}}_1,{\mathbf{p}}_2)$. The bending resistance is a function of the dihedral angle $\varphi$ between the two triangles, which is the initial angle between their respective normal vectors ${\mathbf{n}}_1$ and ${\mathbf{n}}_2$. The constraint aims to restore this angle to its rest value, ${\varphi }_0$.

The constraint function is formulated as: $$C_{\text{bend}}({\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3,{\mathbf{p}}_4)=\arccos ({\mathbf{n}}_1\cdot {\mathbf{n}}_2)-{\varphi }_0\text{\hspace{1em}(33.1.2)}$$ where the normals are computed as: $$\begin{array}{rl}{\mathbf{n}}_1& =({\mathbf{p}}_{2,1}\times {\mathbf{p}}_{3,1})/\Vert {\mathbf{p}}_{2,1}\times {\mathbf{p}}_{3,1}\Vert \\ {\mathbf{n}}_2& =({\mathbf{p}}_{2,1}\times {\mathbf{p}}_{4,1})/\Vert {\mathbf{p}}_{2,1}\times {\mathbf{p}}_{4,1}\Vert .\end{array}$$ The gradients of this function with respect to the four vertex positions $({\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3,{\mathbf{p}}_4)$ are then computed, and the standard PBD projection mechanism is used to derive the position corrections. The stiffness of bending is determined using $k_{bend}$ parameter.

A significant advantage of this formulation is its independence from stretching. Because the angle is defined by normalized vectors, the constraint is invariant to the lengths of the triangle edges.

Isometric Bending

For surfaces that are nearly inextensible, the isometric bending model [Bergou et al. 2006] can be used. This model provides a robust formulation based on the local Hessian of the bending energy.

This model considers a stencil for each interior edge ${\mathbf{e}}_0$ of the mesh, consisting of the four vertices of the two triangles adjacent to that edge, labeled ${\mathbf{p}}_0,{\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3$. The local bending energy for this stencil is defined as a quadratic form: $$E_{\text{bend}}({\mathbf{p}}_s)=\frac{1}{2}\sum _{i,j\in \{0,1,2,3\}}{Q}_{ij}({\mathbf{p}}_i^{\top }{\mathbf{p}}_j)\text{\hspace{1em}(33.1.3)}$$ where ${\mathbf{p}}_s=({\mathbf{p}}_0,{\mathbf{p}}_1,{\mathbf{p}}_2,{\mathbf{p}}_3{)}^T$ is the vector of stencil positions and $\mathbf{Q}\in {\mathbb{R}}^{4\times 4}$ is a constant matrix representing the local Hessian of the bending energy. This matrix depends only on the rest geometry of the stencil and can be precomputed. Its entries are derived from the cotangents of the angles within the two triangles.

The bending constraint is defined directly from this energy: $C_{\text{bend}}({\mathbf{p}}_s)=E_{\text{bend}}({\mathbf{p}}_s)$. Since the energy is quadratic in the positions, its gradient is linear and straightforward to compute: $$\frac{\partial C_{\text{bend}}}{\partial {\mathbf{p}}_i}=\sum _{j\in \{0,1,2,3\}}{\mathbf{Q}}_{ij}{\mathbf{p}}_j\text{\hspace{1em}(33.1.4)}$$ This model is particularly effective for garment simulation where fabric is expected to deform isometrically (i.e., without stretching).