The Isometric Bending Constraint | XPBD

XPBD uses particle constraints to ensure that particles move in a certain manner, relative to other particles. Let's take a deeper look into a common constrant that we'll find in a cloth simulation, the isometric bending constraint.

Keywords: position based dynamics, cloth simulation, lagrange, physics, computer graphics

By Carmen Cincotti  

A big piece of the XPBD puzzle, which we already talked about a few weeks ago, is the satisfaction of particle constraints. Particle constraints force a particle to move in a certain way, relative to the other particles in the system.

How a particle moves back to the string

If you’ve been following me for the past few months, you most likely know of my goal to simulate fabric 🧵. Last week, we talked about distance constraints. This week, we’ll be focusing on isometric bending constraints.

Here are once again those famous articles that we already studied together during an article from a few weeks ago on XPBD that I will reference for the math during this post:

A summary of XPBD

Recall that the XPBD algorithm with small steps is defined as follows:

Algorithm of XPBD with small steps

⚠️ As I said, we have already covered this algorithm thoroughly, and I recommend that you read the article if you are already lost.

For now, we are only interested in the loop that is in lines (6) - (11). As you can see, the function of this loop is to satisfy the constraints of the simulation in order to directly correct the positions of each particle that participates in the constraints.

We just covered the reason why this simulation method is called as Position Based Dynamics - PBD 🎉.

Satisfying the constraints

Our objective during the loop is to calculate the correction vector that a particle must move in order to satisfy, or come closer to satisfying, the constraints.

Δx=ΔλM−1∇CxT\Delta \mathbf{x} = \Delta \lambda \mathbf{M}^{-1} \nabla C \mathbf{x}^T

where Δλ\Delta \lambda is defined as:

Δλ=−C(x)∇C(x)M−1∇C(x)T+α~\Delta \lambda = \frac{-C(\mathbf{x})}{\nabla C(\mathbf{x}) \mathbf{M}^{-1} \nabla C(\mathbf{x})^T + \tilde{\alpha}}

Note that α~=αΔt2\tilde{\alpha} = \frac{\alpha}{\Delta t^2} where α\alpha is the constraint compliance. Recall that α=1k\alpha = \frac{1}{k} where kk is the stiffness of the constraint. C(x)C(\mathbf{x}) is the constraint equation, ∇C(x)\nabla C(\mathbf{x}) is the gradient of ∇C(x)\nabla C(\mathbf{x}) and is therefore a vector. M−1\mathbf{M}^{-1} is a diagonal matrix of inverse masses of the particles participating in the constraint.

Δλ\Delta \lambda, C(x)C(\mathbf{x}), and α~\tilde{\alpha} are scalar values.

Using these equations, we will simulate some constraints with code! 💻

Isometric bending constraint

Let’s now talk about the isometric bending constraint. Here is an example of what that might look like during a cloth simulation:

Bending cloth example

The isometric bending constraint is a bending constraint for inextensible surfaces. Basically, it’s a great bending model for cloth simulations due to the fact that cloth cannot normally be stretched (although that depends on the material).

The constraint consists of four particles that form a configuration of two triangles, which looks like the following:

Four triangles configuration

where each purple circle is a particle at a location xn\mathbf{x}_n. Edges of the two triangles are labeled as ene_n. This configuration is called a stencil.

To calculate the necessary displacement of a particle in order to satisfy the constraint - we first start with the isometric bending constraint equation:

Cbend(xs)=12∑i,jQi,jxiTxj C_{bend}(\mathbf{x}_s) = \frac{1}{2}\sum_{i,j}\mathbf{Q}_{i,j}\mathbf{x}_{i}^T\mathbf{x}_{j}

What 😨?! Don’t worry, we’ll break this equation down little by little.

We start with the matrix Qi,j∈R4x4\mathbf{Q}_{i,j} \in \mathbb{R}^{\text{4x4}}, which is the local constant Hessian bending energy. This resource gives the derivation.

We initialize the matrix Q\mathbf{Q} only once for each configuration of four particles. Recall that the Hessian matrix includes the second order partial derivatives:

Qi,j(f)=∂2f∂xi∂xj \mathbf{Q}_{i,j}(f) = \frac{\partial^2 f}{\partial{x}_i\partial{x}_j}

where A0A_0 and A1A_1 are areas of adjacent triangles. Recall that we can find the area of ​​a triangle if the locations of it’s three points are given with a small linear algebra trick:

Finding area of triangle with cross product

A=12∣AB⃗×AC⃗∣A = \frac{1}{2}|\vec{AB} \times \vec{AC}|

K\mathbf{K} is a vector which consists of the cotangents between the triangle edges cjk=cot⁡∠ej,ekc_{jk} = \cot\angle e_j,e_k:

K=[c01+c04,c02+c03,−c01−c02,−c03−c04] \mathbf{K} = [c_{01} + c_{04}, c_{02} + c_{03}, -c_{01} - c_{02}, -c_{03} - c_{04}]

How to find the cotangent between two vectors 🤔

The cotangent is defined as follows:

cot⁡θ=cos⁡θsin⁡θ \cot \theta = \frac{\cos \theta}{\sin \theta}

We can find cos⁡θ\cos \theta with the dot product:

a⋅b=∣a∣×∣b∣×cos⁡θ∴cos⁡θ=a⋅b∣a∣×∣b∣ \mathbf{a} \cdot \mathbf{b} = |\mathbf{a}| \times |\mathbf{b}| \times \cos \theta \\[6pt] \therefore \cos \theta = \frac{\mathbf{a} \cdot \mathbf{b}}{|\mathbf{a}| \times |\mathbf{b}|}

To find sin⁡θ\sin \theta, we use the cross product:

∣a×b∣=∣a∣∣b∣sin⁡θ∴sin⁡θ=∣a×b∣∣a∣×∣b∣ |\mathbf{a} \times \mathbf{b}| = |\mathbf{a}| |\mathbf{b}| \sin \theta \\[6pt] \therefore \sin \theta = \frac{|\mathbf{a} \times \mathbf{b}|}{|\mathbf{a}| \times |\mathbf{b}|}

We can therefore rewrite the equation for cot⁡θ\cot \theta as follows:

cot⁡θ=cos⁡θsin⁡θ=a⋅b∣a∣×∣b∣∣a×b∣∣a∣×∣b∣=a⋅b∣a×b∣ \cot \theta = \frac{\cos \theta}{\sin \theta} \\[6pt] = \frac{\frac{\mathbf{a} \cdot \mathbf{b}}{|\mathbf{a}| \times |\mathbf{b}|}}{\frac{|\mathbf{a} \times \mathbf{b}|}{|\mathbf{a}| \times |\mathbf{b}|}} \\[6pt] = \frac{\mathbf{a} \cdot \mathbf{b}}{|\mathbf{a} \times \mathbf{b}|}

Let us now go back to this equation:

Cbend(xs)=12∑i,jQi,jxiTxj C_{bend}(\mathbf{x}_s) = \frac{1}{2}\sum_{i,j}\mathbf{Q}_{i,j}\mathbf{x}_{i}^T\mathbf{x}_{j}

x\mathbf{x} is the location of the particles. Finally, we need the gradient ∇Cbend\nabla C_{bend} to find Δλ\Delta \lambda which is defined as follows:

∇Cbend=∑jQi,jxj \nabla C_{bend} = \sum_{j}\mathbf{Q}_{i,j}\mathbf{x}_{j}

We now have all the equations that we need. I find that these equations are difficult to truly grasp without seeing them at work in the code. So, let’s look at some code now!

The code

If you have already seen the code for the distance constraint from the last article, you will notice that the simulation loop does not change. We are just concerned with the bendingConstraint function.

Fortunately, the code is not rocket science. First we need to calculate C, then the gradient - grad, which is ∇C\nabla C . Then we find deltaLagrangianMultiplier which is Δλ\Delta \lambda. Finally, we find Δx\Delta \mathbf{x}, which is the displacement of a particle during a substep in time.

When running the program - we should see the following logged to the console:

BEFORE APPLYING XPBD LOOP ------------------------- POSITION P1: [ '0.0000', '0.0000', '0.0000' ] POSITION P2: [ '0.0000', '1.0000', '0.0000' ] POSITION P3: [ '-0.5000', '0.5000', '0.0000' ] POSITION P4: [ '0.5000', '0.5000', '0.0000' ] VELOCITY P1: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P2: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P3: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P4: [ '0.0000', '0.0000', '0.0000' ] AFTER APPLYING XPBD LOOP ------------------------ POSITION P1: [ '0.0000', '0.0000', '0.0000' ] POSITION P2: [ '0.0000', '1.0000', '0.0000' ] POSITION P3: [ '-0.5000', '0.5000', '0.0000' ] POSITION P4: [ '0.5000', '0.5000', '0.0000' ] VELOCITY P1: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P2: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P3: [ '0.0000', '0.0000', '0.0000' ] VELOCITY P4: [ '0.0000', '0.0000', '0.0000' ]

Great! All that remains is to integrate these constraints into the fabric simulation! See you then!


Comments for The Isometric Bending Constraint | XPBD

Written by Carmen Cincotti, computer graphics enthusiast, language learner, and improv actor currently living in San Francisco, CA.  Follow @CarmenCincotti


Interested in contributing to Carmen's Graphics Blog? Click here for details!