Assignment 5 - Cloth Simulation

Due Wednesday, 11/4 Monday, 11/9 at 23:59:59. You must work individually.


In this assignment, we’re going to simulate cloth with linearly implicit Euler, penalty collision forces, and sparse matrices.

Starting Point

Please download the assignment and go over the code.

Your main task is to complete the function Cloth::step(). I have written a lot of code for you that deal with setting up the scene and rendering the cloth. Study the code carefully to see what is happening.

At each time step, you must fill in \(M, K, v\), and \(f\) with appropriate values, and then solve the linear system

\[ A x = b, \]


\[ \begin{aligned} A &= M - h^2 K\\ b &= M v + h f. \end{aligned} \]

The solution, \(x\), of the linear system contains the new velocities of the particles.

Setting up Eigen

For solving a linear system, we are using a C++ library called Eigen. You’ll have to install this in the same way you installed glm in lab 0. You’ll need to create a new environment variable called EIGEN3_INCLUDE_DIR that points to the Eigen directory. Once you have it installed, you may want to look at this quick tutorial.

Explicit Integration

We’ll start with symplectic Euler, which boils down to the same thing as above but with \(A=M\), which means you don’t need to form \(K\) yet.

  1. In Cloth::Step(), loop through the particles to fill in \(M\) and \(v\), and fill in \(f\) with gravity. (Note the const Vector3d &grav argument of Cloth::step().) Use the block() and segment() functions in Eigen. Remember to not fill in these values for fixed particles. (I.e., only add these values for free particles.) The indices of the particles are stored inside the particle class Particle::i. This is the starting index of the particle into \(M,\) \(v,\) and \(f\). For example, if the index of a particle is 6, then you should store its velocity into the \(v\) vector at indices 6, 7, and 8.

  2. In the same function, loop through the springs to add spring forces to \(f\). Remember to add, not set. Each spring object contains pointers to the two particles attached to the spring. The formula for the spring force is

    \[ \begin{aligned} & f_s = E (l-L) \frac{\Delta x}{l},\\ & f = \begin{pmatrix} f_s \\ -f_s \end{pmatrix}, \end{aligned} \]

where \(E\) is the stiffness constant, \(\Delta x = x_1 − x_0\) is the vector between the two particle positions, \(l = \| \Delta x \|\) is the current spring length, and \(L\) is the spring’s rest length. (Each spring is between particle 0 and particle 1.) This force, \(f_s\), is the force for particle p0 of the spring. The force for the other particle p1 is the negation, \(-f_s\). Again, remember to add this force only to free particles.

  1. Solve the linear system \(A x = b\). There are many linear solvers provided by Eigen. The LDLT solver is a good choice for now: A.ldlt().solve(b). The resulting solution \(x\) contains the new velocities of the particles.

  2. Loop through the particles to extract the velocities from the solution vector.

  3. Loop through the particles one last time to integrate the positions using the new velocities.

With symplectic Euler, you should be able to simulate a \(2 \times 2\) system with h=5e-3 and stiffness=1e1. You should see a cloth that looks like a flat sheet that does not bend. Try changing these values - you’ll see that the simulation quickly becomes unstable. One of the best ways to make it stable is to implement implicit integration, which is our next task.

Implicit Integration

Now implement implicit Euler by filling in and using the stiffness matrix, \(K\). The formula for the \(2 \times 2\) block stiffness matrix for a spring is

\[ \begin{aligned} & K_s = \frac{E}{l^2} \left( \left(1 - \frac{l-L}{l} \right) \left( \Delta x \Delta x^T \right) + \frac{l-L}{l} \left( \Delta x^T \Delta x \right) I \right),\\ & K = \begin{pmatrix} -K_s & K_s \\ K_s & -K_s \end{pmatrix}. \end{aligned} \]

Note that \(\Delta x^T \Delta x\) is a scalar, and \(\Delta x \Delta x^T\) is a \(3 \times 3\) matrix. \(I\) is the \(3 \times 3\) identity matrix.

This \(2 \times 2\) block matrix, \(K\), should be added to the global stiffness matrix into the appropriate rows/cols using the indices of the two particles. Again, remember to not add these values for fixed particles. The negated, off-diagonal values should be added only if both particles of the spring are free.

With the following settings in Scene::Scene():

h = 5e-3;
grav << 0.0, -9.8, 0.0;
int rows = 2;
int cols = 2;
double mass = 0.1;          // total mass of the cloth
double stiffness = 1e1;     // E

the \(A\) matrix and the \(b\) vector should be the following for the very first step.

cout << A << endl;
    0.0254         0    0.0001   -0.0003         0         0
         0    0.0250         0         0         0         0
    0.0001         0    0.0254         0         0         0
   -0.0003         0         0    0.0254         0   -0.0001
         0         0         0         0    0.0250         0
         0         0         0   -0.0001         0    0.0254

cout << b << endl;

A common bug is that the sign of the stiffness matrix \(K_s\) is flipped. Try negating it!

Try increasing the cloth resolution to \(10 \times 10\) or higher. You may have to play around with the time step and stiffness values to make the cloth stable.

You should try running your code in release mode (See lab 0) so that you can simulate more particles.

Sphere Collision

If you uncomment line 91 in Scene::draw(), you’ll see a collision sphere moving, but it won’t be affecting the cloth yet. We need to write some collision code first.

Implement collisions in Cloth::step(), right before the linear solve. You must check all the cloth particles against all the collision spheres (which are also particles) in the scene, which in this case is just a single sphere. For each particle, if there is a collision, add a collision force to the RHS and the corresponding implicit term to the LHS. Unlike cloth spring forces, these collision forces act on a single particle.

Let \(x\) be the position of the particle, and \(x_s\) be the position of the sphere. Then the vector from the sphere to the particle is \(\Delta x = x - x_s.\) Let \(l = \| \Delta x \|\) be the length of \(\Delta x.\) If we subtract \(l\) from the sum of the radii of the sphere and the particle, we obtain the penetration depth: \(d = r + r_s - l.\) If this distance, \(d,\) is positive, then we have detected a collision. The collision force should be computed as:

\[ f_c = c d \frac{\Delta x}{l}, \]

where \(c\) is the collision stiffness constant. The stiffness matrix entry can be approximated simply as:

\[ K_c = c d I, \]

where \(I\) is the \(3 \times 3\) identity matrix.

Sparse Matrices (bonus)

If you’re planning to use cloth simulation for your project, I highly recommend converting your code to use sparse matrices instead of dense matrices. On my computer, for a \(10 \times 10\) cloth, sparse is twice as fast as dense. For a bigger system, the difference will be even more significant. See http://eigen.tuxfamily.org/dox-devel/group__TutorialSparse.html.

You have to modify your matrix filling code and matrix solving code. To fill a sparse matrix, you must first create an array of triplets of the form row, col, val. As an example, the code to create a sparse \(2 \times 2\) identity matrix is

typedef Eigen::Triplet<double> T;
std::vector<T> A_;
A_.push_back(T(0, 0, 1));
A_.push_back(T(1, 1, 1));
Eigen::SparseMatrix<double> A(2,2);
A.setFromTriplets(A_.begin(), A_.end());

To solve using a sparse matrix, we must use a specialized sparse solver. There are many options, but for animation purposes, running Conjugate Gradients for several iterations is often sufficient.

ConjugateGradient< SparseMatrix<double> > cg;
VectorXd x = cg.solveWithGuess(b, x0);

For Conjugate Gradients, it is important to give the initial guess, x0. In our case, x represents the new particle velocities, so you should use the previous particle velocities as the initial guess.

Constructing a sparse matrix is expensive. Therefore, it is not efficient to create the mass and stiffness matrices separately and add them to form the final matrix \(A\). Instead, it is better to create \(A = M - h^2 K\) and fill it as you loop through the particles and the springs. You should also create the right-hand-side vector \(b = M v + h f\) as you loop through the particles and springs.

With sparse matrices implemented, try increasing the resolution to \(30 \times 30\).

Point breakdown

Total: 100 points + 15 bonus points

What to hand in

Failing to follow these points may decrease your “general execution” score.

If you’re using Mac/Linux, make sure that your code compiles and runs by typing:

> mkdir build
> cd build
> cmake ..
> make

If you’re using Windows, make sure your code builds using the steps described in Lab 0.

Generated on Mon Oct 26 15:21:39 CDT 2020