index.php

Main navigation | Main content

These are ungraded assignments intended to help you gain experience with implementing simulation algorithms. Most of them will build on top of the homeworks from previous classes, so that after a while you will have written a moderately sophisticated system that does interesting things.

The starter code takes care of opening a window, creating a grid-shaped mass-spring system, and drawing the masses and springs on the window. To compile it, you will need the FreeGLUT and Eigen libraries.

You are encouraged to look through the starter code even if you are already familiar with OpenGL or you want to use a different programming language, because it provides a simple skeleton for structuring your code. It also includes an example implementation of viscous drag forces, support for dragging points around with your mouse, and examples of reasonable parameter values to start with.

What you should do:

First of all, download the code and make sure you can compile and run it. You should see a window showing a motionless grid of particles.

The starter code has definitions for particles, various forces, and explicit time integration, but most of the important functions are blank. Probably the first things you should implement are the functions related to time integration, i.e.

`ParticleSystem::step`

,`ParticleSystem::getState`

, etc. in`particles.cpp`

, and`explicitEulerStep`

in`integration.cpp`

. Then you should be able to drag individual points around using your mouse, although no spring forces will be acting.Next, implement

`SpringForce::addForces`

, for example according to Witkin and Baraff. (I forgot to add a variable for the damping constant in the`SpringForce`

class, but go ahead and add it in.)The structure should tend to collapse as soon as you pull on it, because there are no shear springs. Change

`makeGrid`

to add springs from particle (*i*,*j*) to particles (*i*+1,*j*+1) and (*i*-1,*j*+1) (as long as they exist). Make sure you give them the right rest length. Should they have the same spring constant as the existing (“structural”) springs?Experiment! Here are some ideas you could try:

- For example, you could see what happens when you changing various parameters, like the resolution of the mass-spring system, the stiffness of the springs, the drag coefficient, and the time step and number of substeps.
- You could add a gravity force that pulls all the points in the negative
*y*direction. - Once you have gravity, you may want to fix two points using an
`AnchorForce`

so that it behaves like a cloth on a washing line. - Or add a
`GroundForce`

that acts upwards on any point whose*y*-coordinate is too low, to prevent the system from falling off the bottom of the screen. - Try simulating a rope or a hair strand by replacing
`makeGrid`

with a function that creates a one-dimensional chain of particles connected to their immediate neighbours. Then add bending springs connecting particle*i*to particle*i*+2; these should be weaker than the structural springs.

Don’t worry if all this seems like a lot! Spend a few hours on it and do however much you can, and in Tuesday’s class we can talk about how far everybody got and what problems you ran into if any.

This one should hopefully compile on Windows, Mac, and Linux. For Windows, some students have provided helpful compilation instructions.

Don’t extract it on top of your existing code because you’ll overwrite any changes you’ve made! Instead extract it somewhere else and copy over the functions you’ve implemented.

Other minor changes in this version:

- Removed
`drawing.cpp`

, moved drawing functions to`particles.cpp`

- Added damping to
`SpringForce`

and`AnchorForce`

- Mouse only grabs particle when it is within 0.2 units
- Changed spring constant in
`makeGrid`

and`mouseForce`

If you find that the starter code doesn’t work for any reason, please send me an email as soon as possible so that I can fix it.

If you’re confused about how to implement the time integration functions, refer to section 4 of “Differential Equation Basics” and section 3 of “Particle System Dynamics” in the Witkin and Baraff notes.

Starter code for implicit integration

Sorry, it turns out that the generic integration framework in the previous code is not quite generic enough to support implicit integration like backward Euler efficiently. Please switch to the new code above.

**Info about the new code:**

The new starter code replaces the

`ODESystem`

base class with`PhysicalSystem`

, which exposes more information like the mass matrix, forces, and Jacobians, which are necessary to perform the backward Euler method as discussed in class. For your convenience, these methods are already implemented in`ParticleSystem`

, but you should look through them and make sure you understand what they are doing.We’ll represent large vectors and matrices with Eigen’s

`VectorXd`

and`MatrixXd`

, which are vectors and matrices of dynamic size allocated at runtime. You can use`VectorXd.segment()`

and`MatrixXd.block()`

to access subvectors/submatrices within them. (See the Eigen quick reference and the example implementations in the code.)Each particle now knows its own index

`i`

, so that if you only have a pointer to the particle you can still write its forces and Jacobians into the right places in the vector/matrix. (See`AnchorForce`

for an example.)The starter code uses dense matrices because they’re easier to work with. Of course, this also means that the code will not be as fast as it would be with careful use of sparse matrices.

**What you should do:**

Copy over the shear springs from your previous implementation of

`makeGrid()`

. Also copy over your implementation of`SpringForce::addForces()`

and make it work in the new style. (See`AnchorForce::addForces()`

for reference.)Implement

`backwardEulerStep()`

and change`ParticleSystem::step()`

to use it. In`backwardEulerStep()`

, you will have to…- obtain the state, mass matrix, forces, and Jacobians from the particle system,
- compute the matrix \(\mathbf A = \mathbf M - \mathbf J_{\mathbf x}\Delta t^2 - \mathbf J_{\mathbf v}\Delta t\) and the vector \(\mathbf b = (\mathbf f^0 + \mathbf J_{\mathbf x}\mathbf v^0\Delta t)\Delta t\),
- solve \(\mathbf A\Delta\mathbf v = \mathbf b\) using the provided
`solve()`

function to obtain \(\Delta\mathbf v\), - and use it to compute the new state \(\mathbf v^1 = \mathbf v^0 + \Delta\mathbf v\) and \(\mathbf x^1 = \mathbf x^0 + \mathbf v^1\Delta t\).

Without the Jacobians for the spring forces implemented, this behaves mostly the same as forward Euler, just slower. (Convince yourself that if

*all*the Jacobians are take as zero, it’s*exactly*the same as forward Euler.) Implement the Jacobians for the spring forces using the formulas in the slides.The reason it seems incredibly slow now is that we’re actually taking 50 time steps per displayed frame (because that’s how small a time step explicit integration needed). Find the

`idle()`

function in`main.cpp`

and change the number of substeps to 1. With implicit integration, each individual time step is more expensive, but you need far fewer of them.One last thing: The backward Euler method is guaranteed to be stable if the Jacobians are negative definite, but when the factor \((1 - \ell/\|\mathbf x_{ij}\|)\) in the spring force Jacobian is negative (i.e. the spring is compressed) this may no longer be the case. One simple hack to deal with this is to replace it with \(\max(1 - \ell/\|\mathbf x_{ij}\|, 0)\).

That should be it! Now you can make any or all of your stiffness and damping coefficients arbitrarily large, and the system should remain unconditionally stable.

Homeworks are back! Here’s a reference implementation of all the homework so far:

This time we will use constraints to model a perfectly inextensible rope, that is, one that cannot be stretched or compressed at all.

Starter code for constrained dynamics

In the starter code, I’ve modeled a rope using a chain of particles connected by springs, with one end attached to a fixed position using an anchor spring. Unfortunately, this makes it really stretchy. We could try increasing the stiffness and switching to implicit integration, but the motion becomes undesirably damped (try it!). Instead, we will replace all the forces with hard constraints, so that the rope becomes inextensible and the end becomes rigidly tied down.

In the set-up function

`makeChain()`

, replace`AnchorForce`

and`SpringForce`

with`PositionConstraint`

and`DistanceConstraint`

respectively. You will have to create*two*instances of`PositionConstraint`

to fix both the*x*and*y*coordinates.Implement

`ParticleSystem::getConstraintValues()`

and`ParticleSystem::getConstraintJacobian()`

. Then in`ParticleSystem::step()`

, change`forwardEulerStep`

to`constrainedForwardEulerStep`

and reduce`substeps`

to 5.If you run the simulation now, the whole chain will fall down because the time integration ignores the constraints. Implement

`projectPositions()`

and`projectVelocities()`

in`integration.cpp`

following the instructions in the comments. Now the one vertex whose position is constrained will remain in place. (For nonlinear constraints,`projectPositions()`

amounts to one iteration of Newton’s method.)Implement

`DistanceConstraint::getValue()`

and`DistanceConstraint::addGradient()`

. You’ll have to work out the gradient of \(\|\mathbf x_1-\mathbf x_0\|-\ell\) by hand.

Now you should have a pretty sturdy length of rope to play with.

**Note:** Our implementation is a lot slower than it could be, because we’re using dense matrices for everything even though most of them are pretty sparse. Working with sparse matrices in Eigen is a lot less straightforward, though, so I decided not to complicate the exercise. If you want to try changing the code to use sparse matrices wherever possible, I’ll leave that for you as an optional exercise.

**An explanation for what \((\mathbf J\mathbf M^{-1}\mathbf J^T)^{-1}\) means and why it shows up:** Suppose we put everything else on hold and just allow the constraints to act for 1 unit of (fictitious) time. Each constraint \(i\) applies a force parallel to its gradient, \(\nabla c_i\), and proportional to what’s known as its Lagrange multiplier, \(\lambda_i\). The total force is therefore \[\begin{align}
\mathbf f &= \lambda_1\nabla c_1+\lambda_2\nabla c_2+\cdots+\lambda_m\nabla c_m \\
&= \begin{bmatrix}\nabla c_1&\nabla c_2&\cdots&\nabla c_m\end{bmatrix}\begin{bmatrix}\lambda_1\\\lambda_2\\\vdots\\\lambda_m\end{bmatrix} \\
&= \mathbf J^T\boldsymbol\lambda.
\end{align}\] Note that \(\mathbf J\) is an \(m\times n\) matrix, where \(m\) is the number of constraints and \(n\) is the number of degrees of freedom.

Over unit time the force \(\mathbf f=\mathbf J^T\boldsymbol\lambda\) produces a change in velocity \(\Delta\mathbf v=\mathbf M^{-1}\mathbf J^T\boldsymbol\lambda\), and let’s say the associated displacement is \(\Delta\mathbf x=\mathbf M^{-1}\mathbf J^T\boldsymbol\lambda\) as well. (If you make it proportional to that with a different factor, you get a different \(\boldsymbol\lambda\), but the final \(\mathbf x\) projection still turns out to be the same.) The change in the constraint value due to this displacement is \(\Delta\mathbf c=\mathbf J\Delta\mathbf x=\mathbf J\mathbf M^{-1}\mathbf J^T\boldsymbol\lambda\).

Therefore, what the matrix \(\mathbf J\mathbf M^{-1}\mathbf J^T\) does is relate Lagrange multipliers \(\boldsymbol\lambda\) to the change in constraint value \(\Delta\mathbf c\) they produce. But we *know* how much we want to change the constraint value, \(\Delta\mathbf c=-\mathbf c\), and what we want to *find* is \(\boldsymbol\lambda\), so…

A nearly identical argument applies to the velocity projection.

Here’s a reference implementation of constrained dynamics.

Starter code for grid-based fluid simulation

This time I’ve added a little more direction in the TODO comments in the code. The starter code provides an implementation of regular and staggered grids, along with some interpolation routines for them, so take a look at how those work. It also includes a simple wrapper around Eigen’s sparse matrix functionality.

The simulator will closely follow Zhu and Bridson’s paper “Animating Sand as a Fluid” (2005). Here’s the order I recommend to do things in:

In

`fluid.cpp`

, first implement the`advection()`

routine to move each particle according to the velocity field on the grid. Use second-order Runge-Kutta, a.k.a. the explicit midpoint method, that is, \(\mathbf x \gets \mathbf x + \mathbf u(\mathbf x + \frac12 \mathbf u(x)\Delta t)\Delta t\). Now you should be able to drag the particles around a bit. (This is Zhu and Bridson sec. 4.2.5.)Implement

`particlesToGrid()`

and`gridToParticles()`

.In

`particlesToGrid()`

the goal is to set each grid velocity to the weighted average of nearby particle velocities. In practice, we will have each particle add its velocity to the nearby grid nodes using`addInterpolated()`

, then divide by the total weight in the end. You can keep track of the total weight by adding \((1,1)\) to`velWeight`

for each particle. (Zhu and Bridson sec. 4.2.2.)In

`gridToParticles()`

, the PIC update rule is \(\mathbf v \gets \mathbf u(\mathbf x)\); the FLIP update rule is \(\mathbf v \gets \mathbf v + (\mathbf u(\mathbf x) - \mathbf u_{\text{old}}(\mathbf x))\). Compute both, and set \(\mathbf v\) to 95% FLIP + 5% PIC. (Zhu and Bridson sec. 4.2.4.)

Now the particles should keep moving after you let go.

Implement

`applyBoundaryConditions()`

. In our current setup of a rectangular grid filled with fluid, this just means setting the normal velocities at the ends of the grid to zero. Remember that because we’re using a staggered grid for velocities,`vel.u`

is an \((m+1)\times n\) grid while`vel.v`

is \(m\times(n+1)\).Finally, implement

`pressureSolve()`

. This is the big one.First, we have to enumerate all the fluid cells, as those are the ones which contain pressure variables. Loop over all the grid cells and count the ones flagged as

`FLUID`

; store their indices in the corresponding location in the`indices`

grid. (Right now this step may seem kind of unnecessary as the whole grid contains nothing but fluid cells, but it will be needed if we later want to simulate a liquid with a free surface.)Second, build the linear system corresponding to the pressure equation \(\Delta t\nabla^2 p = \nabla\cdot\mathbf u\). Suppose \(N\) is the number of fluid cells you counted before. Create an \(N\times N\) sparse matrix

`A`

which will correspond to the operator \(\Delta t\nabla^2\), and an \(N\)-vector`b`

which will correspond to \(\nabla\cdot\mathbf u\). For each fluid cell, compute the corresponding entries of`b`

using the neighbouring velocities, and populate the corresponding row of`A`

by looking at the flags of the neighbouring cells.Third, solve the linear system to get the pressure values, and copy them to the cells of the

`pressure`

grid. At this point, you should be able to see the pressure field in the simulator, even though the flow itself is unaffected.Finally, update the flow velocities according to the formula \(\mathbf u \gets \mathbf u - \Delta t\nabla p\), using the difference of pressure in adjacent cells. Now you have an incompressible fluid simulator!

Section 4 of Bridson and Müller-Fischer’s course notes on “Fluid Simulation for Computer Animation” (2007) may be helpful for understanding the details of the pressure solve. Use the picture on the right to help you remember which pressure nodes and staggered velocity nodes are adjacent in our implementation.

So there’s a quick-and-dirty way to implement liquid simulation, and then there’s the better way which is how most people actually implement it. However, since we’re far enough along into the semester that you should be focusing on your final projects right now, we’ll just do the quick-and-dirty solution, and I’ll leave links about the better technique for future reference.

Change

`flagCells()`

so that only cells that contain fluid particles are marked as`FLUID`

.Change

`init()`

so that the fluid particles don’t fill the whole grid. For example, you could make`i`

and`j`

only go up to`m/2`

and`3*n/4`

for a “breaking dam” type scenario. Also, right now the`i1`

and`j1`

loops create \(2\times2\) particles per grid cell; however, in 2D it seems that using \(4\times4\) works better.Implement a

`GravityForce`

that applies a constant, uniform acceleration to all velocity values.Reduce the time step

`dt`

in`main.cpp`

to something smaller, like`0.01`

.

That should do it.

There are two major limitations to this approach:

We’re doing a binary classification of grid cells into totally liquid or totally air. This will produce stairstepping artifacts whenever the surface is sloped. What we should be doing is to perform a surface reconstruction of the liquid-air boundary and classify some cells as partially containing liquid, along with some notion of where the free surface passes through those cells (Zhu and Bridson 2005, sec. 5). Then we’d have to change the pressure solve to take that information into account (Bridson and Müller-Fischer 2007, sec. 4.5.1).

The velocity field is only defined inside the current liquid volume. In the advection step, particles at the leading edge of the liquid may sample the velocity outside the current volume, and will then get a velocity of zero. One solution is to simply extrapolate the velocity field into the air volume, so that sampling it outside will still give a decent estimate (Bridson and Müller-Fischer 2007, sec. 6.3).

It’s pretty useful to look at example code for this stuff. Robert Bridson has a simple 2D PIC/FLIP implementation on his website; scroll down and look for the link to “simple_flip2d.tar.gz”.