CS 348C Assignment #2: Constrained Dynamics (inextensibility constraints)

Professor: Doug James
Due date: Thursday, Oct 26, 2017 (midnight) 

Inextensible yarn curves from [Kaldor et al. 2010]

In this second assignment, you will extend your assignment #1 code to support hard holonomic (position-based) constraints. Specifically, you will use constraints to implement deformable 2D strands with inextensible (fixed-length) line segments. Why? Well, stretchy spaghetti looks weird, and (as you will see later) exploiting these constraints can speed up other calculations, e.g., of bending forces. You will use your previous collision detection and resolution implementation, and should be able to re-run your prior examples with better length preservation.

Background: As discussed in class, inextensibility constraints will be added for each edge to maintain a fixed length equal to its rest length. For example, we will use constraints to enforce that the squared-length of an edge's inter-particle distance is constant. The formulae for constraints, Jacobian, and constraint matrix system were covered in class.

Implicit constraint direction: There are many ways to discretize and enforce constraints. In this assignment, you will use an implicit velocity-level formulation that seeks to apply an impulse to our predictor velocities (from symplectic Euler) such that the constraint will be enforced at the end of the timestep. The details of this impulse computation have been described in class (Tues Oct 17), and are elaborated on further below; for more details on ICD for cloth simulation, see [Goldenthal et al. 2007].

Combining collisions with hard constraints: As they say, when an immovable object hits an immovable wall, something has to give. Unless you simultaneously resolve both the hard constraints and the collision constraints, they will not both be satisfied. Nevertheless, we can alternate between the two classes of constraints, to improve their satisfaction.  In your implementation, you will (a) solve all inextensibility constraints simultaneously using a matrix solve (details below), and (b) perform the velocity update due to the inextensibility constraint just before the main collision processing loop, but after the velocities have been updated with other forces (gravity, springs, penalty forces).  The inextensibility constraints will be satisfied for all edges, unless they are adjacent to vertices involved in contact impulses. Once you get this working, you can further iterate on both inextensibility and collision constraints (using velocity-level impulse corrections for both), however, this may not be necessary due to our ICD formulation.

How to demonstrate the merits of your simulator:
  1. Simulate a single very long strand of several dozen particles (similar to the one on the A1 webpage).
  2. Report the maximum relative length error (|delta Length|/|length|) for all constrained edges for all time, e.g., as the strand falls on the bottom of the container.
  3. After you've done (1) & (2), re-run the spaghetti factory to demonstrate the effect of inextensibility. There is no need to run for a very long time, however you must identify all strands and enforce constraints on each one. Since meatballs create a constraint loop, they are not handled by the TDMA, and therefore you may simply disable the MeatballFactory for this assignment.
Matrix Solver Details: To greatly simplify your your implementation, you need not support (1) edges with only one pinned vertex, or (2) branching edge structures. In other words, you only scan to find strands in free-flight motion. For each strand, the constraint matrix is a positive definite, tridiagonal matrix, which can be solved efficiently using a linear-time solver such as a classical tridiagonal matrix algorithm (TDMA), such as the Thomas algorithm (wikipedia, code, javacode). Your implementation should identify the vertices in each separate strand (ideally once, at the beginning of the simulation), then assemble and solve each linear system using a TDMA, the resulting impulses are then used to update the particle velocities prior to impulse-based collision processing. (Due to repeated reuse, you should be able to optimize your code to avoid creating/destroying arrays, and simply track strands. Note that the associated constraint matrix and RHS are position dependent, and therefore likely to change at each timestep.)

ADDITIONAL THINGS TO TRY (Optional): Only if you have time, or for a future (final) project, here are some logical next things to explore:
Hand-in using Canvas:  
Start early. Ask questions. Have fun!!! 

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 (or with your partner), and understand what you are writing. Please also list the names of everyone that you discussed the assignment details with.  (You are expected to maintain the utmost level of academic integrity in the course.)