# The Construct

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(…) 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 = writeToGrid(density, constant(0.f), domain);

// 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.

This entry was posted in Uncategorized. Bookmark the permalink.

### 13 Responses to The Construct

1. noiv says:

Is a shader language port also planned?

• brandon says:

Something like this might be possible in GPGPU languages, but I do not think shading languages currently have the flexibility to maintain the node-based representation that The Construct uses to encode mathematical expressions. Somewhat tangentially, I do have some automatic OpenCL code generation stuff working in a limited context, which I will be showing and explaining in a subsequent post.

• burito says:

Nvidia has their NV_GL_shader_program_store extension for GLSL which allows just that. I’ve been doing Sparse OctTree’s with it.

2. val says:

Cool stuff. Have you heard about theano ?

It is a Python library that allows you to define, optimize, and evaluate mathematical expressions involving multi-dimensional arrays efficiently, featuring – transparent use of a GPU, efficient symbolic differentiation and dynamic C code generation among others.

• brandon says:

That sounds pretty rockin’ and totally related to this work. I’ll have to check it out!

3. jelle says:

cool project Brandon. looking at the src on github I couldn’t figure out how the python bindings are generated / compiled, perhaps you could elaborate? finally, I can imagine VTK + Construct to be a very potent combination

• brandon says:

While I have created these Python bindings in a non-public version of this library, I am still busy trying to get them into the public release. They’ll be in the repository “soon”. The previous method used SWIG to generate bindings, but I may use a different to generate the bindings for this public release.

Your note on using VTK with Construct could be on the right track. I currently have my own volume rendering software that I have used, but plugging this into a more general/public rendering solution might be an excellent idea. I also toyed with the idea of adding this system into a tree-based diagram in Blender, so people could create volume content and do volume simulations from within Blender, but it appears that the GUI elements that I would need to use are hard-coded for specific behaviors and there is virtually no documentation.

4. Hi Brandon,

Also, are you familiar with the FEniCS project? At a glance, the main difference between the Construct and the tools provided by the FEniCS project appears to be that the FEniCS tools are designed for the implementation of finite element methods, whereas the Construct allows the implementation of finite difference methods – is this the gist of it? Are there other simulation types that can be implemented using the Construct?

There seems to be a great similarity between the language used in the Construct and the Unified Form Language used in the FEniCS tools (e.g. see http://www.fenicsproject.org/documentation/tutorial/fundamentals.html for a couple of examples) – I guess this is a consequence of them both being used to describe vector calculus expressions – I’d be interested to hear your further thoughts!

• brandon says:

I will be posting my thesis in a more public place soon, but for now, I can send it to you if you can send me an email.

I have only looked shortly at FEniCS before, but it looks like a very related concept. The main difference as far as I can tell is that that system is FEM-centric and bent on producing very accurate results that scientists could use in modeling, whereas Construct is designed with speed/flexibility in mind, rather than accuracy. Most expressions/fields that can be formed in the Construct have no notion of a grid until explicitly “write” them into a grid (if you even need to). As an example, in the next blog post, I’m showing something called gridless advection, when the smoke density is _never_ written to a grid, and so numerical diffusion of the density is gone (at the cost of much-increased computation time).

Other types of simulations are possible, but it will be a few blog posts before that is presented here.

The language that I see at the link you gave is definitely akin to my stuff, but while the FEniCS project does use a very similar language and does solve vector calculus problems (as I attempt to do), that system appears to be centered on solving the weak form of problems (as is typical for FEM). While FEniCS is almost undoubtedly more powerful and general than my software at the moment, I do like the fact that I do not have to derive the weak form of problems in order to get a solution =]

Thanks for your questions and thoughts!

• eric_t says:

Another similar code to look at is OpenFOAM:
http://www.openfoam.com/

It is similar to FEniCS, but uses the finite volume method instead, so you don’t use the weak form of the equations. OpenFOAM has quickly become very popular in the scientific community due to the ease of adding/modifying solvers.

5. Thrax says:

Thanks for posting this! I’m learning a lot from reading it, and the corresponding discussion on reddit as well.. Cheers!

6. maxim says:

Hello Brandon. Thanks for posting.
What about code for volume rendering?
I was wondering do you plan to publish the code to that used in this video: Renderer Test – High Resolution Smoke Sim or 4-Way Comparison between Density/SELMA Rendering and BFECC/SLA Advection.

• brandon says:

I’m so glad you asked. These are exactly the same kinds of tools that I used to create those animations. I am debating currently whether or not there is a rendering system which fits with the Construct nicely. If I can’t find one that works well, I’ll be putting out some simple rendering code — probably an offshoot of my volume rendering tutorial code which is already at GitHub.

In the next post, I’ll show how to do SELMA/Density rendering and how it captures sharp details. I need to add one more operator to this public code before the stable MacCormack advection can be implemented.

Thanks for the question!