Assignment 1: Mass-Spring Systems

**Posted:** Thursday, 2 August 2018

**Updated:** Saturday, 4 August 2018

**Due:** Thursday, 16 August 2018

In this assignment, you will implement a mass-spring system to simulate the dynamics of a sheet of cloth, subject to gravity and other external forces.

Implementing any simulation algorithm involves several components:

- the core physical model: in this case, a collection of particles and springs,
- the time integrator: for example, a procedure implementing the forward Euler scheme, and
- an initialization procedure to set up the system to represent a particular object: in this case, a rectangular arrangement of particles with several types of springs to model a sheet of cloth.

These three components should be cleanly separated in your code. The mass-spring model will store the states of the particles and the parameters of the springs connecting them, and be able to compute the total force on each particle. Keeping the model separate from the time integrator means that you can easily switch between different time integrators whenever you want, without affecting the rest of the code. Similarly, having a separate initialization procedure to create particles and springs means that you can set up different types of scenarios, and even model different kinds of objects other than cloth. See the Witkin and Baraff notes for more about how to implement a time integrator, though you will have to modify their recommended strategy a bit to implement symplectic Euler as mentioned in my time integration slides.

First try implementing a very simple system with a single particle connected by a spring to a fixed point. (You could implement the fixed point as another particle whose position and velocity is simply not updated in the explicit integration step.) Implement a symplectic Euler integrator and run it on this system, and check that the resulting behaviour looks OK.

Verify that your integrator is actually behaving correctly as follows. Initialize the particle with zero initial velocity and no external forces, so that it behaves like a one-dimensional damped harmonic oscillator. Choose a time step \(\Delta t_0\), and run the simulation until a fixed simulation time, say \(t = 1\,\mathrm s\). Report the absolute distance between the computed position and the analytical solution (which you can find using, say, WolframAlpha); this is the error \(e_0\). Repeat the experiment with half the time step, \(\Delta t_1 = \frac12\Delta t_0\), then with half of that, \(\Delta t_2 = \frac12\Delta t_1\), and so on. The sequence of errors \(e_0, e_1, e_2, \dots\) should converge to zero, and the ratio \(e_k/e_{k-1}\) should be approximately \(1/2\) for sufficiently small time steps.

Next, implement integrators for explicit midpoint and fourth-order Runge-Kutta, and verify that they also work correctly. The ratio \(e_k/e_{k-1}\) should tend to \(1/2^2\) and \(1/2^4\) respectively, due to their higher order of accuracy.

Finally, set up the classic test example for cloth simulation: an initially horizontal sheet with two corners pinned to fixed positions, which swings and then back and forth due to gravity. Try to use realistic values in SI units for all the quantities in your simulation. For example, set the acceleration due to gravity to \((0,-9.8,0)\,\mathrm m\,\mathrm s^{-2}\), choose spring lengths so that the total size of the sheet is roughly \(1\,\mathrm m\times1\,\mathrm m\), and choose particle masses to get a reasonable weight for fabric (say \(0.2\,\mathrm{kg}\,\mathrm m^{-2}\)). Spring constants \(k_s\) and \(k_d\) will need to be chosen by trial and error; you will want to use much smaller values for the bending springs than for the others.

Since all three of the above integrators are explicit, you will probably need to take very small time steps to keep the simulation stable. To prevent the animation itself from looking like a slow-motion video, run multiple simulation time steps per frame until the simulated time adds up to the desired time between frames (usually \(1/60\,\mathrm s\) for games or \(1/24\,\mathrm s\) for video). It’s OK if your simulation does not actually run in real time and takes longer than say \(1/24\,\mathrm s\) to compute \(1/24\,\mathrm s\) of animation; you may save the image of each frame and then compile the images into a 24 fps video afterwards.

The specific requirements for the base component of the assignment are as follows:

Implement symplectic Euler, explicit midpoint, and RK4, and verify that they have first-, second-, and fourth-order accuracy as described above. Include the results of the convergence tests in your submitted report.

Set up the example of the rectangular cloth sheet, using at least a \(10\times10\) grid of particles. Submit the video(s) showing results with all three integrators.

Completing the base requirement is worth 60% of the points for this assignment. Each of the following four extensions is worth an additional 20%, so you need to do at least two of them to have the possibility of getting 100% points. Note that the total points for the assignment are capped at 100%. The extensions are independent and not have to be done in sequential order, i.e. you can choose any two or more to implement in your submission, not necessarily the first two.

Your implemented extensions also do not have to be mutually compatible. For example, if your code can do implicit integration without collisions, and can do collision handling with explicit integration, but cannot do implicit integration and collisions simultaneously, you will not lose any points for that.

This is a fairly simple modification to start you off. Wejchert and Haumann’s 1991 paper “Animation Aerodynamics” (posted on Moodle) describes a simple method to model the interaction of objects with wind. In particular, in their section on “Objects in Flows”, they give a formula for normal and tangential forces on a surface in terms of the velocity relative to the wind.

Use their model to add wind forces to your cloth simulator. To apply their formula, you will need to estimate a surface normal at each particle. You can do this by computing two tangential directions along the cloth surface using finite differences, then taking their cross product. Describe how you do this in your submitted report.

Create an animation of a flag blowing in the wind, modeled as a rectangular sheet of cloth with all vertices on one edge constrained. Choose appropriate values for the coefficients \(\alpha^n\) and \(\alpha^t\), and a realistic wind speed in the range of \(5\) to \(15\,\mathrm{m\,s}^{-1}\).

As you’ve probably noticed by now, explicit integration requires small time steps to remain stable, especially if you make the springs strong enough to not look like excessively stretchy cloth. Implement a backward Euler integrator so that you can take large time steps without losing stability.

To do this, you will need to modify your mass-spring model so that it computes not just the spring forces but also their Jacobians with respect to particle positions and velocities. The backward Euler integrator will need these as well as the inertia matrix, which simply contains the particle masses (each one repeated 3 times along the diagonal). Verify that your Jacobian implementations are correct by comparing them with finite differences: for example, compute the forces due to a single spring at some particle positions \(\mathbf x\) and velocities \(\mathbf v\), choose some very small perturbations \(\Delta\mathbf x\) and \(\Delta\mathbf v\) and recompute the forces, and verify that the change is approximately \(\mathbf J_{\mathbf x}\Delta\mathbf x + \mathbf J_{\mathbf v}\Delta\mathbf v\).

For backward Euler, a single Newton step is sufficient for credit in this assignment, so you will only need to perform one linear solve per time step. You have two choices for how to do this.

Construct the inertia matrix and the Jacobians as \(3n\times3n\) matrices, and use a conjugate gradient solver provided by a linear algebra library to do the linear solve. This will end up being very slow and requiring \(O(n^2)\) memory unless you use a sparse matrix format, so you will have to figure out how to do that (see e.g. Eigen’s documentation for sparse matrix manipulation).

Represent the inertia matrix and the Jacobians as functions which, given a vector, return the corresponding matrix-vector product. This can be a bit more intuitive and you don’t have to mess around with sparse matrices, but you may have to implement the conjugate gradient method yourself (see Shewchuk’s Appendix B2).

Describe your implementation choices in your submitted report.

Demonstrate your implementation by creating an animation of a not-too-stretchy cloth simulated with backward Euler at moderately large time steps (at least \(1/60\,\mathrm s\)).

Allow your mass-spring system to react to other objects by adding a simple collision handling scheme. Specifically, implement penalty forces to allow the cloth to collide with simple shapes such as planes, spheres, and boxes. *Self-collision* handling is a lot more complicated, so we will ignore it.

For each “collider” shape that you want the cloth to respond to, you should implement a collision detection function that (i) determines whether a given point \(\mathbf p\) is inside the shape, and (ii) if so, returns the closest point \(\mathbf q\) on the shape’s surface. Then the collision normal can be defined as the unit vector parallel to \(\mathbf q-\mathbf p\). At the beginning of each time step, use this function to check each particle for collision. For each particle that is inside the shape, add a normal force \(\mathbf f_n\) corresponding to a zero-length spring between the particle and the closest surface point. Note that if the shape is moving, the velocity of the surface point will also be nonzero.

If you wish, you may also add a tangential force \(\mathbf f_t\) modeling Coulomb friction with coefficient \(\mu\). This friction force should be opposite to the *relative* tangential velocity of the particle and the surface, and its magnitude should be the minimum of \(\mu\|\mathbf f_n\|\) and the force needed to zero out the tangential relative velocity (so you do not apply excessive friction forces if the particle is already at rest).

Since penalty forces only act when the particle is already inside the other object, you may find it useful to create a collider that is larger than the actual shape that you draw in the animation, so that the visual results do not show undesirable penetrations. Tuning the parameters for the penalty forces will take some trial and error.

For full credit in this component, you must implement collisions with a ground plane, and with at least one object, either a sphere or a box, moving on a scripted trajectory (e.g. with translation \((0,0,\sin t)\)). Describe in your report how you computed closest points, normal forces, and (optionally) tangential forces.

Create an animation showing collisions with both the ground plane and the moving object. See the last few examples under “Animations and Images” on Robert Bridson’s webpage for inspiration.

Many real fabrics show barely any visible stretching. Modeling them as perfectly inextensible and applying a constrained dynamics approach is an efficient way to simulate such materials.

Implement the fast projection algorithm described in Goldenthal et al.’s 2007 paper, “Efficient Simulation of Inextensible Cloth”. This approach fits best with an symplectic Euler integrator. To incorporate this into your mass-spring system, replace the structural springs with inextensibility constraints, integrate the remaining forces over a time step as usual, then perform a projection step that updates the positions and velocities of all the particles. The projection step will require the solution of a sparse linear system. Refer to the discussion of the linear solver in the “Implicit integration” component above for advice.

In your report, give your derivation of the constraint Jacobian (denoted \(\nabla\mathbf C\) by Goldenthal et al.), and describe how you set up and solve the linear system.

Create one or more animations demonstrating inextensibility. Some good examples are the sheet pinned at two corners (compare Goldenthal et al.’s Fig. 1), and a flag under strong wind (Provot’s Fig. 6).

If you are using my C++/OpenGL starter code from Assignment 0, consider switching to **this updated code** which allows you to draw smoothly shaded surfaces using vertex normals, and also includes functions for drawing boxes and spheres to help with the collision handling component of the assignment.

A submission form will be opened on Moodle for this assignment. There, you should submit a zip file containing the following items:

a report in PDF format describing the work you did in this assignment, including the information required by problem descriptions,

a video or set of videos showing the results of your implementation, and

a copy of your implementation code.

Please include in the report a brief description of each of the results shown in your submitted video(s), since it may not be obvious to me what each animation is supposed to demonstrate.