Physics-Based Simulation

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.

Theory: Dihedral Angle Constraint Formulation

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$.

Definition 33.2.1 (Dihedral Angle Bending Constraint). The bending constraint function $C_{\text{bend}}$ for four vertices forming two adjacent triangles is: $$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.2.1)}$$ 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.

Alternative: Isometric Bending Model

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.2.2)}$$ 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.

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.2.3)}$$

This model is particularly effective for garment simulation where fabric is expected to deform isometrically (i.e., without stretching).

Implementation: Bending Constraint Solver

@ti.kernel
def solve_bending_constraints(compliance: ti.f64, dt: ti.f64):
    """Solve bending constraints using distance-based approach"""
    alpha = compliance / (dt * dt)
    for i in range(num_bending_constraints):
        id0, id1 = bending_ids[i, 0], bending_ids[i, 1]
        w0, w1 = inv_mass[id0], inv_mass[id1]
        w_sum = w0 + w1
        if w_sum == 0.0: continue
        
        p0, p1 = pos[id0], pos[id1]
        delta = p0 - p1
        dist = delta.norm()
        if dist == 0.0: continue
        
        grad = delta / dist
        C = dist - bending_lengths[i]
        s = -C / (w_sum + alpha)
        
        pos[id0] += s * w0 * grad
        pos[id1] -= s * w1 * grad

This implementation uses a simplified distance-based approach where bending constraints are treated as distance constraints between non-adjacent vertices. While this doesn't capture the full dihedral angle behavior, it provides a computationally efficient approximation that works well for many cloth simulation scenarios.

We need rest lengths to solve the constraints, so we compute all rest bending lengths once before the simulation:

@ti.kernel
def init_physics():
    """Initialize physics state including rest lengths for all constraints"""
    # Initialize stretching constraint rest lengths
    for i in range(num_stretching_constraints):
        id0, id1 = stretching_ids[i, 0], stretching_ids[i, 1]
        stretching_lengths[i] = (pos[id0] - pos[id1]).norm()
    
    # Initialize bending constraint rest lengths
    for i in range(num_bending_constraints):
        id0, id1 = bending_ids[i, 0], bending_ids[i, 1]
        bending_lengths[i] = (pos[id0] - pos[id1]).norm()

To identify which vertices become bending constraints, we find pairs of triangles that share an edge, then create constraints between the vertices that are not part of the shared edge:

def find_constraint_indices(tri_ids_np):
    """Extract both stretching and bending constraints from triangle mesh"""
    edge_to_tri_map = {}
    for i, tri in enumerate(tri_ids_np):
        for j in range(3):
            v0_idx, v1_idx = tri[j], tri[(j + 1) % 3]
            edge = tuple(sorted((v0_idx, v1_idx)))
            if edge not in edge_to_tri_map:
                edge_to_tri_map[edge] = []
            edge_to_tri_map[edge].append(i)
    
    stretching_ids = list(edge_to_tri_map.keys())
    bending_ids = []
    
    # Find bending constraints (non-adjacent vertices in adjacent triangles)
    for edge, tris in edge_to_tri_map.items():
        if len(tris) == 2:
            tri0_idx, tri1_idx = tris[0], tris[1]
            p2 = [v for v in tri_ids_np[tri0_idx] if v not in edge][0]
            p3 = [v for v in tri_ids_np[tri1_idx] if v not in edge][0]
            bending_ids.append([p2, p3])
    
    return np.array(stretching_ids, dtype=np.int32), np.array(bending_ids, dtype=np.int32)

This creates a network of constraints that resist bending deformation by connecting non-adjacent vertices in adjacent triangles.

stretching_compliance = 0.0  # Very stiff stretching (nearly inextensible)
bending_compliance = 1.0      # Moderate bending resistance

Lower stretching compliance makes the cloth more resistant to stretching, while higher bending compliance allows more bending flexibility, creating softer, more drapeable cloth.