INSTRUCTOR: Adam Bargteil (Office hours: By appointment, NSH 4229)
WEB PAGE: http://graphics.cs.cmu.edu/courses/15-869/
WIKI: http://graphics.cs.cmu.edu/courses/15-869/wiki

For this assignment you will implement a basic fluid simulator. I recommend using the Stable Fluids approach with a staggered grid and buoyancy forces and vorticity confinement (see Visual Simulation of Smoke).

While I expect this assignment to be easier overall than the previous one, the data structures are more complex. Spend some time thinking about how to implement the static grid and how you will encode your boundary conditions. If you have an nxnxn grid of cells then you will have three arrays of velocity components u[n+1][n][n], v[n][n+1][n], and w[n][n][n+1]. The pressures, smoke density, etc will be stored at the cell centers, e.g. p[n][n][n]. Then the six faces around cell i,j,k are u[i][j][k], u[i+1][j][k], v[i][j][k], v[i][j+1][k], w[i][j][k], and w[i][j][k+1].

You will want to write your interpolation function before you begin writing the simulator. For each velocity component, you will in general interpolate the values stored at 8 faces. This can be a bit tricky. To test it I recommend setting the velocities to the half-index notation of position. That is u[i][j][k] = i*dx and v[i][j][k] = j*dx, where dx is the width of a grid cell. This assumes the origin of your coordinate system is at the lower corner (not center) of the 0,0,0 cell. If you interpolate this velocity field at any position x (in this coordinate system), you should get x as the answer.

Once you have interpolation working everything else is relatively easy. You'll also need a traceBack(x,dt) function which traces a point backwards through the velocity field for a timestep dt. Now you have everything you need for advection. If you decide to skip viscosity (a fine decision) all thats left is to compute the divergence, solve the poisson equation and subtract the gradient of the pressure field.

To compute the divergence, for each cell you subtract the velocity of the left (front, down) face from the right (back, up) face. I.e. div(i,j,k) = (u[i+1][j][k] - u[i][j][k] + v[i][j+1][k] - v[i][j][k] + w[i][j][k+1] - w[i][j][k]) / dx.

You can build the laplacian matrix explicitly if you like. Or you can just code the matrix vector product. For each cell, laplacian(p) = (p[i-1][j][k] + p[i+1][j][k] + p[i][j-1][k] + p[i][j+1][k] + p[i][j][k-1] + p[i][j][k+1] - 6 * p[i][j][k]) / dx^2. There could be another constant in there... I forget. When you are near the boundary, this is a little trickier. The key idea is to realize that the gradient across the boundary should be zero (so that we don't change the velocity on the boundary face). So, if the face between cell i-1,j,k and i,j,k is constrained, then we remove the p[i-1][j][k] term from the above and change the 6 to a 5. This essentially copies the pressure from i,j,k to i-1,j,k. To solve the poison equation, you need a linear system solver. This one looks reasonably easy to use: LinearSolver, but I can make no warranty and you should feel free to use whatever you like. Once you've solved the poisson equation you'll need to subtract the gradient of the pressure field from the velocity field. To do this simply do u[i][j][k] -= (p[i][j][k] - p[i-1][j][k]) / dx and the analogs in the other 2 dimensions.

Rendering