The Most Performant Bending Constraint | XPBD

XPBD uses particle constraints to ensure that particles move in a certain manner, relative to other particles. Let's take a look at a bending constraint that will be super performant to create smooth and visually pleasing visualizations.

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

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 a few months, you probably know of my goal to simulate cloth 🧵. For the past few weeks, we’ve been talking about the distance constraint and the isometric bending constraint. This week we will be focusing on a more performant bending constraint.

Why change course? Because the isometric bending constraint is expensive to solve 🥵. Matthias Müller shows us an alternative to simulate bending with great results !

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! 💻

The best performing bending constraint

In this presentation by Matthias Müller, Matthias shows us an alternative to integrate a bending constraint without making complex and expensive calculations. Fortunately, it is very similar to the distance constraint that we already talked about a few weeks ago in this article.

Here is an image that shows exactly what we’re trying to accomplish:

More performant bending constraint

As you can see, there is just one distance constraint between two particles of two adjacent triangles. Of course, this method is a little more complex than a normal distant constraint because you first have to find what triangles are adjacent to each other.

Afterwards, we follow the same process we followed to satisfy the distance constraint. 😮‍💨

We can solve the equations from the last section in order to find the final position corrections, Δxn\Delta \mathbf{x}_n:

C(x2,x3)=∣x2,3∣−dn=x2,3∣x2,3∣Δx2C(x2,x3)=nΔx3C(x2,x3)=−n∴Δx2=−w2w2+w3(∣x2,3∣−d)n∴Δx3=+w3w2+w3(∣x2,3∣−d)n C(\mathbf{x_2}, \mathbf{x_3}) = |\mathbf{x}_{2,3}| - d\\[6pt] \mathbf{n} = \frac{\mathbf{x}_{2,3}}{|\mathbf{x}_{2,3}|}\\[6pt] \Delta_{\mathbf{x_2}}C(\mathbf{x_2}, \mathbf{x_3}) = \mathbf{n}\\[6pt] \Delta_{\mathbf{x_3}}C(\mathbf{x_2}, \mathbf{x_3}) = -\mathbf{n}\\[6pt] \therefore \Delta \mathbf{x_2} = - \frac{w_2}{w_2 + w_3}(|\mathbf{x}_{2,3}| - d) \mathbf{n}\\[6pt] \therefore \Delta \mathbf{x_3} = + \frac{w_3}{w_2 + w_3}(|\mathbf{x}_{2,3}| - d) \mathbf{n}\\[6pt]

The code

Fortunately, we will see that the code for this constraint is almost similar to the distance constraint. Before seeing it, we must first understand how to find two adjacent triangles.

How to find adjacent triangles

Let’s imagine that we load in a mesh to play with that consists of 1,000 triangles. How could we figure which triangles are adjacent to each other? Fortunately, there is a little trick :

Adjacent edges

As you can see, there are two triangles that share an edge (2 and 5). We can exploit this fact! Our steps are as follows:

  1. Create an array where we can search for shared edges between two triangles.
  2. Use this array to find the four points that form two triangles, as in the picture above.

We first create a list that consists of the edges with ids global that are between each pair of particles. It’s our job to label each edge.

For example, imagine we create a list of edges from this image:

Adjacent edges

p_id0 p_id1 ei
0 1 2
1 2 0
2 0 1
0 3 3
3 1 4
1 0 5

You can probably see that the particles p0p_0 and p1p_1 share global edges (2) and (5). To find this pair in the code, we sort the edge array in ascending order by the particle id number in each pair:

p_id0 p_id1 ei
0 1 2
0 1 5
0 2 1
0 3 3
1 2 0
1 3 4

Next, we create an array of neighbors, or neighbors in English, which is initialized with -1s, where -1 means ‘no neighbor’. If we create an array of neighbors correctly, we would see an array of neighbors as follows:

[ne0,ne1,ne2,ne3,ne4]→[−1,−1,5,−1,−1,2][n_{e0}, n_{e1}, n_{e2}, n_{e3}, n_{e4}] \rightarrow [-1, -1, 5, -1, -1, 2]

Here is the code for a function, findTriNeighbors, which will helps us create this neighbor list:

If you have followed my posts up until now, you will know that I only use WebGPU and JavaScript, without a library like ThreeJS. I say this because you will need to access the index buffer that we have already discussed in this article.

This list of neighbors will be essential for the next step.

With our neighbor array, we have to create a list containing the identifier of each particle in the quartet that forms each pair of triangles. It will be used to place a distance constraint between opposing particles as in the image below:

More performant bending constraint

The process of doing so is shown in the image below:

Method to find tri pair ids

We visit each face of our mesh to determine if there is a shared edge (by searching the list of neighboring edges, neighbors). If so, we save the three vertices of this triangle and then we visit the adjacent triangle - where we save the fourth vertex in the configuration.

Here is the getTriPairIds function which does exactly that:

With our list of triPairIds, we are ready to move on to satisfying the constraint! 🎉

The performant bending constraint code

So as you can see, this code is almost the same as the distance constraint code. Just use the two vertices that are not shared instead of just any two adjacent vertices when calculating final position corrections!


Comments for The Most Performant 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!