The third assignment will introduce you to the most common numerical method for simulating, well, almost anything in Computer Graphics, the mighty Finite Element Method.
On all the platforms, we will assume you have installed cmake and a modern c++ compiler on Mac OS X¹, Linux², or Windows³.
We also assume that you have cloned this repository using the --recursive
flag (if not then issue git submodule update --init --recursive).
Note: We only officially support these assignments on Ubuntu Linux 18.04 (the OS the teaching labs are running) and OSX 10.13 (the OS I use on my personal laptop). While they should work on other operating systems, we make no guarantees.
All grading of assignments is done on Linux 18.04
All assignments will have a similar directory and file layout:
README.md
CMakeLists.txt
main.cpp
assignment_setup.h
include/
function1.h
function2.h
...
src/
function1.cpp
function2.cpp
...
data/
...
...
The README.md file will describe the background, contents and tasks of the
assignment.
The CMakeLists.txt file setups up the cmake build routine for this
assignment.
The main.cpp file will include the headers in the include/ directory and
link to the functions compiled in the src/ directory. This file contains the
main function that is executed when the program is run from the command line.
The include/ directory contains one file for each function that you will
implement as part of the assignment.
The src/ directory contains empty implementations of the functions
specified in the include/ directory. This is where you will implement the
parts of the assignment.
The data/ directory contains sample input data for your program. Keep in
mind you should create your own test data to verify your program as you write
it. It is not necessarily sufficient that your program only works on the given
sample data.
This and all following assignments will follow a typical cmake/make build routine. Starting in this directory, issue:
mkdir build
cd build
cmake ..
If you are using Mac or Linux, then issue:
make
Compiling the code in the above manner will yield working, but very slow executables. To run the code at full speed, you should compile it in release mode. Starting in the build directory, do the following:
cmake .. -DCMAKE_BUILD_TYPE=Release
Followed by:
make
Your code should now run significantly (sometimes as much as ten times) faster.
If you are using Windows, then running cmake .. should have created a Visual Studio solution file
called a3-finite-elements-3d.sln that you can open and build from there. Building the project will generate an .exe file.
Why don't you try this right now?
Once built, you can execute the assignment from inside the build/ using
./a3-finite-elements-3d
In this assignment you will get a chance to implement one of the gold-standard methods for simulating elastic objects -- the finite element method (FEM). Unlike the particles in the previous assignment, the finite-element method allows us compute the motion of continuous volumes of material. This is enabled by assuming that the motion of a small region can be well approximated by a simple function space. Using this assumption we will see how to generate and solve the equations of motion.
FEM has wonderfully practical origins, it was created by engineers to study complex aerodynamical and elastic problems in the 1940s. My MSc supervisor used to regale me with stories of solving finite element equations by hand on a chalkboard. With the advent of modern computers, its use as skyrocketed.
FEM has two main advantages over mass-spring systems. First, the behaviour of the simulated object is less dependent on the topology of the simulation mesh. Second, unlike the single stiffness parameter afforded by mass spring systems, FEM allows us to use a richer class of material models that better represent real-world materials.
Part I of this SIGGRAPH Course, by Eftychios Sifakis and Jernej Barbic, is an excellent source of additional information and insight, beyond what you will find below.
The idea of the finite element method is to represent quantities inside a volume of space using a set of scalar basis or shape functions 
where 
is a point inside the space volume. We then represent any quantity inside the volume as a linear combination of these basis functions:
where 
are weighting coefficients. Designing a finite element method involves making a judicious choice of basis functions such that we can compute the 's efficiently. Spoiler Alert: in the case of elastodynamics, these 's will become our generalized coordinates and will be computed via time integration.
For this assignment we will use a tetrahedron as the basic space volume. The reason we work with tetrahedra is two-fold. First, as you will see very soon, they allow us to easily define a simple function space over the volume. Second, there is available software to convert arbitrary triangle meshes into tetrahedral meshes.
Now we are getting down to the nitty-gritty -- we are going to define our basis functions. The simplest, useful basis functions we can choose are linear basis functions, so our goal is to define linear functions inside of a tetrahedron. Fortunately such nice basis functions already exist! They are the barycentric coordinates. For a tetrahedron there are four (4) barycentric coordinates, one associated with each vertex. We will choose 
to be the 
barycentric coordinate.
Aside from being linear, barycentric coordinates have another desireable property, called the Kronecker delta property (or fancy Identity matrix as I like to think of it). This is a fancy-pants way of saying that the 
barycentric coordinate is zero (0) when evaluated at any vertex, 
, of the tetrahedron, 
, and one (1) when evaluated at 
. What's the practical implication of this? Well it means that if I knew my function 
, then the best values for my 's would be 
, or the value of 
evaluated at each vertex of my tetrahedron.
All of this means that a reasonable way to approximate any function in our tetrahedron is to use
where 
are now the tetrahedron barycentric coordinates and 
are the values of at the nodes of the tetrahedron. Because our basis functions are linear, and the weighted sum of linear functions is still linear, this means that we are representing our function using a linear function space.
To apply this idea to physics-based animation of wiggly bunnies we need to more clearly define some of the terms above. First, we need to be specific about what our function will be. As with the particles in the previous assignments, what we care about tracking is the position of each mesh vertex, in the world, over time. For the vertex we can denote this as 
. We are going to think of this value as a mapping from some undeformed space 
into the real-world. So the function we want to approximate is 
which, using the above, is given by
The take home message is that, because we evaluate 's in the undeformed space, we need our tetrahedron to be embedded in this space.
Now that we have our discrete structure setup, we can start "turning the crank" to produce our physics simulator. A single tetrahedron has four (4) vertices. Each vertex has a single 
associated with it. As was done in assignment 2, we can store these nodal positions as a stacked vector and use them as generalized coordinates, so we have
Now let's consider the velocity of a point in our tetrahedron. Given some specific 
, the velocity at that point is
However, only the nodal variables actually move in time so we end up with
Now we can rewrite this whole thing as a matrix vector product
where 
is the 
Identity matrix.
Now that we have generalized coordinates and velocities we can start evaluating the energies required to perform physics simulation. The first and, and simplest energy to compute is the kinetic energy. The main difference between the kinetic energy of a mass-spring system and the kinetic energy of an FEM system, is that the FEM system must consider the kinetic energy of every infinitesimal piece of mass inside the tetrahedron.
Let's call an infinitesimal chunk of volume 
. If we know the density 
of whatever our object is made out of, then the mass of that chunk is 
and the kinetic energy, 
is 
. To compute the kinetic energy for the entire tetrahedron, we need to integrate over it's volume so we have
BUT~ 
is constant over the tetrahedron so we can pull that outside the integration leaving
in which the per-element mass matrix, 
, makes an appearance.. In the olden days, people did this integral by hand but now you can use symbolic math packages like Mathematica, Maple or even Matlab to compute its exact value.
Now we need to define the potential energy of our tetrahedron. Like with the spring, we will need a way to measure the deformation of our tetrahedron. Since the definition of length isn't easy to apply for a volumetric object, we will try something else -- we will define a way to characterize the deformation of a small volume of space. Remember that all this work is done to approximate the function which maps a point in the undeformed object space, , to the world, or deformed space. Rather than consider what happens to a point under this mapping, let's consider what happens to a vector.
To do that we pick two arbitary points in the undeformed that are infinitesimally close. We can call them 
and 
(boring names I know). The vector between them is 
. Similarly the vector between their deformed counterparts is 
. Because we chose the undeformed points to be infinitesimally close and 
, we can
use Taylor expansion to arrive at
where 
is called the deformation gradient. Remember, results from differentiating a 
-vector by another -vector so it is a 
matrix.
Because 
is pointing in an arbitrary direction, , captures information about how any changes locally, it encodes volumetric deformation.
The FEM discretization provides us with a concrete formula for 
which can be differentiated to compute . An important thing to keep in mind -- because our particular FEM uses linear basis functions inside of a tetrahedron, the deformation gradient is a constant. Physically this means that all 's are deformed in exactly the same way inside a tetrahedron.
Given we can consider the squared length of any 
Like the spring strain, 
is invariant to rigid motion so it's a pretty good strain measure.
The potential energy function of a tetrahedron is a function that associates a single number to each value of the deformation gradient. Sadly, for the FEM case, things are a little more complicated than just squaring (but thankfully not much).
Like the kinetic energy, we will begin by defining the potential energy on an infinitesimal chunk of the simulated object as 
where 
is called the *strain energy density function. Mostly, we look up strain energy density functions in a book. Material scientists have been developing them for many years so that they mimic the behaviour of realistic materials. For this assignment you will use the well established, Neo-Hookean (its better than Hooke's Law because its new) strain energy density for compressible materials. This model approximates the behaviour of rubber-like materials. There are many formulations on that page, and we'll use the following:
where 
.
The total potential of the tetrahedron can be defined via integration as
Typically we don't evaluate potential energy integrals by hand. They get quite impossible, especially as the FEM basis becomes more complex. To avoid this we typically rely on numerical quadrature. In numerical quadrature we replace an integral with a weighted sum over the domain. We pick some quadrature points 
(specified in barycentric coordinates for tetrahedral meshes) and weights and evaluate
However, for linear FEM, the quadrature rule is exceedingly simple. Recall that linear basis functions imply constant deformation per tetrahedron. That means the strain energy density function is constant over the tetrahedron. Thus the perfect quadrature rule is to choose as any point inside the tetrahedron (I typically use the centroid) and as the volume of the tetrahedron. This is called single point quadrature because it estimates the value of an integral by evaluating the integrated function at a single point.
The per-element generalized forces acting on a single tetrahedron are given by
and the stiffness is given by
These can be directly computed from the quadrature formula above. Again, typically one uses symbolic computer packages to take these derivatives and you are allows (and encouraged) to do that for this assignment.
For a tetrahedron the per-element forces are a 
vector while the per-element stiffness matrix is a dense, 
matrix.
Extending all of the above to objects more complicated than a single tetrahedron is analogous to our previous jump from a single spring to a mass-spring system.
The initial step is to divide the object to be simulated into a collection of tetrahedra. Neighboring tetrahedra share vertices. We now specify the generalized coordinates of this entire mesh as
where 
is the number of vertices in the mesh. We use selection matrices (as we did in assignment 2) which yield identical assembly operations for the global forces, stiffness and mass matrix.
Because you can compute all the necessasry algebraic operators (
, 
, 
) you can use your linearly-implict Euler code from assignment 2 to integrate your FEM system. To see the limitations of this approach, run a3-finite-elements-3d arma and press N. Interacting with the armadillo will almost immediately cause your simulation to explode. This is because our bunny was secretly gigantic! So its effective stiffness was very low. This armadillo is much smaller (less than 1 meter tall). Even though it uses the same material parameters, its effective stiffness is much higher. Linearly-implicit Euler just cannot handle it, and so ... Kaboom!
To fix this you will implement the full backward Euler integration scheme which will solve the discrete time stepping equations
The position of the mesh is updated using 
.
To solve for 
it is useful to notice that solving the update equation above is equivalent to the optimization problem
We are going to solve this minimization problem using Newton's method.
Newton's method computes the local minimum of an objective function by iteratively solving a sequence of quadratic minimizations. Let's look at the process for a single iteration (say the iteration). At this iteration we already have the guess for our new velocity, denoted 
. The goal of the 
iteration is to compute a small change in the current velocity estimate,
, that reduces the objective function.
This is done by solving a local, quadratic optimization of the formula
where 
is the Gradient of the backward Euler integration cost function, and 
is the Hessian. In this case we have 
and 
The solution to this problem is given by 
. You then compute 
.
Repeating the process of constructing and solving this quadratic optimization is at the heart of Newton's method.
Unfortunately the computed by Newton's method can be overly ambitious causing the algorithm to never find a local minimum. To avoid this problem, robust optimization schemes give themselves the option of taking a fractional newton's step. One simple way to find a good step size is to use backtracking line search. Line search is so named because it searches in the direction computed by Newton's method (along the line) but looks for a value of the full energy function that guarantees, for instance, sufficient decrease. Backtracking line search gets its name because it initially tries to take a full Newton Step, and if that step is flawed, divides the step by some ratio and then checks the result of this new step. This procedure is repeated until either a suitable step is found, or the step length being checked goes to zero.
The final component of this assignment involves making our FEM simulation look a bit more appealing. FEM calculations can be slow, especially if we want real-time performance. This often limits us to the use of relatively low resolution simulation meshes. To compensate for this we can adopt the concept of skinning from computer animation.
Let us define two different meshes. The first is our simulation tetrahedral mesh with vertex positions given by 
. The second is going to be a higher resolution triangle mesh, for display purposes. We will assume that our display mesh is completely enclosed by the simulation mesh when then simulation mesh is undeformed.
Our goal is to transfer the motion of the simulation mesh to the display mesh. Remember that the FEM discretization defines the motion of our object, not just at the vertices, but everywhere inside the mesh (via the basis functions). Specifically, for any point in the undeformed space I can reconstruct the deformed position as long as I know which tetrahedron contains .
This gives us a simply algorithm to deform our display mesh. For each vertex in the display mesh (
), find the tetrahedron, 
, that contains then move that vertex to position 
. Here 
and 
are the basis function and generalized coordinates for the containing tetrahedron.
This can be expresses as a linear operation
where 
are the deformed vertex positions of the display mesh and 
is the skinning matrix which contains the appropriate basis function values.
For this course most functions will be implemented in .cpp files. In this assignment the only exception is that time integrators are implemented in .h files. This is due to the use of lambda functions to pass force data to the time integration algorithms. Finite element derivatives are both tricky and tedious to do correctly. Because of this you DO NOT have to take these derivatives by hand (unless you want to show off your mad skills). You can use a symbolic math package such as Maple, Mathematica or Matlab. Note: You can not use automatic differentiation, only symbolic math packages.
Running a3-finite-elements-3d will show a coarse bunny mesh integrated with linearly-implicit Euler. Running a3-finite-elements-3d arma will show a coarse armadillo mesh integrated using backward Euler. For each mesh you can toggle between the simulation and skinned meshes by pressing S. You can toggle integrators by pressing N. Note: the linearly-implicit integrator WILL explode when used with the armadillo mesh.
Some of the functions from assignment 2 are reused here. If you correct implementation errors in those functions, your grades for the previous assignment will be updated to reflect that.
Evaluate the linear shape functions for a tetrahedron. This function returns a 4D vector which contains the values of the shape functions for each vertex at the world space point x (assumed to be inside the element).
Piecewise constant gradient matrix for linear shape functions. Row 
of the returned matrix contains the gradient of the shape function.
Compute the Neo-Hookean strain energy density.
Compute the first Piola-Kirchoff (the gradient of the strain energy density with respect to the deformation gradient). You can use a symbolic math package to do this (but don't just look it up on the internet please).
Compute the hessian of the strain energy density with respect to the deformation gradient. You can use a symbolic math package to do this (but don't just look it up on the internet please).
Compute the kinetic energy of a single tetrahedron.
Single point quadrature for a constant strain tetrahedron (CST).
Compute the potential energy of a single tetrahedron. Note: you will need both psi_neo_hookean.cpp and quadrature_single_point.h to do this.
Compute the gradient of the potential energy of a single tetrahedron. Note: you will need both dpsi_neo_hookean_dq.cpp and quadrature_single_point.h to do this.
Compute the hessian of the potential energy of a single tetrahedron. Note: you will need both d2psi_neo_hookean_dq2.cpp and quadrature_single_point.h to do this.
The potential energy of a non-zero rest length spring attached to two vertices of the mesh. Use your code from the last assignment.
The gradient of the spring potential energy. Use your code from the last assignment.
Compute the dense mass matrix for a single tetrahedron.
Assemble the full mass matrix for the entire tetrahedral mesh.
Assemble the global force vector for the finite element mesh.
Assemble the global stiffness matrix for the finite element mesh.
Build the skinning matrix that maps position from the coarse simulation mesh to the high resolution rendering mesh.
Use your code from the last assignment
Use your code from the last assignment
Use your code from the last assignment
Implement Newton's method with backtracking line search. Use the following parameter values: alpha (initial step length) = 1, p (scaling factor) = 0.5, c (ensure sufficient decrease) = 1e-8.
Using your Newton's method, implement a fully implicit solver. To ensure reasonable performance, use a maximum of five (5) iterations.