# 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:
• Strain measures how much the material has been deformed, usually measured from some "rest" configuration.  The symbol is usually e or epsilon.  Units are usually in inches of displacement per inch of material.  For Hooke's law on a one-inch spring, strain is the end displacement x.  Strain can be simulated (for example, integrated from particle forces), or measured with either accurate calipers or a strain gauge.
• Stress is the force the material exerts to fight back against strain.  The symbol is usually the Greek letter sigma σ. The units for stress are in pounds of force per square inch of material.  In 1D material, stress is a scalar. In 3D material, stress can be different in different directions, but it can always be represented by a 3x3 matrix, the "Cauchy stress tensor", which maps each direction to the force the material exerts along that direction.
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 a purely elastic material, the stress-strain curve is linear, like a spring: push on it, and it deforms (or deform it, and it pushes); let it go, and it returns to its original shape and stops pushing.  This linearity is captured in Hooke's Law: F=-kx in physics terminology; or σ=Ee (stress is stiffness time strain) for engineers.   The E there isYoung's modulus, which measures the stiffness of an elastic material; this varies from a few thousand psi for rubber to a few million psi for metals.  In general, strain is a combination of uniaxial and shear.  Of course, under enough strain all real materials break down and go nonlinear at some point--that is, no material completely follows Hooke's Law (it's more a guideline, really.)
• By contrast, in a plastic material, like wet clay, the stress-strain curve is highly nonlinear--enough strain causes the material to yield and permanently deform.  Many materials are work hardening, displaying an increase in stiffness upon plastic deformation.
• A brittle material, like dried clay, fractures into pieces instead of deforming.  The stress-strain curve breaks down entirely.  Most metals are elastic for small strains, plastic for some distance outside that, and then brittle.
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 */			nodes[n].F+=apply(stress,elements[e].nodes[n]);	}		// 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.

## Implementation

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.