## 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:
• Constraints: Support pinned-vertex constraints in your constraint solver. You can do this by adding additional constraints, in more than one way.
• Constraints: Support branching or looping structures, e.g., meatballs. This will require a better matrix solver, which motivates...
• Numerical Linear Algebra: Try a different matrix solver, such as Cholesky factorization (direct solver) or the preconditioned conjugate gradient (PCG) method (iterative solver).
• Numerical Linear Algebra: Instead of a direct matrix solver, try an iterative relaxation scheme (such as Jacobi or Gauss-Seidel) that essentially relaxes just one constraint at a time. Iterative until you reach a fixed number of iterations, or achieve a desired relative-error tolerance.
• Constraint Drift: How well are your inextensibility constraints being satisfied? Do you notice a difference if you enable/disable the edge spring forces from A1?
• Modeling: Try isometric bending forces that exploit the fact that edge lengths are constant when deriving the bending forces. Note this 2D optimization is analogous to the 3D cloth bending forces for which isometric deformations (with inextensible triangle edges) lead to fast quadratic bending forces [Bergou et al. 2006].  Note that when the inextensibility constraint is satisfied, there should be no difference in your bending forces except that they are much faster. One further benefit is that an implicit integrator is more easily implemented given the simpler derivative formula.
• Modeling: Simulate a furry creature. You can ignore collisions for real-time speed, but you will need to fix the end position and orientation of each hair.
Hand-in using Canvas:
• Document: A very brief written document (in txt or PDF) that details what you did, your findings, any problems, and who you discussed the assignment with.
• Include any comparisons or analysis/results that demonstrate the functioning behavior of your implementation---primarily the "largest relative length error" information.
• Citation:
• Include the names of people that you discussed the assignment details with.
• Include citations of external resources (books, webpages, etc.) that you used in your R&D.
• Code: Documented software implementation derived from the Java starter code, or based on an analogous framework in another language, e.g., C++.
• Results: Your creative simulation artifacts, videos, images, etc.
• We may run your software to evaluate your implementation, and/or ask you to provide a demo, wowever please submit video that demonstrates the requested features.
• Submit as a group: If you are working with a partner, be sure to form a two-person group, and submit your zip file as a group.
• Late submission: Keep in mind this is a short assignment which should take less time than other assignments. Please submit on time so that we can continue to the next assignment.
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.)

References: