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. 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).
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.
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). Of the two versions the authors suggested, I implemented the most basic one, which doesn't handle very low Poisson's ratios but looks nice otherwise. It has a tendency to break down catastrophically under tension, especially with low step counts or low compliance. Linear triangle and tetrahedron elements (T3/T4) exhibit severe volumetric locking with the "Uniform" Tri/Tet Pattern. Linear quad elements (Q4) eliminate most of the locking, but still act unnatural at higher Poisson's ratios. Linear 3D hexahedron elements (H8) do quite a bit better, and quadratic elements (T6/T10/Q9/H27) can handle very high Poisson's ratios before you start seeing locking effects.
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 ๐ผโ - 3 (i.e. tr(๐นแต๐น) - 3) term. On the positive side, this eliminates the locking in Q4/H8 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. Internally, this option always attempts to maintain full incompressibility, though I added an offset in the Neo-Hookean constraint so that the motion will match the Pixar method for Poisson's ratios close to 0.5.
Serial: Serial mixed method. Only here for demonstration purposes, as it's worse than the previous option. The previous option solved the Neo-Hookean and volumetric constraints simultaneously per element using the matrix form of XPBD. This option uses two regular XPBD steps in sequence, one for the Neo-Hookean and one for the volumetric. For a given compliance, requires many more steps per second to converge, so it's significantly less efficient.
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. Generally works best when Vol. Compliance is set to 0.
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 T3 element type, only integrates the hyperelastic invariants, and then calculates the energy using the summed invariants. This greatly reduces locking (this same idea would apply to the T4 element, but isn't implemented). Per constraint, this method is significantly more expensive than constraining per-element, though it generally requires fewer Steps/Second to converge. For T3/T4 elements, per vertex constraints are probably a better way to go than per-element.
V Mix: A version of the Mixed Subregion constrain implemented per-vertex. This is pathologically slow for higher order elements, but for T3/T4, it might be the best choice.
V Skin: Yeoh Skin implemented per vertex. Like V Pix, should be implemented with pre-integrated invariants for T3/T4, but currently isn't. At least that makes it a good point of comparison to see the advantage of the pre-integrated invariant technique.
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 T3/T4 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.
Vol. Compliance: For mixed energy functions, controls the relative compliance of the volume preservation constraint. Generally, you want this as low as possible. Only provided as a slider because once it goes below the default 0.005, can cause strange oscillations at certain combinations of Steps/Second and Compliance.
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