Finite Element Method (FEM) Simulation

CS 482 Lecture, Dr. Lawlor

Our goal is to be able to simulate solid objects, like wood, metal, and bone, as they deform, bend and break.  This is actually a problem from materials science, which uses a thick pile of specialized terminology:
Stress and strain can sometimes be related via a stress-strain curve, with the strain on the X axis, and stress on the Y axis.  This is a somewhat misleading curve, because stress also depends on the material's internal state, so the stress-strain curve can become multi-valued.  In general, stress-strain behavior is usually modeled via a Constitutive Equation, such as:
For liquids, the terminology is a bit different.  Most liquids respond purely elastically to stresses that are equal in all directions (called hydrostatic stress), but flow in response to shear.  A viscous material, like honey, pushes back strongly in response to shear stress.  As many solids approach their melting point, they often display a form of viscous flow called creep; at normal temperatures ice creeps rapidly, causing glaciers to flow downhill.  Non-Newtonian materials respond differently depending on the rate at which the deformation is applied.   For example, ketchup is thixotropic, getting less viscous at higher strain rates; this is why ketchup squirts easily, but takes a long time to flow down the bottle.  By contrast, cornstarch is shear thickening--it gets stronger at higher strain rates, so in a large vat of cornstarch, one person can float in it while another person walks over it (quickly!).

In real life, material stresses often depend somewhat on all the above factors--real materials are elastoviscoplastic and subject to creep and fracture.  Luckily, most of the time only one or two of these phenomena are important at once!

Practice--Finite Elements

So that's what really happens--each point in the material examines its local strain and internal properties, and then figures out how much stress to apply to its neighboring points.  To discretize this, we usually just take an element of finite, nonzero size, like a triangle or tetrahedron, calculate the strain over that element, compute the stress based on that strain, and apply that stress to element's boundary points ("nodes").  Like any discrete computation, this isn't exact, but it's close enough to use to design working machines.

In code, a typical finite element simulation looks a lot like our little spring system:
	// zero out node forces
for (n=0;n<nodes.size();n++) nodes[n].F=vec3(0);

// add element forces to nodes
for (e=0;e<elements.size();e++) {
mat3 strain = compute_strain(e); /* looks at node positions */
mat3 stress = constitutive_model(strain,e); /* e.g., elastoplastic */
for (n=0;n<elements[e].numnodes;n++) /* apply our stress to each node */

// move nodes based on net force
for (n=0;n<nodes.size();n++)
vec3 A=nodes[n].F/nodes[n].m; /* F=mA */
nodes[n].V+=dt*A; /* dv/dt = A */
nodes[n].P+=dt*V+0.5*dt*dt*A; /* dp/dt = V */
All the good stuff above is hidden inside the constitutive_model call.


I've written several FEM implementations in 2D and 3D, and they work, but they tend to freak out in one situation or another, such as large rotations, large deformations, etc.  This is a common bug.  For engineering type simulations, you don't care what happens if your bridge gets tied into a knot (it's not a bridge by that point anyway); but for games, or simulations where large deformations occur, you really want something reasonable to happen. 

I'm excited by Irving, Teran, and Fedkiw's 2004 paper, which proposes a simple but robust matrix-based approach for keeping track of the rest configuration and computing strains, from which you can then calculate stresses according to any constitutive response you like.  He's got a bunch of amazing simulation videos on his site.  The one downside is the dense thicket of matrix math required to disentangle object deformation from simple rotation.  This matrix math is both difficult to write and slow to run: in 2004, around 1 second per frame (!) for just 11K elements.  A paper in 2005 by Nesme, Payan, and Faure recommended a QR decomposition to find the rotation portion of the matrix, and so achieve 30fps on small meshes of a few thousand elements.

See this matrix notation refresher for DirectX-style row major matrices, as used by BABYLON.Matrix.