University of Utah School of Computing
CS 5967/6967:
Physical Simulation for Computer Animation
INSTRUCTOR:    Adam Bargteil (Office hours: By appointment, MEB ???)
WEB PAGE: http://www.eng.utah.edu/~cs6967/
TIME: M/W 1:25-2:45
PLACE:     WEB 122
UNITS: 3

Assignment 2

For this assignment you will implement a basic deformable body or cloth simulator. Though you are given the option here, I strongly encourage you to do deformable bodies unless you have a good reason to prefer cloth.

There are a few reasonable options:
  • A standard non-linear strain finite element simulation as in O'Brien et al. [1999]. This approach will require painfully small timesteps.
  • Use implicit integration as in Baraff & Witkin [1998]. This approach requires the construction and solution of a large linear system.
  • Use linear finite elements and "stiffness warping" as in Mueller et al. [2002] (see also Mueller & Gross [2004] for a better element-based approach). This will also require solving a large linear system, but the system will remain constant and the method will be lightning fast (i.e. interactive).
  • Use invertible finite elements as in Irving et al. [2004] (see also Teran et al. [2005] for a few extra details). This will be somewhat more robust and stable than the standard approach, but not as fast and stable as the implicit options (unless you make this one implicit, which is kind of tricky).
Just to help get you started here is an overview of the basic non-linear simulator.
  1. Initialization
    1. Read input file and create lists of nodes and elements. Setting positions of the nodes and indices of each elements nodes.
    2. Loop through the nodes and set their mass and velocity to zero
    3. Loop through the elements and compute
      • volume
      • Beta
      • Sum of area weighted normals adjacent to each node
    4. Add 1/4 of the mass (volume * density) of the element to each of the element's nodes.
  2. For (time=0; time<end_of_time; time+=dt)
    1. Loop over the nodes and set force accumulator to zero (this is also a fine time to add external forces, e.g. gravity)
    2. Loop over the elements
      1. Compute X from the world positions of the element's nodes
      2. Compute F = X * Beta
      3. Compute Green's strain epsilon = F^T*F - I
      4. Compute 2nd Piola-Kirchhoff stress S = lambda * epsilon_kk * delta_ij + 2 * mu * epsilon_ij
      5. Compute 1st Piola-Kirchhoff stress P = FS
      6. Compute forces for each node in the element f_i = P * 1/3 * sum of area weighted normals (in reference configuration) of faces incident to node i.
      7. Add forces to each node's force accumulator
    3. Loop over nodes and integrate forward in time by updating positions and velocities with an Euler step.
    4. Handle Collisions
      1. Loop over the nodes checking to see if the node has collided with anything. I suggest starting with just a plane and ignoring self-collisions at first.
      2. If there has been a collision
        1. Project the node onto the surface of the obstacle
        2. Set the normal component of its velocity to the velocity of the obstacle (e.g. zero for a stationary obstacle). Don't do this if the node was already moving away from the obstacle (unlikely).
        3. Apply friction
Implicit integration will require changes to step 2.3. Stiffness warping will require precomputing K_e for each element during initialization and computing best fitting rotations for each element at each timestep (by performing a polar decomposition). Additionally, forces are computed through a few matrix multiplies (f_e = R_e * K_e *(R_e^-1 *x - x_0)) and F, strain, and stress need never be explicitly computed. Invertible finite elements will require an eigen-decomposition for each element with stress related to the diagonal component of F.

If you choose to go the implicit route, a simple linear solver can be found here and some other useful code (including 3x3 matrix and sparse matrix classes) is here.
Here are some meshes to get you started. A line that begins with a "v" gives the 3D position of a node in the mesh. A line that begins with a "t" gives the indices (starting from zero) of the nodes of the element. The original star mesh had problems, here is a fixed version. The bunny may still have problems. I've just added some new meshes. If you are frustrated by small, slow timesteps, I recommend trying some of these meshes. The new meshes were generated with NetGen and should all be high quality. I used this code (pipe the .vol file to stdin in and pipe the output to the .mesh file, e.g. "netgenToMesh < foo.vol > foo.mesh") to convert from NetGen's format to my usual .mesh format. Feel free to make your own models with netgen.

Example Meshes

Rendering

Its a good idea to render just the surface mesh. Here is some sample code to extract the surface mesh. Its ripped out of my source file, so it will probably not "just work." There is also code to read the input file here.

Debugging Hints

Take a look at the various quantities you're calculating. If your object is just floating in the air (i.e. no collisions, no external forces), F should be nearly the identity. The forces should be nearly zero. Etc. If you apply gravity, the forces should be equal to gravity. If you're simulation is exploding, try decreasing the timestep. Additionally, try applying known deformations, e.g. multiply all the positions by 2 and see that your object shrinks back to its rest state.