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.