Hi guys,

I've been working on a small step XPBD-based implementation of FEM off and on for the last several months, and I think I finally have something that someone might find interesting. Here's a WASM demo showing the current state of things: https://jak-xyz.github.io/xpbd-fem/. At some point I had delusions of grandeur of building an interactive walkthrough to show off the features, but unfortunately right now you're just presented with a wall of controls and an even bigger wall of text explaining them.

I think the two most interesting things about the demo are that it attempts to tackle incompressible FEM and it shows off higher-order elements (like quadratic quads/hexahedra). Until you peek at the demo this next sentence won't make too much sense, but the performance of the Sub energy function with Q4 and H8 element types seems promising. You can keep the number of steps/second pretty low and still maintain good stability and volume preservation, especially with 2-3 extra volume constraint passes.

Interestingly, I've found the WASM version running in Chrome to perform indistinguishably from native version. I have not attempted to add WASM SIMD support, but with an older version of the code I was getting a full 8x speedup with a trivial structure of arrays AVX2 implementation in a native build. The one issue I had to work around is that small step XPBD tends to require double precision floats and, in fact, single precision floats were insufficient for storing nodal position values. The constraints themselves, on the other hand, are calculated in normalized element space, which is not so sensitive, so the vast majority of calculations can be on regular floats.

The code is here: https://github.com/jak-xyz/xpbd-fem, though I abandoned all attempts at tutorializing anything, so it's probably only readable to someone with a pretty good handle on FEM to begin with. That said, someone implementing their own XPBD FEM may find it useful to look at the code for normalizing things like compliance and damping against element size.

In the not-too-distant future I intend to add support for higher-order triangle and tetrahedron element types. I may also investigate storing additional values than just positions per node. I briefly tried storing volume per node, and that fixed hourglassing in the fully integrated mixed energy function, but I didn't really dive into its performance characteristics, so I might do that later.

Finally, I'd like to give a shout out to Matthias Müller and his SIGGRAPH(?) 2008 course notes on linear, co-rotational FEM (https://cgl.ethz.ch/Downloads/Publicati ... enotes.pdf). I've attempted to learn about FEM several times in the past, but always bounced off. Coming from a game development background, this is still the clearest intro to FEM I've seen, and it was enough to let me get a toe-hold into the wider subject.

## Small Step XPBD FEM Demo

### Re: Small Step XPBD FEM Demo

Found a sudden burst of inspiration to post an update to the demo ;p. The 3 major additions are a few new element types (second-order triangles and first- and second-order tetrahedrons), Stanford armadillo meshes for each of the element types, and a new set of energy functions that constrain per-vertex instead of per-element. The demo is still at https://jak-xyz.github.io/xpbd-fem/. If you've loaded the page before, you may need to press the "Reset Settings" button to clear your settings cache.

In terms of the new element types, there's not much to say, in isolation.

Adding the armadillo meshes was a learning experience, however. Now, it started gently enough with the tri and tet meshes. For the tri armadillo, I followed an old blog post by Miles Macklin. For tets, I was able to use TetGen without much fuss. The quad mesh was a bit more of a journey, but I eventually got gmsh to spit one out. The boundary could still use some refinement, but it seems like getting good-quality automatically-generated 2D quad meshes is tractable. The hex "armadillo," on the other hand, is where things really took a turn. Hex meshing is a much more difficult problem than tet meshing, and given that, the progress people have made on it is impressive. I tried a variety of automatic meshers, and I got farthest with the PolyCut toolchain from the University of British Columbia here: http://www.cs.ubc.ca/labs/imager/tr/2018/HexDemo/. It was able to generate some nice meshes, but not at a low-enough resolution for a browser demo. In the end, I ended up just hand-modeling the armadillo-adjacent creature you see in the demo.

Now, up to this point, the H8 linear hex element seemed like by far the most promising direction, since it's locking-resistant and it's the most performant element on a per-vertex basis by a large margin. And I still think if you're in a situation where a hex mesh is available, it's the thing to use. For instance, adding a hex-based flesh layer to a skinned character seems possible. For arbitrary volume meshes, though, hex mesh generation is too much of an impediment, so I decided to look a bit more into getting tets up to snuff.

Like with triangles, if you use the basic technique of constraining the energy function per-element on linear tets, you get pretty bad volumetric locking. Quadratic tets don't lock too much, but they are vastly more costly. Eventually, inspired by this paper: https://graphics.pixar.com/library/Inco ... /paper.pdf, I tried moving the constraint from the element to the one-ring around each vert. Each constraint is now quite a bit more costly, but you only have about 1/6th as many (in the case of tets, at least). Also, with the Gauss-Seidel solver, you need fewer Steps/Second to converge. With my current implementation, performance is a wash for tets between using element-centric constraints and one-ring-centric constraints, with the big tie-breaker being that one-ring constraints don't lock. All other element types are at least a bit slower. For the very simple mixed Neo-Hookean energy function (represented as "V Mix" in the demo), that's all there is to it: integrate each of the energies of the surrounding elements, divide by 4 (to represent the portion of the tet that belongs to the given vert), and you're set. For more complex energy functions, like Pixar or Yeoh Skin, you have to be a bit more careful. Integrating the energy function per-element still leads to locking. Instead it seems like you can just integrate the hyperelastic invariants I1 and J, and then use the weighted sum of those to calculate the energy function. I hacked in this technique for the T3 triangle version of "V Pix," though T3 "V Skin" and both T4 "V Pix" and "V Skin" still integrate the full energy function per-element, and thus lock.

Looking towards the future, I'm thinking of ditching the higher-order elements. It was educational to implement them, but I haven't been able to get them performance-competitive with the first-order elements on a per-vert basis. Even first-order tets are about 3x slower to achieve convergence compared to first-order hexes, but I think if I gave them a specialized pathway I could at least make up some of that difference. Along those lines, I might investigate WASM SIMD, since browsers are starting to roll out support for that, and it would be a good way to keep me honest about the complications of parallelizing one-ring constraints vs. simple element-constraints.

In terms of the new element types, there's not much to say, in isolation.

Adding the armadillo meshes was a learning experience, however. Now, it started gently enough with the tri and tet meshes. For the tri armadillo, I followed an old blog post by Miles Macklin. For tets, I was able to use TetGen without much fuss. The quad mesh was a bit more of a journey, but I eventually got gmsh to spit one out. The boundary could still use some refinement, but it seems like getting good-quality automatically-generated 2D quad meshes is tractable. The hex "armadillo," on the other hand, is where things really took a turn. Hex meshing is a much more difficult problem than tet meshing, and given that, the progress people have made on it is impressive. I tried a variety of automatic meshers, and I got farthest with the PolyCut toolchain from the University of British Columbia here: http://www.cs.ubc.ca/labs/imager/tr/2018/HexDemo/. It was able to generate some nice meshes, but not at a low-enough resolution for a browser demo. In the end, I ended up just hand-modeling the armadillo-adjacent creature you see in the demo.

Now, up to this point, the H8 linear hex element seemed like by far the most promising direction, since it's locking-resistant and it's the most performant element on a per-vertex basis by a large margin. And I still think if you're in a situation where a hex mesh is available, it's the thing to use. For instance, adding a hex-based flesh layer to a skinned character seems possible. For arbitrary volume meshes, though, hex mesh generation is too much of an impediment, so I decided to look a bit more into getting tets up to snuff.

Like with triangles, if you use the basic technique of constraining the energy function per-element on linear tets, you get pretty bad volumetric locking. Quadratic tets don't lock too much, but they are vastly more costly. Eventually, inspired by this paper: https://graphics.pixar.com/library/Inco ... /paper.pdf, I tried moving the constraint from the element to the one-ring around each vert. Each constraint is now quite a bit more costly, but you only have about 1/6th as many (in the case of tets, at least). Also, with the Gauss-Seidel solver, you need fewer Steps/Second to converge. With my current implementation, performance is a wash for tets between using element-centric constraints and one-ring-centric constraints, with the big tie-breaker being that one-ring constraints don't lock. All other element types are at least a bit slower. For the very simple mixed Neo-Hookean energy function (represented as "V Mix" in the demo), that's all there is to it: integrate each of the energies of the surrounding elements, divide by 4 (to represent the portion of the tet that belongs to the given vert), and you're set. For more complex energy functions, like Pixar or Yeoh Skin, you have to be a bit more careful. Integrating the energy function per-element still leads to locking. Instead it seems like you can just integrate the hyperelastic invariants I1 and J, and then use the weighted sum of those to calculate the energy function. I hacked in this technique for the T3 triangle version of "V Pix," though T3 "V Skin" and both T4 "V Pix" and "V Skin" still integrate the full energy function per-element, and thus lock.

Looking towards the future, I'm thinking of ditching the higher-order elements. It was educational to implement them, but I haven't been able to get them performance-competitive with the first-order elements on a per-vert basis. Even first-order tets are about 3x slower to achieve convergence compared to first-order hexes, but I think if I gave them a specialized pathway I could at least make up some of that difference. Along those lines, I might investigate WASM SIMD, since browsers are starting to roll out support for that, and it would be a good way to keep me honest about the complications of parallelizing one-ring constraints vs. simple element-constraints.

### Re: Small Step XPBD FEM Demo

The excellent paper Miles Macklin posted yesterday (A Constraint-based Formulation of Stable Neo-Hookean Materials) elucidated some things about the relationship between constraints and energy in XPBD that I had been fuzzy on, and I couldn't resist jamming the corrected mixed Neo-Hookean energy function into the demo. The updated version is still at https://jak-xyz.github.io/xpbd-fem/, and I moved the previous version to https://jak-xyz.github.io/xpbd-fem/v2. The biggest issue with the previous version is that I was constraining the full Neo-Hookean energy, instead of the square root. The irreducible Pixar Neo-Hookean energy function was similarly squared compared to what it should be, and I went ahead and fixed that as well. With the fixed energy function, the material is more permissive of large extensions and Poisson's ratio is respected in most modes. The "Mixed" energy function (which is a fully integrated mixed energy function) also now works correctly with quadratic elements.