Home

Assignment 5 — Implicit Cloth

Due Thursday, 11/9 at 23:59:59. You must work individually.

Goals

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

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.

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 deals with setting up the scene and rendering the cloth. Study the code carefully to see what is happening.

The most important function is Cloth::step(), which you will implement later in this assignment. 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, \] where \[ \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 the Cloth

First, add some code into Cloth::Cloth() to set up the particles and the springs. Currently, the cloth is hard-coded to be 2x2, with only horizontal (x) and vertical (y) springs.

Explicit Integration

We’ll start with symplectic Euler, which boils down to solving \(Ax=b\) with \(A\) set to just \(M\). This 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<3,3>(i,j) and segment<3>(i) 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 to \(f\), rather than setting the values of \(f\). (In other words, we want += rather than =.) 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.

  3. 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.

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

  5. 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 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;
             0
       -0.0012
             0
             0
       -0.0012
             0

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 run your code in release mode (See lab 0) so that you can simulate more particles. Press the t key while simulating to show the amount of time the step function takes.

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.

To compute the collision normal, we use the positions of the particles and the sphere. Let \(\vec{x}\) be the position of the particle, and \(\vec{x}_s\) be the position of the sphere. Then the vector from the sphere to the particle is \(\Delta\vec{x} = \vec{x} - \vec{x}_s\). Let \(l = \| \Delta\vec{x} \|\) be the length of \(\Delta x.\) The collision normal is then \(\vec{n} = \frac{\Delta\vec{x}}{l}\). To obtain the penetration depth, we subtract \(l\) from the sum of the radii of the sphere and the particle: \(d = r + r_s - l\). If this distance, \(d,\) is positive, then we have detected a collision. The collision force is computed using the depth and the normal: \[ f_c = c d \vec{n}, \] 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

Convert your code to use sparse matrices instead of dense matrices. On my computer, for a \(30 \times 30\) cloth, sparse is noticeably faster than dense. For a bigger system, the difference will be extremely significant. See Eigen’s sparse matrix tutorial.

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());

Not everything needs to be sparsified.

To solve using a sparse matrix, we must use a specialized sparse solver. There are many options, but for animation purposes, running Conjugate Gradients (CG) for several iterations is often sufficient. The following code solves \(A\vec{x} = \vec{b}\) with CG.

ConjugateGradient< SparseMatrix<double> > cg;
cg.setMaxIterations(25);
cg.setTolerance(1e-6);
cg.compute(A);
VectorXd x = cg.solveWithGuess(b, x0);

For CG, 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.

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

Parameter Exploration

Try changing the stiffness and time step parameters. Write a short paragraph in the README discussing the following:

Point breakdown

Total: 100 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
> ./A5 <SHADER DIR>

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

Generated on Sat Nov 4 17:16:54 CDT 2023