Houdini VBD

Houdini Vertex Block Descent

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

There\’s an OpenCL version for performance, and an old VEX version to show how it works. Both are included in the HIP files.

I ported everything from TinyVBD, some bits from AVBD, Gaia and some ideas from the papers.

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

Todo

• Steal from TinyVBD
  • Mass-spring energy definition (for strings)
  • Accelerated convergence method (section 3.8)
• Steal from NVIDIA Warp
  • StVK energy definition (for cloth)
• Steal from Gaia
  • Neo-hookean energy definition (for tetrahedrons)
  • Damping
  • Boundary collisions
  • Self collisions
  • External collisions
  • Friction
• Steal from AVBD
  • LDLT decomposition to improve speed
  • SPD hessian approximation (section 3.5) to improve stability
  • Hard constraints
• Touch grass

What\’s Vertex Block Descent?

VBD is very similar to Vellum. You might call it Vellum 2.

Vellum uses a technique called XPBD (Extended Position-Based Dynamics). XPBD uses constraints to simulate soft body behaviour. Constraints are solved in parallel workgroups (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 bendy but stiff in terms of edge lengths. In XPBD this is simulated with distance constraints. Distance constraints try to preserve their rest length based on stiffness. 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 and vice versa, 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. They also run over each point rather than each constraint, meaning less workgroups (colors) overall. The Graph Color node allows workgroups for points as well as prims, so it works both for VBD and XPBD.

Here\’s a quick comparison between VBD and XPBD:

The most important part of VBD is the energy definition, but no one seems to agree on this.

I\’ve seen many different energy definitions, including mass-spring (used by TinyVBD and AVBD, but removed from Gaia), StVK (used by NVIDIA Warp) and neo-hookean (used by Gaia).

Currently mass-spring and neo-hookean are supported, StVK is coming soon.

How does Vertex Block Descent run?

Ignoring collisions, VBD is really just 3 steps. These steps are nearly identical to XPBD apart from the constraints.

1. Integrate the positions

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

// First-order integration
v@pprevious = v@P;
v@P += v@v * f@TimeInc * v@gravity * f@TimeInc * f@TimeInc;

2. Apply the constraints

The core idea of VBD is updating the position based on force elements and a hessian.

If these are correct, moving the position should reduce the overall variational energy.

Caution

This should be run in workgroups based on graph coloring!

If points move while their neighbours access them (like if running in sequential order), it breaks the assumption used by VBD:

We adjust each vertex separately, assuming the others remain fixed

This causes growing error each iteration, leading VBD to explode much more than usual.

vector force = 0;
matrix3 hessian = 0;

// Add influences to force elements and hessian
accumulateInertiaForceAndHessian(force, hessian); // Influences due to mass and inertia
accumulateMaterialForceAndHessian(force, hessian); // Influences due to constraints (eg mass-spring or neo-hookean)
accumulateDampingForceAndHessian(force, hessian); // Influences due to damping
accumulateBoundaryForceAndHessian(force, hessian); // Influences due to boundaries (eg floor planes)
accumulateCollisionForceAndHessian(force, hessian); // Influences due to collisions

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

3. Update the velocities

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

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

Why does stiffness have a limit?

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 in the \”Advanced\” tab and disabled by default, as it tends to explode with high values, but worth a try.

AVBD adds hard constraints which should resolve much faster, but I haven\’t implemented this yet.

Why do collisions not work sometimes?

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

Why does it explode randomly?

This solver used to explode every 5 seconds, but now it\’s much better. Explosions are a common issue with VBD.

Like mentioned, VBD involves updating the position based on force elements and a hessian:

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 new AVBD paper uses an approximation to make the hessian positive semi-definite. This massively improves the stability.

AVBD Q&A

There\’s a new paper called Augmented Vertex Block Descent (AVBD). It adds many improvements to VBD.

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

Missing accelerated convergence

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

Energy definition used

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

Previous velocity definition

Note

I noticed this while looking into Vellum. Vellum uses 4 variables to track the previous 2 values of position and velocity:

• @pprevious (@P 1 substep ago)
• @plast (@P 2 substeps ago)
• @vprevious (@v 1 substep ago)
• @vlast (@v 2 substeps ago)

Vellum sets all of these at the start of each substep. They\’re needed for 1st and 2nd order integration.

However, TinyVBD and AVBD set @vprevious in a different place. I thought this was a typo, but turns out it\’s not.

Hello, 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.