WIP of Vertex Block Descent (VBD) in Houdini. It runs natively without plugins, as god intended.

I ported it to OpenCL based on all official references (TinyVBD, Gaia, NVIDIA Warp, AVBD-2D).

Now supports stuff from Augmented Vertex Block Descent (AVBD) too! Note this currently affects AVBD constraints only.

As well as the OpenCL version, there's an old VEX version to show how it works. Both are included in the HIP files.

Thanks to Anka He Chen and Chris Giles for making these open source with permissive licenses!

• Mass-spring constraints (for strings)
• Accelerated convergence method (section 3.8)
• Hessian direct inverse

• VBD integration (inertia, inertia and acceleration, adaptive)
• Neo-hookean constraints (for tetrahedrons)
• General damping

• Ground collisions
• Ground friction
• Ground damping

• Spring constraints (currently missing 6 DOF support)
• Joint constraints (currently missing 6 DOF support)
• Hard constraints (adaptive stiffness for AVBD constraint types)
• Dual constraint updates
• Breaking constraints
• Hessian LDLT decomposition
• SPD hessian approximation (section 3.5)

• Vellum integration (1st order, 2nd order)
• Self and external collisions (very buggy, to be replaced soon)
• Stopped, permanent and animated pin support
• Custom forces support
• Graph coloring

• Hessian approximation using eigenvalue projection

• Steal from NVIDIA Warp
  • StVK energy definition (for cloth)
• Steal from Offset Geometric Contact or AVBD
  • Self collisions
  • External collisions
• Steal from AVBD
  • Motor constraints
  • Sticking contacts
  • Post stabilization
  • 6 DOF support (packed prims)
  • Hard constraints for collisions
  • Update all existing constraints to AVBD
• Steal from Stable Cosserat Rods
  • Cosserat rod energy definition
• Steal from myself
  • SDF collisions
  • Add bounce for collisions
  • Stop everything exploding all the time
  • Convert from HDALC -> HDA
• Touch grass

You heard it here first! The new Otis solver in Houdini 21 uses Vertex Block Descent!

They made some nice improvements for stability, including several new hessian approximations. It also runs in chunks of 16 for some reason. They even added new collision handling based on how VBD did it!

Currently it's specialized for muscle and tissue, so it might be a while until it gets support for rigid bodies (packed prims). I'm planning to add support for this soon, so we'll see who wins!

VBD is very similar to Vellum. The main differences are constraints and collisions.

Vellum uses a technique called XPBD (Extended Position-Based Dynamics). XPBD uses constraints to simulate soft body behaviour. Constraints are solved in parallel worksets (colors) in OpenCL for better performance. Colors are groups of constraints that aren't directly connected.

Cloth is a good example of a soft body. It's easy to bend but hard to stretch. In XPBD this is simulated with distance constraints. Distance constraints try to preserve their rest length. When you stretch or squash a distance constraint, it pulls the points towards the middle until they reach their rest length again. Since shortening one constraint makes others longer, it's an iterative process. It propagates over several iterations until everything converges to the target length.

VBD constraints are similar, but they're defined in terms of energy instead. The goal is reducing overall variational energy by reducing local energy per point. VBD constraints run over each point rather than each prim, meaning less worksets (colors) overall. The Graph Color node allows worksets for points as well as prims, so it works both for VBD and XPBD.

The image above is for mass-spring energy, based on the rest length like in XPBD. Rather than directly using the rest length, the energy contributions get added to a force gradient and a hessian matrix. Once all the contributions are added, the position gets updated.

vector force = 0;
matrix3 hessian = 0;

// For each connected constraint (and extra forces)
force += constraintGradient;
hessian += constraintHessian;

// Reduce the variational energy of the system
v@P += force * invert(hessian);

Here's a comparison between Vellum and VBD:

The most important part of VBD is the energy definition, which depends on the constraint type. Here's a few:

• Mass-spring (in TinyVBD but removed from Gaia)
• StVK (in NVIDIA Warp)
• Neo-hookean (in Gaia)
• Spring (in AVBD, mass-spring but with a different formula)
• Joint (in AVBD)
• Motor (in AVBD)

I've implemented most of these as different constraint types, so you can connect them together like in Vellum.

AVBD is an extension to VBD mainly to improve stiffness. It adds hard constraints, an extension to regular constraints meaning the stiffness gets changed adaptively.

Stiffness is stored on the prims, so both the points (primal elements) and prims (dual elements) must be updated. As you might expect, looping over both points and prims is around 2x slower. Luckily I found you can merge the dual solve logic by tracking how many points are updated. Once the last point is updated, it's safe to dual solve the prims.

Adaptive stiffness currently only affects AVBD constraints. Eventually I'll rewrite the other VBD constraints to use it too.

AVBD also includes rigid bodies in a hacky way. Instead of solving each point of the rigid body, they solve the entire body as one point with 6 DOF (like a packed prim). This allows each point to have a rotation, included in the hessian for better results.

I think this is misleading, it gives a false impression of stiffness and goes against the design of VBD, but I'll eventually include packed prims and their rotations in the hessians. For now AVBD constraints solve translation but not rigid rotation.

AVBD also includes a SPD hessian approximation which greatly improves stability, used on all constraints by default. However it causes instability for neo-hookean constraints, so it's optional for them.

Personally I think dual updates could help, but not in the form implemented by AVBD. Vellum (XPBD) already uses dual updates for everything, with much greater success since it updates both sides of the constraint at a time (see stiffness limits). Even though stiffness gets affected by dual updates in AVBD, it only gets applied in primal updates. In other words, AVBD runs 2x slower than VBD without the benefits of dual updating you'd expect from XPBD.

Despite the hype, sadly not. It's near impossible for VBD to be faster than XPBD.

On paper VBD is faster because it runs over points, which have less graph colors than prims. Graph colors control the number of worksets, meaning how many points can be processed at the same time.

While it's true points have less graph colors, for all constraints in VBD each point loops over its connections to compute their energy contributions. This adds tons of extra operations, often 2x more than XPBD!

Even worse, these operations are expensive matrix calculations. The worst culprit is neo-hookean constraints. Tetrahedral meshes are densely connected, and each tet computes a 9x9 matrix!

Like with Vellum (XPBD), stiff objects are limited by the number of constraint iterations and substeps. The more constraint iterations and substeps, the more accurately stiff objects are resolved.

VBD also has accelerated convergence method meant to improve convergence for stiff constraints. It's named "Improve Convergence" in the Advanced tab and disabled by default, as it tends to explode with high values. AVBD also adds dual solving meant to improve stiffness. This is used for all AVBD constraints.

Personally I think Vellum has a huge advantage over VBD for stiffness, since it updates prims instead of points. For distance constraints this means both sides of the edge. You can easily match the target length just by spacing the points to the exact distance. After one iteration, you can guarantee the edge is the right length. This is much harder with VBD, since it only updates one point at a time.

I'm still working on collisions, they aren't been implemented properly yet (apart from ground collisions).

VBD solves collisions as soft constraints, meaning collisions get added onto the force and hessian like everything else.

In practice this means other forces can overpower collisions. For example, stiffer materials than the ground can penetrate it. This can be fixed by increasing the stiffness of the ground, or reducing the stiffness of everything else.

AVBD adds hard constraints which should prevent this from happening, but I haven't implemented this for collisions yet.

VBD involves updating the position based on force elements and a hessian matrix:

v@P += force * invert(hessian); // force and hessian depend on the energy definition, eg mass-spring or neo-hookean

invert(hessian) is very unstable, so everyone tries to bandaid it in various ways. The VBD paper uses the determinant of the matrix:

if (abs(determinant(hessian)) > 1e-7) { // if |det(H𝑖)| > 𝜖 for some small threshold 𝜖
  v@P += force * invert(hessian);
}

This helps, but it also explodes when the values gets too large (for example with very stiff constraints).

The AVBD paper uses an approximation to make the hessian positive semi-definite. This massively improves stability so it's used on all constraint types by default, except neo-hookean which is optional. Neo-hookean often fully degenerates when it's enabled.

The recent Augmented Vertex Block Descent (AVBD) paper adds many improvements to VBD.

I asked the authors about some differences I noticed. They responded with lots of useful information. Thanks guys!

Hi Chris, In the original VBD paper and in TinyVBD, they used an acceleration method to improve convergence (Section 3.8). I noticed in AVBD there's no mention of this method. Was it causing too much instability? Thanks!

Hi, yeah we ended up not using the acceleration from VBD as it was in general kind of unstable and difficult to tune, even with the original VBD method. It would be interesting to explore other acceleration methods as future work though.

-Chris

No, we haven't looked into acceleration for AVBD.

-Cem

Hi Chris, I was wondering what type energy you used for constraints? There were multiple used in the VBD paper, including mass-spring, StVK and neo-hookean. It looks like you used mass-spring energy. Is this correct, or did you use neo-hookean? Thanks!

Hello, so you are correct, in our demos we only used a simple spring energy for the deformable examples, as we weren't focused on rehashing what the original VBD paper showed. However, in AVBD, you can use any energy that works in VBD, such as the ones you mentioned. This is because AVBD is purely an extension of VBD. The only thing to keep in mind with those more complex energy types, is that you need to be careful about how you solve each block since their hessians can be indefinite.

In general, you can follow the same pattern that AVBD uses for constraint energies. That is, decompose the hessian into an SPD part and a non-SPD part, then use the diagonal lumped approximation proposed in the paper for the non-SPD part. Hope that helps!

-Chris

No. The AVBD tests we have are for contacts and joints. VBD already covers soft bodies. AVBD makes no changes to that.

-Cem

Hi Chris, I was wondering if the order of these 2 lines is correct?

body->prevVelocity = body->velocity; 
if (body->mass > 0) 
   body->velocity = (body->position - body->initial) / dt;

It seems like the previous velocity should be set after it gets recalculated, instead of before.

I saw the same code in TinyVBD, but I believe it is a mistake. The opposite code is present in Gaia.

The current code is correct (and probably in TinyVBD as well), since we use prevVelocity to compute an acceleration estimate during the adaptive warmstarting at the beginning of the step:

float3 accel = (body->velocity - body->prevVelocity) / dt;

If we switched the order as suggested, then this acceleration would always be zero, and the adaptive warmstart would not help.

Hi Chris, I've been looking through the code for AVBD 2D and want to ask about the design. The way it's implemented seems strange, with many differences to the original paper. For example I was expecting cubes to be implemented as 4 points connected by hard constraints, but instead they're represented as a unique class of object (rigid).

To me it reads more like a rigid body solver with ideas from VBD on top. For example the loop of the solver goes over rigid bodies. This differs from the AVBD paper, which loops over each point. The code for springs also includes rotation, which is more a rigid body concept.

I'm wondering how this applies to a cloth sim for example. Would you create a rigid body for each vertex of the cloth? In this case, what role does the rotation have?

AVBD and even VBD are not tied to the representation of the degrees of freedom (ie particles, rigid bodies, etc). You can use any representation that you wish. In AVBD we mainly showcased rigid bodies as these were the easiest to show examples of hard constraints. However, you could certainly create hard constraints between particles in AVBD. Also, collisions solved using AVBD are applicable to both particles and rigid bodies. Practically, when solving a particle vertex, you'd be solving a 3x3 system. With rigids, a 6x6 system.

For a cloth sim, usually you'd use particles, and each AVBD vertex maps to one vertex of the rendered mesh - no reason to use rigid dofs for this. In general, just use whichever dof representation makes sense for what you are simulating - AVBD / VBD support all of them!

Hi Chris, I'm currently working on rewriting the logic for 6 DOF bodies based on the AVBD paper. I was hoping you could provide more information about the math used to translate between hessian and quaternion representations, as I couldn't find much information about it in the paper.

In my implementation I used the same conventions as Vellum (XPBD) in Houdini. Using those conventions, @w is the angular velocity (vector) and @orient is the angle (vector4 / quaternion).

This means equation 21 looks like this:

orient = normalize(orient + timeinc * 0.5f * qmultiply((quat)(w, 0.0f), orient));

And equation 20 looks like this:

w = (2.0f * qmultiply(orient, qinvert(orientprevious)).xyz) / timeinc;

Both of these seem to work correctly, and give matching results to other rigid body solvers in Houdini. Where it gets foggy is translating between hessians and quaternions. It would be amazing to get more information about how this works in 3D.

In AVBD-2D, the bottom section (position->z) and the sides of the hessian contain angular components:

H = {
    dxx.row[0].x, dxx.row[0].y, dxr.x, // dxr is angular
    dxx.row[1].x, dxx.row[1].y, dxr.y, // dxr is angular
    dxr.x,        dxr.y,        drr    // dxr and drr are angular
};

With 6 DOF, I'm guessing this becomes a 6x6 matrix where the corners are a mix of linear and angular, then angular in the bottom:

H = { linear_hessian, linear_angular_hessian
      linear_angular_hessian, angular_hessian };

If this is correct, I'm wondering how you reduced the 3x3 hessians (linear_angular_hessian and angular_hessian) into a single quaternion (@orient in my case) ?

The 2D version solves both the linear and angular components using the same equation. In other words, it does this for everything:

position += force * invert_LDLT(hessian);

I'm not sure how this works when you use quaternions. You can no longer simply add, since the angle is a vector4 and the hessians are matrix3:

orient += angular_force *  invert_LDLT(angular_hessian) // Doesn't work, angular_hessian is 3x3 and orient is vector4

I'm guessing the way this was handled is by flattening the 3 hessians to a vector representation:

vector vector_representation; // Something involving angular_hessian
orient = normalize(orient + 0.5 * qmultiply((quat)(flat_vector3_representation, 0), orient));

Is this correct? If so, I'm very curious about the exact logic used here.

So in 3D, your intuition is spot on. You'll first build a 6x6 system like:

H = { linear_hessian, linear_angular_hessian
      linear_angular_hessian, angular_hessian };

Then you solve for the 6d update:

dx = force * invert_LDLT(hessian)

Then:

position += dx(0, 3)
orient = normalize(orient + 0.5 * qmultiply((quat)(dx(3, 6), 0), orient));

Actually in our code, we just override the += operator for quaternions to perform the above quaternion update.

Conceptually, when computing dx, for the angular part we are really getting a rotation vector, which is then used to update the actual orientation representation. The quaternion math there is just converting the rotation vector update to a quaternion update.

Note that you don't have to use quaternions for your angular representation (which would change the angular update rule). For example, if your orientation was also using rotation vectors as the representation, then you can directly add the orientation update using plain +=. (I wouldn't recommend rotation vectors in practice since they have various issues, just illustrating the point)

Another representation that can work is a 3x3 rotation matrix. This would also have a different update rule for += (Rodriguez formula).

We are currently working on a 3d version of the 2d demo which should provide a good reference implementation. But hopefully that helps you in the meantime!

Forgot to mention, checkout this paper's appendix (Geometric Stiffness for Real-time Constrained Multibody Dynamics by Andrews et al.) for some derivations of the hessian for various constraint types.

You'll see references to lambda, you just replace those with the force F = K * C + lambda.

VBD's high level design is simple, it's really just 3 steps. AVBD adds another 2 steps.

Add the velocity to the position (same as Vellum). VBD uses a warmstarting strategy to scale the gravity term below.

v@pprevious = v@P;

// First-order integration, same as Vellum
v@v += v@gravity * f@TimeInc;
v@inertia = v@P + v@v * f@TimeInc;
v@P = v@inertia;

AVBD adjusts the stiffness based on lambda and penalty. They get dampened by alpha and gamma before constraint solving to prevent explosions.

// Warmstart the dual variables and penalty parameters (Eq. 19)
@lambda *= alpha * gamma;

// Penalty is safely clamped to a minimum and maximum value
@penalty = clamp(@penalty * gamma, PENALTY_MIN, PENALTY_MAX);

// If it's not a hard constraint, we don't let the penalty exceed the material stiffness
@penalty = min(@penalty, f@stiffness);

This is the hard part. The core of VBD is moving the position based on a force gradient and a hessian matrix.

When the hessian doesn't explode, moving the position reduces the overall variational energy.

In AVBD, lambda and penalty change the results of accumulateMaterialForceAndHessian().

vector force = 0;
matrix3 hessian = 0;

// Add contributions to force elements and hessian
accumulateInertiaForceAndHessian(force, hessian); // Contributions due to mass and inertia
accumulateMaterialForceAndHessian(force, hessian); // Contributions due to connected constraints (eg mass-spring)
accumulateDampingForceAndHessian(force, hessian); // Contributions due to damping
accumulateBoundaryForceAndHessian(force, hessian); // Contributions due to boundaries (eg ground plane)
accumulateCollisionForceAndHessian(force, hessian); // Contributions due to collisions

v@P += force * invert(hessian); // Reduce the variational energy of the system

Clamp lambda and penalty to prevent them exploding again. The C variable depends on the constraint type.

// Use lambda as 0 if it's not a hard constraint
float lambdaTmp = isinf(stiffness) ? lambda[i] : 0;

// Update lambda (Eq 11)
lambdaTmp = lambda[i] = clamp(penalty[i] * C[i] + lambdaTmp, fmin[i], fmax[i]);

// Update the penalty parameter and clamp to material stiffness if we are within the force bounds (Eq. 16)
if (lambdaTmp > fmin[i] && lambdaTmp < fmax[i]) {
    penalty[i] = min(penalty[i] + beta * abs(C[i]), min(PENALTY_MAX, stiffness));
}

Update the velocities based on the change in position (same as Vellum).

// First-order velocities
v@v = (v@P - v@pprevious) / f@TimeInc;