CS 567 Assignment #2: Implicit integration of stiff dynamical systems
Professor: Doug James

Due date: Friday, March 2 (before midnight).

Videos: Please submit videos (described below) that demonstrate the requested (and any additional) functionality of your system. Details on video generation are here.

In this assignment you will add implicit integration to your mass-spring particle system, modifying your first assignment code appropriately.  For simplicity, you are only required to integrate a single undamped hair model implicitly. Using your implementation, you will be able to simulate a long strand of stiff hair falling on the ground and colliding with the walls of the cell, as well as interact with it interactively (include a final video of this).

Assignments Steps: Here are the steps you need to address:

1. Procedural hair generation:  Automate the construction of a hair primitive parameterized by the number of particles, N--so you don't have to draw it by hand each time.  Make it so that your simulator can easily create a long hair with a moderate N value, e.g., N=50-500.  With stiff spring forces, you will notice that your explicit integrator becomes quite inefficient and/or unstable for such models.
   
2. Backward Euler Integrator: You will integrate the hair's stiff extension and bending spring forces using backward Euler as in the Baraff course notes (Implicit Methods for Differential Equations).  The equation that you will need to set-up and solve, hereafter called A*dv=b, is given in (4-6)--note that f in (4-6) is actually the particle acceleration vector. For simplicity we will ignore spring damping forces, and therefore the acceleration function f will only depend on the particle positions, x. Therefore you will only need to compute the function, f0, and it's sparse Jacobian matrix, df/dx. How and where you accumulate the Jacobian is something you should think about.
   
3. Spring force derivatives:  The first step is to add sparse derivative support to the Force objects used to model a single hair (SpringForce2Particle, SpringForceBending).  You can modify the Particle object to accumulate sparse force Jacobians as well forces (.f), keeping in mind that the Jacobian is a sparse matrix with nonzero 2x2 blocks. You could store/access this sparse block row vector in a variety of ways, including just using a HashMap<Integer,Matrix2d> accumulator--keep in mind that the Jacobians are symmetric. (If you go this route, you'll probably want to make a helper object that handles get/set/acc/multiply operations involving the indexed Matrix2d information and global double[] vectors.) You could then modify the Force interface to optionally support apply_forceAndJacobian() to accumulate both force and Jacobian information. Alternately, you can accumulate the sparse derivatives in a generic sparse matrix container from an external library, e.g., MTJ. It's really up to you.
4. Conjugate Gradients: Once you can compute the quantities in equation (4-6), then you can solve the system, A*dv=b. Since A is symmetric and positive definite we will use the Conjugate Gradient (CG) iterative method.
   
• Note that it is not required to formally assemble the matrices df/dx, or A=(I -  h^2*df/dx) in a sparse matrix container in order to solve the matrix equation: you only need to evaluate b=h(f0 + h*df/dx*v0), then provide CG with a way to evaluate matrix-vector products Az, for some vector z.  For example, here is a simple java interface for a matrix-vector multiply function that could be implemented by some object (that uses your ParticleSystem to appropriately accumulate force/Jacobian information, and then repeatedly/quickly multiply by it):
• public interface MatrixMultiplyInPlace {
      /**
       * Multiplies a double[] and puts the result in a preallocated double[].
       * @param x  Vector to be multiplied by matrix.
       * @param Ax Matrix-vector multiply result.
       */
      public void multiply(double[] x, double[] Ax);
  }
• Note that you can implement this interface to multiply by the total spring force df/dx, as well as to implement A by composition.
  
• Feel free to use any implementation of the CG algorithm that you would like, however it is quite easy and instructive to "roll your own":  simply type in the implementation from Shewchuk's paper (page 50... it's not very painful ;). To manipulate double[] vectors you'll need some kind of simple java BLAS (add/acc/dot/axpby/twoNorm), and here is a class (djames.matrix.DoubleVector) which has functions you can borrow.
• Accuracy is something that you can adjust to make your simulation go faster, but be careful since for large relative errors you'll start to see artifacts. For interactive applications, people often only take a fixed number of iterations per time-step to avoid big delays, but this isn't guaranteed to be accurate.  Another problem is that the A matrix may not be positive for large time-steps (doh!), and so CG may fail to converge.
  
5. Preconditioned CG (PCG):  In practice, CG can converge a lot faster if you use almost any reasonable preconditioner (see Shewchuk, page 51--it's only ~two extra lines which isn't too painful). For example, just using the inverted block diagonal can be a big help. Try to use a simple preconditioner, but it is not required. You can have it implement  MatrixMultiplyInPlace for use in PCG (reverting to CG if the preconditioner is null).
6. Additional external forces: You can add additional forces to the system, but you need not integrate them implicitly (phew!). For example, user interactions can be modeled using an explicitly integrated spring force, as can gravitational forces, e.g., using forward Euler. Just remember which Force objects are explicit and implicit in your "advanceTime()" method.
   
7. Collisions and other constraints: Collisions and other constraints, e.g., Pin constraints, will be handled using a velocity-level filter, as (will be) described in class. You should be able to pin your models, as well as have them collide with the walls of the unit computational cell.
   
8. Create Artifacts: See if you can do something interesting. For example, can you model a head of hair?  a furry beast?  some long grass or reeds?  Let your imagination run wild.
9. Videos: Please include videos illustrating the behavior of your system for a single long strand of hair, and if possible its interactive performance for a moderately sized system, e.g., N=100----how big a strand can you do at 10 FPS?
   

On collaboration and academic integrity: You are allowed to collaborate on the assignments to the extent of formulating ideas as a group, and derivation of physical equations. However, you must conduct your programming and write up completely on your own, and understand what you are writing. Please also list the names of everyone that you discussed the assignment with.  (You are expected to maintain the utmost level of academic integrity in the course. Any violation of the code of academic integrity will be penalized severely.)

Hand-in using CMS: Please submit a short write-up (detailing what you did, your findings, and who you discussed the assignment with, etc.), as well as your Java implementation, and any artifacts, videos, etc. Submit videos in a compressed format.

Enjoy!!

Copyright Doug James, February 2007.