A quasi-static solver means that the solver is almost static. In this case that means the solver is independent of time but deforms a geometry to satisfy some constraints. This can be beneficial if no dynamics are required, for example in modeling and sculpting applications where a modeler doesn’t care about deformation over time, or for post-simulation sculpting. This post is about how one can implement a quasi-static soft body solver in Houdini using OpenCL and position based dynamics (PBD).
From left: rest geometry(blue), animated geometry(red), solved geometry(green) at 500 constraint iterations and 1^10 stiffness.
Introduction
I did this setup over a year ago, and I struggled quite hard before I finally got it working. I had previously dabbled with cloth PBD so I felt familiar with the idea, but since this method is much more math heavy it took me quite some time to understand the concept. And because I’m cursed with bad memory I don’t really remember much from when I first set it up. So I figured, let’s write it down and perhaps it will stick better, and perhaps someone else can benefit from this as well!
The idea is to implement the Strain Based PBD paper for soft bodies. You should definitely read the paper for a better explanation of how it works from a math perspective.
Basic PBD
PBD is short for Position Based Dynamics and is a simple idea of how we can simulate physical properties of a geometry. The method involves creating “constraints” which will define how the geometry is allowed to deform. There has been lots of PBD constraints developed by many people, everything from cloth and hair to tissue and fluids. The most basic one is probably the “edge” constraint. Consider a line segment with rest length alpha, if the line segment stretches to a new length then we should try to shorten it until the length matches alpha again. This relationship between the rest length alpha and the current length is the constraint, and we solve it by updating the points x1 and x2 of the line segment. Here’s how that can look:
C(x1,x2) = |x1 −x2|−\alpha
And then the point deltas become:
\Delta p_1 = - \frac {w_1}{w_1+w_2} C(x_1,x_2) \frac{x_1-x_2}{|x_1-x_2|}\Delta p_2 = + \frac {w_2}{w_1+w_2} C(x_1,x_2) \frac{x_1-x_2}{|x_1-x_2|}Where w1 and w2 are the inverse masses of the particles.
If we had two line segments with points x1, x2, x3, where x2 connect the two segments, and we first solve for seg1, x2 will deform, then we solve for seg2 and x2 will deform again, meaning both lines pull on x2, so only the latter will satisfy the constraint. To get around this issue, we can just solve the constraints for the two lines multiple times, and they will eventually satisfy their respective constraints. The general idea of PBD is therefore: solve a bunch of constraints iteratively over and over again until we’re happy.
PBD is a great method for simulating deformable materials since the solver is generally quite easy to implement and it can easily run on GPUs. If you didn’t know, Vellum in Houdini is purely PBD and a lot of constraints are implemented for it (take a look in <your_houdini_install>/houdini/ocl/sim/pbd_constraints.cl for some inspiration). For a deeper understanding I suggest reading the original PBD publication.
Strain Based PBD
Strain based constraints are independent off edge directions, instead the constraint is based on material coordinates. In practice, that means we create a single matrix for each tetrahedra or triangle which inscribes the local deformation, and the constraint will then try to limit the deformation to fit its rest configuration (similar to the example above, but solve for each axis independently). This effectively allows us to simulate anisotropic materials, since we can have a separate stiffness for each axis of the matrix. In this case, we also simulate the shear separate from the stretch, which means we have even more control of the material.
Implementation
Here’s what we’re aiming for:
- Be able to deform any(within reason) tetrahedral mesh and get valid result
- Control over anisotropic rest length
- Control over anisotropic stretch
- Control over anisotropic shear
- Volume preservation
- Fit the XPBD framework
For this, we need two inputs, our rest geometry defining the rest configuration of each tetrahedra, and the animated geometry. It is important that they have the same topology.
Initialize Precomputed Variables
Let’s start by creating some precomputed variables. Namely the Qinv material matrix attribute explained in eq. 6 of the paper, the stiffness parameters, the rest volume and some other parameters we’ll need. This all goes in a wrangle on the non-animated input running over Primitives:
matrix3 computeTetMaterialPosition(const int pts[]; float volume) {
vector x0 = point(0,"P",pts[0]);
vector x1 = point(0,"P",pts[1]);
vector x2 = point(0,"P",pts[2]);
vector x3 = point(0,"P",pts[3]);
vector x = x1 - x0;
vector y = x2 - x0;
vector z = x3 - x0;
// Compute material matrix Q
matrix3 Q = set(x, y, z);
// Compute rest volume
volume = dot(z, cross(x, y));
// Since we only care about Q^-1, invert it here
return invert(Q);
}
// Used by OpenCL in order to know what points belong to each tet
i[]@__pts = primpoints(0,@primnum);
// The Q^-1 material matrix
3@__Qinv = computeTetMaterialPosition(i[]@__pts, f@__volume);
// Lagrange multiplier for each axis of Q and for volume constraint
resize(f[]@__L, 10);
// Set stiffness
if (chi("uniform_stretch")) {
v@__stretch_resist = chf("stretchf");
} else {
v@__stretch_resist = chv("stretchv");
}
if (chi("uniform_shear")) {
v@__shear_resist = chf("shearf");
} else {
v@__shear_resist = chv("shearv");
}
f@__volume_resist = chf("volume");
// Set restlength
if (chi("uniform_restlength")) {
v@__rls = chf("restlengthf");
} else {
v@__rls = chv("restlengthv");
}
f@__volume *= dot(v@__rls, v@__rls) / 3.0f;
On this same thread, we can also place a graphcolor node which will be used to set the OpenCL workgroups.
While on the topic, let’s also set some flags that can be sent to OpenCL. In this case, we just want some toggles for enabling/disabling the constraints. In a new detail wrangle:
s@__flags;
// Enable tet constraints
if (chi("solve_tets")) {
s@__flags += "-DTET_CONSTRAINTS ";
}
// Enable volume constraints
if (chi("solve_volume")) {
s@__flags += "-DVOL_CONSTRAINTS ";
}Now all constraint values should be initialized. We can copy them over to our first input using an attribcopy node, copying all of our __* primitive and detail attributes.
On the first input thread, we can initialize the mass for our points. This can be driven by a mask just as well, but for now I set all the point masses to be 1. For any points outside the point group __group, it will default to 0:
if (inprimgroup(0, "__group", @primnum)) {
f@__mass = 1; // Particle mass
}Setting up OpenCL
We can now use these attributes in an OpenCL sop node to deform the geometry. Houdini 20 differs quite a bit when it comes to the OpenCL sop node compared to earlier Houdini versions. I will use Houdini 20 but without the @-bindings and other new stuff.
Options
Firstly we can specify the workset attributes we created from the graph coloring which will tell OpenCL which tetrahedra will run in the same local work group.
Bindings
Then we can set our bindings. In this case, we only need P and L to be writeable, all the others can be left on only readable. Also, make sure to set the correct size of your vectors:
Kernel
Then on the Kernel tab, we can specify the __flags attribute we created before.
Generated Code
Lastly, we can click Generate Kernel in the Generated Code tab to get a default kernel function definition with all our parameters created for us. The result will look something like:
kernel void QuasiStaticStretch(
int worksets_begin,
int worksets_length,
float timeinc,
int P_length,
global fpreal * restrict P ,
int pts_length,
global int * restrict pts_index,
global int * restrict pts ,
int stretch_resist_length,
global float * restrict stretch_resist ,
int shear_resist_length,
global float * restrict shear_resist ,
int volume_resist_length,
global float * restrict volume_resist ,
int L_length,
global int * restrict L_index,
global float * restrict L ,
int M_length,
global float * restrict M ,
int restlength_length,
global float * restrict restlength ,
int rest_volume_length,
global float * restrict rest_volume ,
int Qinv_length,
global float * restrict Qinv
)
{
int idx = get_global_id(0);
if (idx >= worksets_length)
return;
idx += worksets_begin;
}This can be copied in to the Kernel Code parameter and everything should be running fine without errors.
The Kernel
We’ll want two functions for our two constraints, the stretch/shear and the volume constraints. We can start with the easier one just to get comfortable 🙂
Volume Constraint Solve
Just as we did when setting up the rest volume, we need to calculate the deformed volume for a tetrahedra
V = \frac{(p_1-p_0) \times (p_2-p_0) \cdot (p_3-p_0)}{6}where p0..3 are the points of the tetrahedra. Formally, this is the triple product of 3 edges of the tetrahedra, or the determinant of the matrix P, divided by 6, which forms the volume function for a tetrahedra.
Following the example from before, the constraint formulation is then:
C(p_0, p_1, p_2, p_3) = V - V_r
where Vr is the rest volume. We can also find the derivatives of each point:
\begin{align*}
\nabla p_0 = (p2-p1) \times (p2-p3) \\
\nabla p_1 = (p2-p0) \times (p3-p0) \\
\nabla p_2 = (p0-p1) \times (p3-p1) \\
\nabla p_3 = (p1-p0) \times (p2-p0) \\
\end{align*}We can then define the deltas for each point:
\begin{align*}
\Delta p_k = -C(p_0, p_1, p_2, p_3)\frac{w_k}{w_{sum}}\nabla p_k \\
\end{align*}wsum is the sum of the inverse masses of the points. However, In this case, we can scale it by the squared norm of the point derivatives:
w_{sum} = \sum_{k}w_k|\nabla p_k|^2This is functionally fine, however, we’ve left the stiffness parameter alpha out. One could multiply the point deltas by alpha, as long as alpha is within the range [0,1]. However, a more modern implementation of PBD makes use of the XPBD framework, which is a much more stable way to modify the stiffness parameters. The issue with regular PBD is that the stiffness is dependent on the number of iterations we use to solve the constraint, whereas XPBD respects the result of any previous solve iterations, and modifies the strength of the constraint thereafter. See figure 6 of the XPBD paper for a visual demonstration.
The basic idea is to use the inverse stiffness coefficient, and modify how much we’ll displace each point for every solve iteration, rather than doing the naive approach demonstrated above. This is achieved using a Lagrange multiplier. We will reformulate what we’ve done slightly, in order to fit this new method (I took inspiration from Vellum on this part, since the XPBD paper goes a bit over my head). This is the final function for solving our volume constraints:
#define EPS 1e-10f
#include <matrix.h>
static fpreal safediv(fpreal a, fpreal b) {
return select(0.0f, a / b, b > 0.0f);
}
void solveVolumeTet(const float timeinc,
const int ptidx,
const int prmidx,
global int* pts,
global fpreal* P,
global fpreal* tmp,
const float rest_volume,
const float stiffness,
global float* L,
global float* M) {
int pt0 = pts[ptidx];
int pt1 = pts[ptidx + 1];
int pt2 = pts[ptidx + 2];
int pt3 = pts[ptidx + 3];
// Invert stiffness for XPBD
fpreal alpha = safediv(1.0f, stiffness);
if (alpha == 0.0f) {
return;
}
// Read result from the shear constraint solve.
fpreal3 p0 = vload3(pt0, P);
fpreal3 p1 = vload3(pt1, P);
fpreal3 p2 = vload3(pt2, P);
fpreal3 p3 = vload3(pt3, P);
// Invert point masses
fpreal w0 = safediv(1.0f, M[pt0]);
fpreal w1 = safediv(1.0f, M[pt1]);
fpreal w2 = safediv(1.0f, M[pt2]);
fpreal w3 = safediv(1.0f, M[pt3]);
// Current volume
fpreal volume = dot(p3-p0, cross(p1-p0, p2-p0));
// Our point derivates
fpreal3 dp0 = cross(p2-p1, p2-p3);
fpreal3 dp1 = cross(p2-p0, p3-p0);
fpreal3 dp2 = cross(p0-p1, p3-p1);
fpreal3 dp3 = cross(p1-p0, p2-p0);
// Calculate Wsum
fpreal wsum = w0 * dot(dp0, dp0) +
w1 * dot(dp1, dp1) +
w2 * dot(dp2, dp2) +
w3 * dot(dp3, dp3);
if (wsum < EPS) {
return;
}
// l is the lagrange multiplier of this constraint
fpreal l = L[prmidx * 10 + 9];
// Compute C just as before
fpreal C = (1.0f / 6.0f) * (volume - rest_volume);
// Not strictly necessary, but seem to give a bit more stable result
alpha /= timeinc * timeinc;
// dL is the derivate of the constraint for this iteration
fpreal dL = (-C - alpha * l) / (wsum + alpha);
// Add dL back on to L where the next iteration will read our result
L[prmidx * 10 + 9] += dL;
// Modify our points.
p0 += dL * w0 * dp0;
p1 += dL * w1 * dp1;
p2 += dL * w2 * dp2;
p3 += dL * w3 * dp3;
// Store points back in to the temp array
vstore3(p0, pt0, P);
vstore3(p1, pt1, P);
vstore3(p2, pt2, P);
vstore3(p3, pt3, P);
}We can then update the kernel function to call solveVolumeTet for each tetrahedra:
#ifdef VOL_CONSTRAINTS
solveVolumeTet(timeinc, pts_index[idx], idx, pts, P, rest_volume[idx], volume_resist[idx], L, M);
#endif
As demonstrated here with a simple tube scaled to 0.5 in y, the constrained geometry(green) will bulge out in order to maintain its volume.
Stretch / Shear Constraint Solve
Without the necessary foundation in mathematics I had a hard time fully grasping this constraint at first, and I’m still not sure I fully do. I will try to explain it to the best of my ability, though!
With that said, the idea of this constraint is to solve the stretch and shear of a tetrahedra. By solving them separately, we have a lot of control over the deforming material. We previously created the Qinv matrix inside the wrangle. We’ll use it to calculate the deformation matrix F
F = PQ^{-1}P is calculated just as Q but from the animated geometry, i.e from the material coordinates of the tetrahedrons. These three matrices are explained in eq. 5, 6, 7 of the paper. In eq. 8 we’re introduced to the Green St Venant strain tensor G
G = F^TF-I
The off-diagonal entries of G will be our shear values, and the diagonal our stretch values, and I is the identity matrix describing the rest length scale. With this info, we can write this in the form of two PBD constraints:
\begin {align*}
& C(p_0, p_1, p_2, p_3) = S_{ii} - s_i^2 \\
& C(p_0, p_1, p_2, p_3) = S_{ij}
\end {align*}Where S is FTF. The first constraint formulation is the stretch constraint where si2 is the rest stretch scale (previously the identity matrix), usually set to 1. This is showed in eq. 9 & 10 of the paper.
Multiplying a matrix’s transpose by the matrix itself (FTF) will always result in a symmetric matrix (Sij = Sji), hence we only need to care about the three diagonals and the three top right off-diagonal entries: i, j ∈ {1, 2, 3} and j ≥ i.
These are the basic constraint formulations, however, some more steps can be taken in order to increase the convergence rate. The first one is to linearize the stretch constraint which is explained in section 3.4 and eq. 14.
We can also modify the shear constraint since Sij will not only consider shear, but also length of the principle axes. This is because Sij is the dot product between vectors Fi and Fj, and the dot product will of course respect their lengths. We can therefore normalize this variable to ignore any stretching. This is explained in section 3.5 and eq. 15. Our updated constraints will look like this now, for stretch and shear respectively:
\begin {align*}
& C(p_0, p_1, p_2, p_3) = \sqrt{S_{ii}} - s_i \\
& C(p_0, p_1, p_2, p_3) = \frac{1}{|f_i||f_j|} S_{ij}
\end {align*}Particle Update
Next we need the gradients for each point w.r.t. Sij. Considering what FTF actually means, we can better understand how these gradients are calculated.
S_{ij} = f_i \cdot f_j = (Pc_i) \cdot (Pc_j)Here, c are the column vectors of Q-1 and f are column vectors of F. So we practically multiply the inverted “rest” matrix in to the deformed space P and compute the dot product between the axes. If Sij = 0, then the axes i and j are orthogonal, and the tetrahedra is not subjected to any shearing across that plane. If Sii = 0 then there’s no stretch along the i axis. The gradients of the particles for each value of Sij are then:
\nabla S_{ij} = [\nabla p_1, \nabla p_2, \nabla p_3]S_{ij} = f_jc_i^T + f_ic_j^T \\ci is transposed in order to define how the vectors should be multiplied, vector-vector multiplication isn’t generally defined, but in this case, we follow the definition of abT which results in a matrix where the columns are our point gradients for the constraint at Sij. Since we use p0 as our origin of the coordinate system, i.e. no strain will be measured from it, its gradient can be obtained as the negative sum of the other gradients:
\nabla p_0S_{ij} = -\sum_{k=1}^{3} \nabla p_kS_{ij}
Since we normalized the shear constraint value, we also need to take that in to consideration for the gradients of the shear constraints. This is explained in eq. 36. The updated gradient of the shear constraint become:
\nabla S^\prime_{ij} = \frac{1}{|f_i||f_j|}\nabla S_{ij} - \frac{|f_j|^2 f_i c_i^T + |f_i|^2 f_j c_j^T}{|f_i|^3|f_j|^3}S_{ij}Now we can find the projection vectors for each point of the tetrahedra w.r.t. Sij:
\Delta p_k = -\lambda w_k \nabla p_kS_{ij}The λ in the function is the scale amount for our stretch and shear projections respectively:
\begin{align*}
&\lambda = 2\frac{\sqrt{S_{ii}}-s_i} {\sum_k w_k |\nabla p_k S_{ii}|^2} \sqrt{S_{ii}}, \\
&\lambda = \frac{{S_{ij}}} {\sum_k w_k |\nabla p_k S_{ij}|^2}
\end{align*}Just as for the volume constraint, we divide the constraint value by the sum of all inverse particle masses of the tetrahedra. This is the original formulation of the constraint demonstrated in the paper. All the particle projections are explained in the appendix. However, since we want to fit in the XPBD framework, we need to reformulate this slightly.
Reformulate for XPBD
We keep track of the Lagrange multipliers λ for each stretch/shear value of the exes of the constraint:
\Delta \lambda _{ij} = (-C(p_0, p_1, p_2, p_3) - \alpha_{ij} \lambda_{ij}) / (w_{sum} + \alpha_{ij})Here, wsum is again the sum of the inverse masses, and 𝛼 is the stiffness parameter: 𝛼 ∈ ℝ≥0. The updated projection vector for each particle is then:
\Delta p_k = \Delta\lambda_{ij} w_k \nabla p_kS_{ij}Then at each iteration of the solver, we update the Lagrange multipliers with their respective deltas.
Implementation
Here’s the final implementation of the constraint inside the OpenCL node of Houdini. One important note if anyone is implementing this for Houdini 19: SideFx seems to have reversed the order for matrix vector multiplication in their matrix.cl header. If using Houdini 19 or lower, use mat3Tvecmul instead of mat3vecmul.
void solveStrainTet(const float timeinc,
const int ptidx,
const int prmidx,
global int* pts,
global fpreal* pos,
global fpreal* Qinv,
global float* rl,
global fpreal* stretch_stiffness,
global fpreal* shear_stiffness,
global float* L,
global float* M) {
int pt0 = pts[ptidx];
int pt1 = pts[ptidx + 1];
int pt2 = pts[ptidx + 2];
int pt3 = pts[ptidx + 3];
fpreal3 p0 = vload3(pt0, pos);
fpreal3 p1 = vload3(pt1, pos);
fpreal3 p2 = vload3(pt2, pos);
fpreal3 p3 = vload3(pt3, pos);
// inverse masses
fpreal w0 = safediv(1.0f, M[pt0]);
fpreal w1 = safediv(1.0f, M[pt1]);
fpreal w2 = safediv(1.0f, M[pt2]);
fpreal w3 = safediv(1.0f, M[pt3]);
// Load Qinv into Q
mat3 Q;
mat3load(prmidx, Qinv, Q);
// Compute P from material coordinates x, y, z
fpreal3 x = p1 - p0;
fpreal3 y = p2 - p0;
fpreal3 z = p3 - p0;
mat3 P = {x, y, z};
// In order to access Q by index, flatten the matrix
fpreal QT[9] = {Q[0].x, Q[1].x, Q[2].x,
Q[0].y, Q[1].y, Q[2].y,
Q[0].z, Q[1].z, Q[2].z};
// Unroll loop for some performance
__attribute__((opencl_unroll_hint(3)))
for (int i = 0; i < 3; i++) {
for (int j = i; j < 3; j++) {
// Deformation gradient F = PQ^-1
fpreal3 fi = mat3vecmul(P, Q[i]);
fpreal3 fj;
if (i == j) {
fj = fi;
} else {
fj = mat3vecmul(P, Q[j]);
}
// S = F^TF
fpreal Sij = dot(fi, fj);
fpreal wi, wj, s1, s3 = 0.0f;
// Initialize some variables for shear constraints
if (i != j) {
wi = length(fi);
wj = length(fj);
s1 = 1.0f / (wi * wj);
Sij *= s1;
}
// d stores the gradients of each point
fpreal3 d[4];
d[0] = (fpreal3) {0.0f, 0.0f, 0.0f};
for (int k = 0; k < 3; k++) {
// eq. 30, 31
fpreal3 grad = fj * QT[k*3+i] + fi * QT[k*3+j];
if (i != j) {
// eq. 35
fpreal3 gradn = wj*wj * fi * QT[k*3+i] + wi*wi * fj * QT[k*3+j];
grad = s1*grad-(gradn/(pow(wi, 3) * pow(wj, 3))) * Sij;
}
d[k+1] = grad;
d[0] -= grad;
}
fpreal wsum = w0 * dot(d[0], d[0]) +
w1 * dot(d[1], d[1]) +
w2 * dot(d[2], d[2]) +
w3 * dot(d[3], d[3]);
if (fabs(wsum) < EPS) {
continue;
}
// dL = update to lagrange multiplier, C is contraint term
fpreal dL, C, alpha = 0.0f;
if (i == j) {
// Diagonal - stretch. eq. 35
alpha = safediv(1.0f, stretch_stiffness[prmidx * 3 + i]);
if (alpha <= 0.0f)
continue;
fpreal s = sqrt(Sij);
C = 2.0f * s * (s - rl[prmidx * 3 + i]);
wsum *= 2.0f * s;
}
else {
// off-diagonal - shear. eq. 34
alpha = safediv(1.0f, shear_stiffness[prmidx * 3 + i+j-1]);
if (alpha == 0.0f)
continue;
C = Sij;
}
// XPBD term
// Previous lagrange value
int lidx = prmidx * 10 + i+j*3;
fpreal l = L[lidx];
// Update lagrange multiplier
//alpha /= timeinc * timeinc;
dL = (-C - alpha * l) / (wsum + alpha);
L[lidx] += dL;
// update points with their respeicive deltas
p0 += dL * w0 * d[0];
p1 += dL * w1 * d[1];
p2 += dL * w2 * d[2];
p3 += dL * w3 * d[3];
}
}
// Store points back to tmp for write back kernel
vstore3(p0, pt0, pos);
vstore3(p1, pt1, pos);
vstore3(p2, pt2, pos);
vstore3(p3, pt3, pos);
}
In the kernel function, we then call this function if the constraint is enabled
#ifdef TET_CONSTRAINTS
solveStrainTet(timeinc, pts_index[idx], idx, pts, P, Qinv, restlength, stretch_resist, shear_resist, L, M);
#endifWe can now modify all the parameters however we want from the wrangle node. Ideally, turn all of it in to an HDA and promote all the parameters! For setting groups to be deformed, just set the point mass to 0 for any point not in the group.
Result
Here’s a simple test case running 100 constraint iterations with stretch and shear stiffness = 10 and volume stiffness = 1e+10. This was run over a geometry with 953 tetrahedra and took about 0.2 s / frame.
Conclusion
I remember struggling when first implementing this, so it was fun to go back and try and re-learn what I had done. I even realized some misstakes I had made before which I now addressed. I’m sure there’s more though, but this implementation does what I expect it to do, at least!
I am getting some jittering with really high stiffness coefficients, which I find hard to adress while still maintaining the correct material properties but since this is quasi-static deformer, it’s unfortunately prone to jittering, and I decided to leave it here. However, turning this in to a dynamic solver, which solves over time, is really straight forward! Just run the constraint inside a SOP solver, or implement it with a Gas OpenCL node in DOPs, and update the point deltas for each frame. This solves a lot of the jittering we get when running it statically, It also converges much better! But that’s not the topic of this post.
The performance isn’t fantastic either, and I’m sure with a little bit of love, this could be improved significantly! I think for example that creating a dedicated kernel for multiplying some of the matrices could be achieved by splitting this up in to multiple OpenCL nodes inside a compiled block in Houdini. I might explore that at some point, but for now, I just have to live with a slow deformer.