Octrees are fundamental in many computer graphics and simulation applications, but I found myself writing one for the first time last weekend, so I thought I’d share the basic ideas and show some code demonstrating how it could be used. The code for this post is available on Github for those eager to jump the gun, but below is an explanation of the semi-interesting parts of the code.

What is an Octree?

An octree is tree built over a region of space, allowing for fast queries for nearby objects. As the name implies, every interior node in the tree has 8 children. The 8 children are determined by the physical center point of their parent node. Pictured above is a simple octree. In 2D, there is a very similar data structure called the Quadtree, which is the same thing but with only 4 children per node since it is in 2D. (You can imagine analogous data structures for higher dimensional spaces with $2^d$ children per node.) A great visualization of an octree is available here.

Typically, we will store objects in an octree so that we can have fast queries for questions like “what points lie withing a certain bounding box”, or equivalently, which points are within a distance $r$ from some query point. There are other uses that we will explore in a later article on particle-based fluid dynamics.

The idea of storing objects is pretty simple in an octree. The actual data that you will eventually want to query will be stored in leaf nodes, and at least for the purpose of this article, interior nodes of the tree will not store data. Oftentimes we will just store one point in a leaf node, but it is easy to handle more than one point in a leaf. We should also notice that

How can we build an Octree?

Building the octree simply requires $n$ calls to $insert$, one for each point to be inserted into an empty tree. The logic for inserting into the tree is relatively straight forward. There are precisely three cases which must be considered for inserting points into an octree.

1. The node is an interior node to the tree (has 8 children).
   

   Since we never store data in an interior node of the tree in this article, we will find out which of the eight children the data point lies in and then make a recursive call to insert into that child.

2. The node is a leaf (no children/not split) and has no data assigned.
   

   This is the easiest! We’ve ended up in a small region of space with no data currently assigned and no children, so we will simply assign this data point to this leaf node and we’re done!

3. The node is a leaf (no children/not split), but it already has a point assigned.
   

   This is slightly more complicated. We are at a leaf but there’s something already here. Since we only store one point in a leaf, we will actually need to remember what was here already, split this node into eight children, and then re-insert the old point and our new point into the new children.  Note: it’s entirely possible that this will happen several times during insert if these two points are really close to each other. (On the order of the logarithm of the space separating them.)

This is hopefully pretty simple. I’ve tried to comment the code for insertion well so that this is as understandable as possible.

How can I determine which child of a node contains a point?

Since we use this several places, it’s important to know why/how this works, but it’s really simple. Remember that every node has a center point (which I called the “origin” in the code) and that all of the child nodes meet at this center point. There are eight children, which is exactly three bits of data. The children are assigned numbers in a simple way so that we can perform this operation quickly. If the point’s x-value is greater than the center point’s x-value, then the third bit of the children index will be set, otherwise it will be zero. Similarly if the y-value of the point is greater than the center point’s y-value, than the second bit is set, and similarly for the z-value. This basically means we do three comparisons and a few bitwise ORs in order to get the child index.

How can I quickly find all the points in a query box?

The main point of making this octree was so we could build the tree with tons of points in it before we start rendering/simulating something/etc and then quickly make many calls to “getPointsInsideBox” during execution. This is also pretty straight forward. The input to this algorithm is a axis-aligned bounding box defined by it’s minimum and maximum corner points and a list which we will insert our results onto. Like when we were building the tree, we can imagine this as a number of “cases” and just look at what we should do in each case. We begin by looking at the root:

1. If this node is a leaf node and we have data, see if the point is contained within the query box. If it is, then insert it into the results list.
2. If we’re at a leaf node and there’s nothing stored here, then we won’t insert a result for this node.
3. We must be at an interior node in this case. We will look at each of this nodes eight children. If the child’s bounding box lies outside the bounding box, then there is no need to check that child or any of it’s descendants in the tree. In any other case, we will recursively look into the child node (step 1).

Note: If the size of the querying box is large with respect to the whole space the tree encompasses, then you will visit many (or, in the worst case, all) nodes in the tree. When you are looking at this many nodes, it may actually be faster to use another method. In the case where the result is all the points, you can see how the octree actually performs worse than brute-force checking every node manually! But, octrees give huge gains when the query is a small region of space compared to the size of the tree’s space, as the code demonstrates.

Example Code

The example code at github is meant to educate rather than be extremely optimized. I hope there is a good amount of explanation between this post and the code comments. If not, please leave a comment here!

For those curious, the opening picture depicts an octree used for quickly simulating fluids using “vortons”, which will be the subject of the next article

The Basics

The wave equation is a classic equation in partial differential equations (PDEs). It has a number of uses in physics and is often used to explain energy transfer through space. It is an example of a time-dependent equation, which simply means that the “solution”, or state of the system, changes over time. The PDE explains how the solution takes place. In the wave equation that we’ll simulate, we will use the notion of the Laplacian, which is a linear operator. The practical implication of the “linearity” of this problem is that there are good analytical insights into the problem and how it behaves, fast ways to solve it, and many different types of ways to go about solving it.

In this article, I’ll assume the reader has a decent understanding of calculus and linear algebra, but I’ll also try to explain things in as intuitive a way as I can..

The classical wave equation has the form

$$\frac{\partial^2}{\partial t^2}u = c \nabla^2 u$$

where $u$ is the “solution” which possibly changes value as we move around in space, and as we vary time. This means that $u$ is more properly written as a function of space and time, like $u(x,t)$, but for brevity, we’ll drop the extra notation. In order to simulate waves, we will need a way to describe the Laplacian $\nabla^2$, and prescribe a way to advance our current solution at a point in time to the next point in time, so we can simulate and animate this process.

Waves on a Mesh

We can find descriptions of wave simulations all over the web, so why is this article any different? Most of the articles on the web assume that space is “flat”. In 2D, the term “flat” just refers to simulating the wave as something that moves on a flat plane, rather than on rolling hills for example. What this usually means programming wise, is that there is a regular (equally-spaced) grid of points and we calculate how the value of $u$ changes over time at each point in the grid in order to generate waves. While this is perfectly valid, it doesn’t allow us to do things like simulate waves on the surface of a sphere, which is obviously curved everywhere! We’ll see that there is actually a very nice method of simulating waves on a triangular mesh (surface) based on the idea of “flux”.

The Laplacian, Divergence Theorem on a Triangular Mesh, and Flux

The Laplacian is a fundamental linear operator in PDE Theory and physics. It has deep connections to the notion of curvature. It actually has a very simple meaning though: It is simply the sum of all the second (non-mixed) partial derivatives. A slightly more rich definition, which probably sounds more confusing but is actually the one we need in this article, is that it is the divergence of the gradient. One simple (“layperson”) way to think about the Laplacian is the following: If the value is constant around a point, then the Laplacian there will be zero. This is like a flat plane, no curvature means zero Laplacian. If the value is change, but the rate at which it’s changing is a constant, then this will also be zero. This is like a hill with no valleys or crests. However, if there is a crest/valley, then there will be a non-zero Laplacian value there. It’s also worth noting that this value is signed depending on if the “bump” is locally pushing outwards or inwards relative to neighboring values.

Anyway… We mentioned that the Laplacian has a definition as the divergence of the gradient:

$$ \nabla^2 u \equiv \nabla \cdot (\nabla u)$$

which means that our wave equation must also be equivalent to

$$ \frac{\partial^2}{\partial t^2} u = \nabla \cdot (\nabla u) $$

The divergence theorem tells us that if we take a region of space and integrate the divergence of something in that region, it will be equal to the net flow across the boundary of that region:

$$ \int_\Omega \nabla \cdot \mathbf{g} = \int_{\partial\Omega} \mathbf{n} \cdot \mathbf{g} $$

Another way to say this is that if we imagine the net effect of adding up all the “compression” and “expansion” occurring within a space, that will be equal to adding up how much stuff crosses the boundary of that space. Practically though, this will lend us good insight into a way to solve the problem.

In our problem, we will assume that we have a value of $u$ known on each triangle in a manifold mesh. If there’s a question about what manifold means here, you can either google, ask in a comment, or basically imagine almost any model you would make by subdividing a sphere in Blender/Maya/etc and pulling things around. One fun thing about this method is that the mesh should be topologically manifold, but the immersion doesn’t have to be. What that means is we can do strange things like make surfaces that are manifold but don’t have any way to exist in reality, like a Klein Bottle, and actually simulate waves on them.

Since we have a value of $u$ on each triangle, we will want to know how to change that value over time, but in order to do that, we will want to know how to compute the Laplacian at each triangle. We will utilize our last result from the Divergence Theorem, which told us that the integral of the Laplacian was equal to adding up the gradient across the boundary. Integrating the entire wave equation over a single face $\Omega_i$

$$ \frac{\partial^2}{\partial t^2} \int_{\Omega_i} u_i = \frac{\partial^2}{\partial t^2} |\Omega_i| u_i = \int_{\partial\Omega_i} \mathbf{n} \cdot \nabla \mathbf{u} $$

$$\frac{\partial^2}{\partial t^2} u_i = |\Omega_i|^{-1} \int_{\partial\Omega_i} \mathbf{n} \cdot \nabla \mathbf{u}$$

To illustrate the “normals” for a the boundary of a face:

So now we need to calculate the gradient of the solution $u$ at the “boundary” of the triangle. The boundary of the triangle is just formed by the edges!. One simple estimator for the gradient across an edge is to calculate the difference between the value at the triangle across the edge $u_j$ and our value $u_i$ and to divide by the shortest distance between them.

Calculating this distance seems simple enough, right? Just get the centers of the triangles $C_i$, $C_j$ and compute $|C_i – C_j|$, right? No. The “straight line” distance is not correct, since our notion of distance now depends on the path lying on the surface. That means that if we’re given two triangles that share an edge, we need a path that goes from one triangle center to the other and that path must stay on the triangles. It’s “obvious” that this path must cross the edge at some point along the path. On each triangle, a portion of the path must go from the center to this special point on the edge. Since the triangles are themselves flat though, we know that this portion of the path will be a straight line. This results in the following: The minimum distance between the two centers is made by drawing a straight line from one center to a special point on the edge, and then continuing onward from that point on the edge to the center of the second triangle!

So this is actually pretty simple, but we now need to know what we should choose for center points, and also what the point on the edge, $e^*$ is that makes the total path length $|C_i-e^*|+|C_j-e^*|$ as small as possible. The proof is simple, but if our triangles have a property known as being “well-centered”, then they will have a special center called the circumcenter which lies inside the triangle. If all of the triangles have this property, then we can use them as the center points $C_i$, $C_j$, and the edge point to choose for the middle of the path is simply the middle of the edge. This is pretty simple to show with basic trig/geometry. Later perhaps. Side note: Another way to derive this would be to fix one of the triangles, compute a rotation that rotates the neighbouring triangle into the same plane, compute the straight line between them, and then find where this intersects the edge joining the two.

This is in fact all we need to compute the integral along the boundary of a triangle! For each triangle, we will compute the Laplacian there as the sum of gradients across the edges of the triangle. A natural result of this is that if triangles along the boundary of the surface (if there is a boundary at all) will have less than three edges contributing to the value of the Laplacian.

So, procedurally, we will need to precompute a few things in order to calculate the Laplacian: the circumcenters, the area of triangles, edge midpoints, etc. which can all be calculated pretty easily. From then on, the Laplacian at a triangle is computed as the sum of gradients across the edges:

For each Triangle Ti:
  L[i] = 0 # Laplacian at Ti
  For each Triangle Tj neighboring Ti:
    Let eij be the edge shared by Ti and Tj, having vertices a and b
    Let Ci be the circumcenter of Ti, Cj the circumcenter of Tj
    Let Ai be the area of Ti

    # midpoint of eij
    E = .5 * (eij.a + eij.b) 

    # Compute the length of the path joining the centers
    D = |Ci - E| + |Cj + E| 

    # Compute the flux across this edge
    F = (uj - ui) * |eij| / D 
    L[i] += F / Ai

While this will work, it is a little slow… Neighboring triangles will compute the same value twice, just with opposite signs, so in practice they can be computed simultaneously, i.e. one per pair of neighbouring faces. Also, I found it much nicer to encode this into a matrix: Since $u_j$ and $u_i$ are the only unknowns at runtime and they only appear linearly in the expression for the Laplacian (this is expected since the Laplacian is linear!), we can write the Laplacian as a matrix $\mathbf{L}$, so the laplacian at each triangle is computed by just multiplying $\mathbf{Lu}$. The source code will make this more obvious! I have created a very small but working sparse matrix class for doing these calculations.

By the way, if face areas are neglected while you’re forming the Laplacian, the result may still “work”, but you’ll notice funky behaviour when a wave front crosses a triangle of a different size from it’s surrounding neighbours!

Obligatory Demo

A demo of this idea is viewable at GitHub! I highly recommend you run this code in Chrome for all of the typical reasons.

The source code is too!

If you have any questions or comments, please do write below!

In the first post introducing The Construct, I showed how the library could be used to create a simple fluid simulation with little effort. The visual effect wasn’t really that awesome for a couple of reasons. We’ll extend that same “project” in this post and show how we can implement a couple of ideas and papers from Computer Graphics literature with little additional effort.

Artificial Viscosity

One of the main problems with the Semi-Lagrangian advection (SLA) approach is that, being only an approximate solution, it introduces error into the system. When we advect quantities that lie in a grid around using SLA, quantities are interpolated from nearby points and this leads to a kind of diffusion of the velocity, density, or whatever happens to be advected. This turns out to be a difficult problem to solve in general and there isn’t really any easy and efficient silver bullet to fix it.When the velocity field is advected using SLA, it causes a loss of energy and vorticity (which is a measure of spinning and rotation in the fluid). Since this spinning motion is the lively thing that we would like to capture in the computer graphics setting, we would like a way to mitigate this loss of vorticity.

One technique, called Vorticity Confinement, seeks to mitigate this effect by essentially measuring where vorticity is being lost, and to inject rotational motion back into the flow. While it’s not a perfect solution, it is pretty easy to implement and does give a more lively flow. Mathematically, the procedure of computing Vorticity Confinement is also pretty straightforward.

We first calculate the vorticity by taking the curl of the velocity

$\B \omega = \nabla \times \B u$

In The Construct, this is simply computed as

VectorField omega = curl(velocity);

Once we have the curl, we compute the the gradient of the magnitude of this vorticity, which will be a vector field pointing towards locations of high vorticity. Mathematically we have

$\B \eta = \nabla |\B \omega| $

$\B N = \B \eta / |\B \eta|$

which is handily translated to

omega = writeToGrid(omega, constant(Vec3(0,0,0)), domain);
VectorField eta = grad(length(omega));
VectorField N = eta / (constant(.00001) + length(eta));

Just a couple of notes here. First, we had to write omega to a grid before we created eta, for the simple reason that $\B \eta$ involves a gradient of the curl, which would require analytically taking two derivatives of the velocity, which isn’t possible in The Construct. By writing it to a grid, we can take another derivative without worry. Secondly, when we compute N, we had to add a small epsilon to the denomintator since it is possible (and even likely) that some places will have zero velocity/curl which would create a divide by zero. This just ensures that we do not divide by zero when we normalize this vector field.

Finally, the actual Vorticity Confinment force pushes perpendicular to this gradient like a paddle wheel churning up turbulence:

$\B f_{VC} = \epsilon h (\B N \times \B \omega)$

where $\epsilon$ is a user parameter adjusting the strength of the confinement forcing, and $h$ is the side length of a voxel in the grid. This means that the entire Vorticity Confinement force comes down to just a few more extra lines of code, which for the most part, mirror the mathematics almost exactly as described in the paper:

VectorField VorticityConfinement(
 VectorField velocity, Domain domain, float epsilon) { 
   VectorField omega = curl(velocity);
   omega = writeToGrid(omega, constant(Vec3(0,0,0)), domain);
   VectorField eta = grad(length(omega));
   VectorField N = eta / (constant(.00001f) + length(eta));
   return constant(epsilon * domain.H[0])) * cross(N,omega);
}

Applying this force to the original simulation gives us these results:

The result is pretty subtle, but if you looks closely, you will notice that there are finer features, and more importantly that the motion of the fluid does not dampen as quickly as the vanilla simulation.

Rayleigh-Taylor Instability

One interesting scenario that we can compute demonstrates a phenomenon called the Rayleigh-Taylor Instability. This instability refers to a shift to chaotic motion after a symmetry is broken across an interface of two varying densities. Essentially you can imagine a “heavy” material resting in a perfectly balanced way on top of a less dense material. So long as there is no net force across the boundary, the system will remain in a quasi-stable state. The moment a slight disturbance destroys this perfect symmetry, there is a net force which causes the motion to become chaotic.

We can mimic this scenario by create two densities, a red one on top that is heavy, and a blue one on the bottom which is lighter.

auto red = mask(dot(identity(), constant(Vec3(0,1,0))));
auto blue = mask(dot(identity(), constant(Vec3(0,-1,0))));

We only need to augment the force now so that red pushes downward, and blue pushes upward, and have each of these advect along the velocity field independently, since we have two different densities now.

for(...) {
  red = advect(red, velocity, delta_t);
  red = writeToGrid(red, constant(0.f), domain);
  blue = advect(blue, velocity, delta_t);
  blue = writeToGrid(blue, constant(0.f), domain);

  VectorField force = (blue - red) * constant(Vec3(0,1,0));
  ...
}

One small extra thing we can do is add a small bit of asymmetry to one of the densities so that the turbulent transition will occur more quickly. After red is initialized, we can just add a small bit of extra density to break things up:

VectorField x = identity(), N = constant(Vec3(0,0,1));
ScalarField extra = mask(constant(.05f) - length(x-N*dot(x,N)));
red = red + extra;

Putting these things together produces the following result:

Gridless Advection

This next example is rather curious in that we will only need to remove code in order to see its effect. Strictly speaking, we do not need to write the density into a grid at each time step. We have been writing the density into a grid at each time step because this means that the subsequent time step will only need to read from a grid to get a value for the density. If we do not write any of the density to a grid at each time step, we will construct a larger and larger expression. Writing density to a grid makes things fast, but we do pay a price. Because later reads from that grid require interpolation, diffusion takes place, which makes the density “smear” out as it flows. If we avoid writing density to a grid, then this diffusion is removed!

Removing those lines is easy, but be prepared for long render times. Two examples of the results of this are below, rendered from an angle with a version of my own renderer (approximately 2 hours to render the higher-resolution video):

SELMA and the Characteristic Map

Gridless Advection was really quite nice visually, but takes forever to compute, because each evaluation at render-time requires following back a number of times along the velocity field. Fine details that don’t diffuse out are one feature of the Gridless Advection approach that is useful in certain circumstances and isn’t too hard to capture approximately by making use of a construction called the Characteristic Map.

The mapping function phrases the movement of density in a special way. Instead of moving density around on a grid, we have a “mapping function” which changes at each time step that maps where material “came from” after a period of flow. The paper explains the details behind this map, but implementation follows directly from the fact that a map at time $n$ is based on the previous map:

$\B X_{n+1}(\B x) = \B X_n(\B x – \B u \Delta t)$

with the initial condition $\B X_0(\B x) = \B x$. With this map, the density at time $t$ is given by a combination of the map at time $t$ and the original density only:

$\rho_n(\B x) = \rho_0(\B X_n(\B x))$

In the Construct,  this looks like

ScalarField initial_density = ...;
VectorField X = identity();
...
for(...) {
  // Update the mapping function by advecting coordinates.
  // Please see the paper for details
  X = advect(X, velocity, dt);
  X = writeToGrid(X, identity(), dt);

  ...

  // Now the density is given by mapping the density through
  // the Characteristic Map
  auto render_density = warp(initial_density, X);
  render(render_density, ...);
}

A comparison of these techniques is shown below:

On the left is the vanilla traditional rendering of density, in the middle is the CM-based rendering of the same density, and on the right is the CM-based rendering of a different density (at no additional cost). As you can see, the CM-based approaches have much less diffuse appearance to the density.

 More to Come

More examples are coming! If you have questions or comments, please leave them here.

In this post, I am introducing a framework that I wrote as part of my MSc degree. If you have any questions or comments, or would like to help develop this project further, please comment or contact me!

What the heck is “The Construct”?

The Construct is a high-level C++-based library for doing 3D “volumetric” computations, with Python bindings in the works. These computations allow us to do things like, say, smoke simulation in a very small amount of code. This is possible because concepts in this library are expressed in a very natural and mathematical language.

At it’s base, The Construct allows us to express mathematical ideas that take place in 3D. For instance, consider the following function:

$f(\B x) = \frac 1 {1 + \B x \cdot \B x}$

There are a ton of ways we could implement something like this in any number of languages; I mean, come on, right? It’s just one dot product and then a divide. Simple stuff.

Let’s say though that I instead asked you for all the partial derivatives of a more complex expression:

$\nabla g(\B x) = \nabla \left( \frac{ \B u(\B x) \cdot \B v(\B x)}{ 1 + \B x \cdot (\B u(\B x) – \B v(\B x)) } \right) = ?$

This is a bit more involved! Evaluating just $g$ takes only a little bit more effort to work to code,  but evaluating its gradient would require us to first work out the mathematics behind it, using all the appropriate quotient and several applications of product rules from multi-variable calculus! If you’re trying to test out new ideas, you really don’t want to have to waste time on these kinds of things when you can avoid it.

It turns out that we can compute $\nabla g(\B x)$ very easily in The Construct:

auto g = dot(u,v) / (1.f + dot(x, u-v))
auto gradg = grad(g);

After we have created an object like this, we can use the syntax

std::cout << gradg.eval(Vec3(1.2, 3.7, -2.9)) << std::endl;

to evaluate this expression exactly at some point in space. What’s more, if I change g to be a multiplication instead of a division, I don’t have to go through the process of deriving the gradient of the new expression again!

Basics

There are 3 major types of objects in The Construct: ScalarFields, VectorFields, and MatrixFields. These correspond to functions/expressions which can be evaluated at a 3D position and return, respectively, a real number, a 3D vector, or a 3×3 matrix of values.

There are a large number of supported operations on these fields: addition subtraction, scalar multiplication, dot product, cross product, outer product, matrix-vector products, spatial gradients, line integrals, and many more. What this large library of operations allows us to do is build up tools for creating volumetric content (which can be later rendered or saved to disk).

One of the cool things about how this system is organized is that it doesn’t presuppose any notion of a grid. Many mathematical tools allow operations on grids, but assume for rendering or other purposes that content is on a grid. While objects in this library can be made into grids (and as we’ll see later, may require being “stamped” into a grid), there are essentially no requirements on data being gridded. This is why calculations are exactly evaluated. One of the only major restrictions (which can partially be circumvented), is that if we can not use grids to represent data, then derivatives can only be taken once. For example, it is impossible to compute the Hessian of a scalar field exactly without using grids in this framework.

Let’s see some examples!

So this is all “neat”, but how is it actually useful? In the development of this library, virtually all of these tools were designed with a degree of flexibility and also speed in mind, as I used them for experimenting new techniques for simulating and using fluid dynamics of smoke for computer graphics purposes. In the next section, we’ll show how to actually write a simple fluid simulation that you can mess around with and easily extend in approximately no time flat. We’ve already demonstrated how lots of algebraic operations can be done. We’ll go ahead and show a few more examples though, some of which we’ll need in the next section.

constant()

The simplest fields in the library are constant ones: Independent of where they’re evaluated, they return a single scalar, vector, or matrix. They are often not explicitly needed since C++ can provide an implicit construction of these constant fields from the types float,Vec3,Mat3 for scalars, vectors, and matrices respectively.

ScalarField C1 = 4.5f;
VectorField C2 = Vec3(1,2,3);
VectorField C3 = C2 / constant(1.3f);
C3 = C2 / 1.3f; // Equivalent to previous line

identity()

Another of the simplest built-in functions is identity(), which returns a vector field that performs a really simple task: when evaluated, it returns a vector equal to the location passed in! This is useful when we need to build expressions involving the location the are evaluated at in the function itself.

warp()

Another function that is fairly simple but incredibly important is the “warp” function. This function takes some field $f(\B x)$ of any type, and a vector field $\B g(\B x)$ and returns a new field which returns a value from $f$, but evaluated at $\B g(\B x)$, which is exactly function composition. So, mathematically, if h=warp(f,g), h will be evaluated like $h(\B x) = f(\B g(\B x))$.

mask()

mask(…) takes one argument, a scalar field, and returns a new scalar field which takes on the value 1 when the original field is positive, and 0 elsewhere. This is useful in “selecting” areas of space to fill with a value or to apply an effect to.

// Create a field which is the signed distance to 
// the surface of the unit sphere
ScalarField sphere = 1.f - length(identity());

// Evaluates to one inside the sphere, zero elsewhere
ScalarField unit_ball = mask(sphere);

Domain() and writeToGrid()

Fields can exist as grids in The Construct, but their existence is transparent: once they have been created, they act just like regular fields. Grids are defined by the “lower left” and “upper right” corners of the axis-aligned bounding box enclosing the grid, as well as the resolution:

int res[3] = {64, 128, 64};
Vec3 bmin(-1,-2,-1), bmax(1,2,1);
Domain domain(res[0], res[1], res[2], bmin, bmax);

We can later use this Domain when we need to write a field to a grid. The writeToGrid(…) takes three arguments: First, the field which we would like to write to a grid, second a field which should be evaluated in the event that something tries to access data outside the bounding box of this grid, and last, the actual domain which defines the grid. As an example:

Domain domain(64, 64, 64, Vec3(-1,-1,-1), Vec3(1,1,1));
VectorField x = identity();
ScalarField f = 1.f / (1.f + dot(x,x));

// Write f to a grid, return 0.f if we call eval outside the grid...
ScalarField g = writeToGrid(f, constant(0.f), domain);

// Interact with gridded fields just like regular ones!
ScalarField h = (f + g) * .5f;

Writing a Simple Fluid Simulator in One Minute

Fluid simulation is typically a complex process to model and simulate. In engineering and mathematics circles, it has spurred research and software development for the last 50 years, with whole careers made in forming better and more accurate tools for understanding flow. In computer graphics however, we typically care only that things look nice or realistic, and so a number of simplifications are often made at the expense of accuracy. There is a kind of popular paper by Jos Stam, which outlines a method for simulating fluids which has been adopted … extensively, over the last decade for simulating smoke and water. This technique is fairly easy to follow and understand, but isn’t really in the scope of this article. If what I’m about to show you is confusing, then suffice it to say that most implementations written in “typical” programming languages like C, Java, tend to be much longer (JavaScript!), harder to debug and maintain, and are more difficult to add features to.

The Stable Fluids approach really has two parts: advection and pressure projection. The advection procedure basically moves material and velocity around in space, while the pressure projection step insures that the velocity field is divergence free, leading to the nice vortex and swirling motion. Moving material around typically is done  by using Semi-Lagrangian advection. Mathematically, this process means we should sample a field along lines dictated by the velocity field backwards in time in order to see where material flows to. It looks like this:

$Q(\B x, t + \Delta t) = Q(\B x – \Delta t \B u(\B x, t), t)$

So, if you repeat this process many times, you can get your various quantities like density of smoke (or, $Q$ in this case) to flow around in the velocity field. We can easily package this up using the library with the following code:

template<typename T>
Field<T> advect(Field<T> Q, VectorField velocity, float deltat) {
   return warp(Q, identity() - deltat * velocity);
}

What we’ve done here is make a one line function which corresponds exactly to one equation! It takes in some field Q (which could be a scalar,vector, or matrix field, hence the templating), a velocity field, and a time step. The function returns a new field of the same type that is essentially evaluated in a new place. The only two mysteries for new people here are what the “warp” and “identity” functions mean, but these were explained in the last section.

Typically, this advect() would have taken around a dozen lines of code in other languages and would have to deal with things like interpolation on grids and other minutia not really pertinent to what we’re trying to do here. Instead we have a short, easy to understand single line which we could easily play around with and get new behavior. The tiny details have been wrapped up for us within The Construct.

The other part of fluid simulation, the pressure projection which ensures that the velocity is divergence free, is (in a bit of a cheating way) precisely one line of code. Pressure projection involves solving a linear algebraic problem in many unknowns, dependent on the simulation mesh and divergence of the velocity field.

velocity = divFree(velocity, constant(0.f), domain);

this divFree(…) function takes in the original velocity field, a second argument which is used in determining boundaries inside the flow, and a domain definition which is used when the velocity is written to a grid (which is required, since this process can not be done without being discretized). This expensive step, which is actually where most of the processing is taking place, has a couple of assumptions on boundary conditions etc. that have been assumed for the time being. In a later release, more of this will be abstracted. For the time being though, this is a very convenient one-line solution to the divergence-free projection problem.

With all of these things together, the actual fluid simulation code is very, very small. We first define a simulation domain (along with some resolution, here it is 128^3 voxels), a density which is initially a solid sphere of material, a velocity (which is initially at rest), and a time step size. We then go through a hundred time steps, advancing the smoke/density by letting it flow along the velocity at that point in time, and then updating the velocity in “Stable Fluids” way. We then send this off to a routine which will actually render a picture for us for each time step so we can see what’s going on. The force term that actually makes things move around here is termed a “buoyant” force: it has an upward direction and its magnitude is proportional to the amount of density/smoke at a spot. This causes the smoke to travel upward and collide against the edge of our domain, where it curls around.

// Create a cube-shaped domain to do simulation in
Domain domain(128, 128, 128, Vec3(-1,-1,-1), Vec3(1,1,1));

// Create a source to pump in density
auto density = mask(.8 - length(identity()));

// Initially, everything is at rest: zero velocity
auto velocity = constant(Vec3(0,0,0));
ScalarField delta_t = .1f;

for(int T=0; T<100; ++T) {
  // Move smoke around in the flow
  density = advect(density, velocity, delta_t);
  density = writeToGrid(density, constant(0.f), domain);

  // Advect velocity 
  velocity = advect(velocity, velocity, delta_t);

  // Apply force upward, proportional to the smoke
  velocity = velocity + delta_t * Vec3(0,1,0) * density;

  // Enforce divergence-free condition
  velocity = divFree(velocity, constant(0.f), domain);

  render(density);
}

This code is essentially duplicated in the GitHub project as an example (along with the ultra-basic volume rendering and image output code) so that you can run it yourself. The actual result looks like this:

More to come

I hope you have found this cool. I’ve spoken with a number of people who have always been interested in toying with fluid simulation stuff but not been sure on how to code the under-the-hood parts. I think this kind of software can make this fun and useful for exploration. This is really only one example of the most simple kind of fluid simulator that you can build with these tools. I will be making a follow up post (edit: complete!) which shows how to implement several related papers which have come out in the last few years using these same tools (typically in just a few extra lines!).

Also, before I go, I should remind you that there is quite a bit more to be added to this project (including a number of utilities from my thesis work that have not been made public yet). I hope as I add new functionality, more and more people can use these tools.