Show explainer text

This web page demonstrates the finite element method (FEM) implemented with small step XPBD (paper). The jumping off point for this implementation was Position-Based Simulation of Continuous Materials (paper), which suggests directly constraining non-linear hyperelastic potential energy functions. I've extended it to support more advanced element types than just triangles and tetrahedra and I've experimented with mixed method constraints, as well. Recently, A Constraint-based Formulation of Stable Neo-Hookean Materials (paper) was published, which shows how to more precisely implement energy constraints with XPBD, and I have incorporated those suggestions. You can click/drag on the elements in the viewport above to move them around, as well as tweak many aspects of the simulation with controls down below (people on smaller screens may want to zoom out a little).

Notes on the controls:

You can enable a second instance of the simulation by selecting a non-null element type in the right hand column, which should make it easy to compare and contrast different settings. The column of buttons between the two sets of controls are link toggles. When the link is turned on, the corresponding right-hand control will dim and take on the value of the left hand control. Press 'L' to toggle most of the links at once. If the simulation ever becomes catastrophically unstable, 'R' will reset it.

Control descriptions:

Element Type: The type of the individual finite elements. All shapes other than T3 and T4 are non-linear, which is an issue for many other integration methods, but is no sweat for XPBD. T6, T10, Q9 and H27 are higher order quadratic elements (and therefore especially non-linear).

Unbalanced configurations of the tetrahedral elements (T4 and T10) experience significant non-physical ghost forces. This is especially obvious if you set the Shape to "Line", the Tri/Tet Pattern to "Uniform" and set all forms of drag and damping to 0. Switching the Tri/Tet Pattern to "Mirrored" and then choosing a balanced Shape like "Beam M" should eliminate most of the chaotic movement.

Energy Function: Controls the constraint used by the finite elements. These are worth describing in detail.

Pixar: A simple and effective energy function outlined in Pixar's Stable Neo-Hookean Flesh Simulation (paper). It has a tendency to break down catastrophically under tension, especially with low step counts, low compliance or high Poisson's ratios. First order elements (T3/Q4/T4/H8) all exhibit minor "locking" where movement is artificially restricted for numerical reasons. At nearly incompressible Poisson's ratios (>0.49), the locking becomes quite severe. Quadratic elements (T6/T10/Q9/H27) do better, but even they lock somewhat at very high Poisson's ratios.

Reduced: Fully reduced integration. By performing numerical integration with fewer samples than necessary, we can see locking is gone! Unfortunately, so is any semblance of structural stability. This is probably not a practical mode, but I thought it was interesting for demonstration purposes. (This is only enabled for the quad and hex elements.)

Sel: Selectively reduced integration. Only uses reduced integration on the volume-preserving portion of the energy function. Fully integrates the Neo-Hookean ๐ผโ‚ (i.e. tr(๐นแต€๐น)) term. On the positive side, this eliminates the locking in quad/hex elements and the algorithm is significantly faster, as ๐ผโ‚ can be factored into O(๐‘›ยฒ) pre-integrated terms (where ๐‘› is the number of nodes in an element), and performing reduced integration on the volume-preserving component requires a fraction of the samples compared to full integration. On the negative side, some amount of "hourglassing" is introduced, manifesting as zig-zag artifacts.

Mixed: Mixed method. Instead of trying to solve each element with a single constraint function, breaks the solution into separate incompressible Neo-Hookean and volumetric constraints. Compared to the previous irreducible energy functions, incompressible Neo-Hookean with simple J = 1 volume constraints are quite stable, and converge to full stiffness with relatively low step counts. This is especially true at very high Poisson's ratios, where the Pixar energy functions required many Steps/Second to converge.

Serial: Serial mixed method. While all other mixed methods solve deviatoric and volumetric constraints simultaneously with the matrix form of XPBD, this option uses two regular XPBD steps in sequence. This option is somewhat historical, as in previous versions of this demo simultaneous solve was a huge advantage. With the updates from A Constraint-based Formulation of Stable Neo-Hookean Materials the differences are now more subtle.

Sub: Optimized mixed method with separate volume constraints for each subregion of the higher order shape functions. T6 uses 3 volume constraints, T10 and Q9 use 4 and H27 element uses 8. The Neo-Hookean is calculated with the prefactored integral optimization and each volume subregion is calculated with a single quadrature point. A very good energy function for the Q4 and H8 shapes, but the O(n^3) cost of solving 9 constraints simultaneously significantly hurts the performance for H27s.

YeohRubber: Replaces the Neo-Hookean constraint of the mixed method with a Yeoh constraint tuned approximately for rubber. Less stable than Neo-Hookean and requires high step counts to resolve properly. Fairly similar response to Neo-Hookean, but it resists extreme extension a bit more.

YeohSkin: Yeoh constraint tuned approximately for skin. When given enough steps to properly resolve, displays a markedly different response than Neo-Hookean. It's very loose under low amounts of extension, and then quickly becomes very stiff as it gets stretched.

V Pix: A version of the Pixar energy function, but instead of constraining per-element, constrains per-vertex by summing the energy integrals of all surrounding elements. For the linear tri and tet elements (T3/T4), only integrates the hyperelastic invariants, and then calculates the energy using the summed invariants. This greatly reduces locking, though does not eliminate it completely. Per constraint, this method is significantly more expensive than constraining per-element, though it generally requires fewer Steps/Second to converge.

V Mix: A version of the Mixed Subregion constraint implemented per-vertex. This is too slow for higher order elements, but for T3/T4, it might be the best choice.

V Skin: Yeoh Skin implemented per vertex.

Shape: The shape of the aggregate FEM block. I've done most of my testing using Beam L, Box L and Armadillo. Armadillo loads the Stanford Armadillo mesh or something like it. I was unable to get an automatic mesher to work for hex meshes, so for H8/H27 you get my hand made Armadillo-Bear.

Tri/Tet Pattern: The pattern of the orientations of tri and tet elements. "Uniform" most closely resembles the distribution you'd get in a generated mesh, and also unfortunately has serious locking problems for first order elements. "Mirrored" is mainly thrown in out of interest, as it significantly reduces locking, though it has subtle artifacts of its own.

Steps/Second: Affects the fidelity of the simulation. Has linear CPU cost. Depending on the combination of energy function, compliance and how much stress the simulation is under, may need to be increased for the simulation to converge. When trying out new settings, it's a good idea to crank this value to as high as your machine can handle, and then lower it to just above the point where the sim breaks down to get best performance.

Volume Passes: Controls how many additional volume-only constraint solves are run after the regular energy function constraints. Only useful when going for incompressibility, but in certain cases might allow significantly turning down the Steps/Second, saving CPU time. Especially useful for Q4 and H8 elements, as their volume constraints are relatively cheap. Currently not implemented for vertex constraints.

Gravity: Gravitational constant in meters per second. The sim is quite "zoomed in," which is why the default gravity is so low.

Compliance: Inverse of stiffness. Thanks to XPBD, it's independent of timestep, and I normalized it against element size, as well. Should be within a constant multiple of physical compliance, but I haven't yet sought out a source of ground truth to calibrated it.

Rayleigh Damp: Rayleigh Damping as described in the original XPBD paper. A bit tricky to tune, as it depends on compliance and has a narrow range where it looks good. Also, while the damping calculations themselves are very cheap, even moderately high damping requires very high Steps/Second, which can be a performance issue.

PBD Damp: The non-rigid body mode damping described in the venerable Position Based Dynamics paper. Similar to Rayleigh damping, though a bit more aggressive. Useful because it does not depend on compliance or step count, though, on the other hand, it's cost is not de minimis.

๐‘ฃ *= (1-x) Drag: Damps all motion non-physically. Basically what it says on the tin, though corrected to be independent of timestep.

Poisson's Ratio: As discussed in the energy function section, controls how volume preserving the sim is. All of the energy functions are geared towards high Poisson's ratios, and the non-Pixar ones want it very close to 0.5, indeed. (Though setting Poisson's ratio to exactly 0.5 will cause a divide by zero, and all energy functions will explode.)

Orientation: Which way you want the block to hang.

Lock: Left and Right lock the corresponding side of the block in place. Rotate will rotate the right side of the block around the block's midpoint, which can be useful when examining volumetric locking.

Left-Right Separation: When both the Left and Right locks are selected, allows you to test stretching and compression.

Wonkiness: Perturbs the verts to demonstrate support for non rectilinear elements. Values above 0.5 are out of spec, but you can try setting it higher to see how the sim fails.

Time Scale: Allows you to slow down the passage of time in the simulation. Useful to observe how things explode when your settings aren't quite right.

Hide explainer text