Box2D Fluid (Part 1)

Porting the JBox2D Liquid Demo to C# and XNA.

Source: FluidPort_1.zip
Github: klutch/Box2DFluid

This is the first in a series of posts that will explain how to add a fluid simulation to your Box2D game. I will be using C#, XNA and the Farseer Physics Engine.

Disclaimer: The main focus of this series is to show how to implement a fluid simulation that performs well and looks “good enough”. This isn’t meant to be a rigorous explanation of fluid dynamics, physics, or parallelization. I barely understand a lot of this stuff, and I’m sure there are better sources of information available on the web. What I’ll be presenting is the best I could come up with after a month or so of experimenting with fluids in Box2D. I make no promises that this is the best approach. It’s just the approach that has worked the best for me.

Having said that, please leave a comment if you have any questions, find any errors, or can think of a better way of doing something.

Introduction

During the early stages of development for Loder’s Fall, I decided that I wanted to add some fluid to the game. One of the most impressive solutions I found was JBox2D’s Liquid Demo. I like the JBox2D demo because it does more than just toss a ton of small Box2D circles on top of each other and call it “water”. I’m not educated enough to tell you exactly what type of simulation the demo uses, but it’s obvious that it handles particle forces outside of the Box2D engine in a way that gives the simulation a very nice, fluid appearance.

The basic idea of how the fluid simulation works is for the particles to adjust their positions and velocities based on the relative position between themselves and their neighboring particles. The particles have an “ideal” distance that they like to maintain. If particles are too close, they are pushed away from each other. To prevent the simulation from checking the distance between every single particle, a spatial partitioning scheme is used to simplify the search for neighboring particles. Basically, space is broken up into a grid of cells that contains particles. Any given particle’s neighbors can then be found by getting the particles from neighboring grid cells.

The main drawback of the JBox2D demo (for my purposes, at least) is its use of Box2D bodies as particles. While using bodies simplifies collision handling by letting Box2D resolve all of the collisions, it causes performance problems when trying to use too many particles. As far as I can tell, there is a lot of overhead involved with contact generation in Box2D, and because of that Box2D is not very well suited for particle simulations. Unfortunately, handling collisions outside of Box2D is hard. However, handling collisions externally seems to be the best way of running the simulation if performance is a concern.

Another drawback of the JBox2D demo is the grid it uses for spatial partitioning. The grid is limited in size, which means particles will only act as liquid when they’re inside the region of space that the grid occupies. Loder’s Fall is a kind of side-scroller/platformer game, so it was necessary for me to replace the demo’s fixed grid with one that existed anywhere fluid did.

Enough of that, let’s get on to coding.

Porting the JBox2D Demo

The source for the demo is available at the JBox2D repository.
I’ll be using Microsoft Visual Studio 2010 to compile the port.

The FluidSimulation Class

To start I’m going to create a new XNA 4.0 ‘Windows Game’ project in Visual Studio called FluidPort, and create a new class called FluidSimulation:

public class FluidSimulation
{
    private SpriteBatch _spriteBatch;
    private float _scale = 32f;

    public FluidSimulation(SpriteBatch spriteBatch)
    {
        _spriteBatch = spriteBatch;
    }
}

I’ll create an instance of FluidSimulation in Game1 ’s Initialize() method.

protected override void Initialize()
{
    base.Initialize();

    _fluidSimulation = new FluidSimulation(_spriteBatch);
}

Since we’re passing a SpriteBatch instance into the FluidSimulation, we have to call base.Initialize() first, otherwise the SpriteBatch instance is null.

The Particle Class

The demo stores the liquid as an array of Box2D bodies. I want to run the whole simulation outside of Box2D, so I’m going to go ahead and create a Particle class too.

public class Particle
{
    public Vector2 position;
    public Vector2 velocity;
    public bool alive;
    public int index;

    public Particle(Vector2 position, Vector2 velocity, bool alive)
    {
        this.position = position;
        this.velocity = velocity;
        this.alive = alive;
    }
}

Now lets start the actual porting. I’m going to handle particle creation a little bit differently than the demo, so we don’t need the boxWidth and boxHeight variables, which were used to define the initial position of particles in the demo. I’m also going to be using a different spatial partitioning scheme, so we don’t need the fluidMin and fluidMax variables. Right now, all we really need is

Some of those can be constants.

public const int MAX_PARTICLES = 1000;
public const float RADIUS = 0.6f;
public const float VISCOSITY = 0.004f;
public const float DT = 1f / 60f;
private int _numActiveParticles = 0;
private Particle[] _liquid;
private List<int> _activeParticles;

The _liquid array needs to be initialized and populated. This isn’t an array of active particles, it’s just an array of pre-initialized particles that can be activated/deactivated and used quickly. We’ll know which particles are active, because their index will be stored in the _activeParticle list.

The _liquid array and _activeParticles list need to be initialized, so add the following to the FluidSimulation constructor:

_activeParticles = new List<int>(MAX_PARTICLES);
_liquid = new Particle[MAX_PARTICLES];
for (int i = 0; i < MAX_PARTICLES; i++)
{
    _liquid[i] = new Particle(Vector2.Zero, Vector2.Zero, false);
    _liquid[i].index = i;
}

Since we’re skipping most everything that has to do with hashing (because we’re going to implement it differently), we can go straight to porting the applyLiquidConstraint() method.

The applyLiquidConstraints() Method

The comment at the top of applyLiquidConstraints states:

/*
 * Unfortunately, this simulation method is not actually scale
 * invariant, and it breaks down for rad < ~3 or so.  So we need
 * to scale everything to an ideal rad and then scale it back after.
 */

idealRad is then defined and assigned a value of 50f. A multiplier is also defined, and given a value of idealRad / rad . Since the value of these variables is never going to change, I’m going to move them out of the applyLiquidConstraint method, and make them constants in FluidSimulation:

public const float IDEAL_RADIUS = 50f;
public const float MULTIPLIER = IDEAL_RADIUS / RADIUS;

Next, the java code defines a couple of arrays for holding the change in position:

float[] xchange = new float[liquid.length];
float[] ychange = new float[liquid.length];

Instead of creating new arrays every time the applyLiquidConstraints method is called, I’m going to define a member variable that can be reused. I’m also going to make it an array of Vector2 instead of two arrays of floats. Add the following member variable to FluidSimulation:

private Vector2[] _delta;

and initialize it in the FluidSimulation constructor:

_delta = new Vector2[MAX_PARTICLES];

Next, the java code scales the particle’s position and velocity by the multiplier:

float[] xs = new float[liquid.length];
float[] ys = new float[liquid.length];
float[] vxs = new float[liquid.length];
float[] vys = new float[liquid.length];
for (int i=0; i<liquid.length; ++i) {
    xs[i] = multiplier*liquid[i].m_sweep.c.x;
    ys[i] = multiplier*liquid[i].m_sweep.c.y;
    vxs[i] = multiplier*liquid[i].m_linearVelocity.x;
    vys[i] = multiplier*liquid[i].m_linearVelocity.y;
}

Our code’s going to look a little different since we’re not using Box2D to get the position and velocity. Also, instead of creating a bunch of arrays every call, lets add some more member variables. Add the following to FluidSimulation:

private Vector2[] _scaledPositions;
private Vector2[] _scaledVelocities;

Initialize them in the FluidSimulation constructor:

_scaledPositions = new Vector2[MAX_PARTICLES];
_scaledVelocities = new Vector2[MAX_PARTICLES];

In applyLiquidConstraints(), we need to loop through the active particles and calculate the particle’s scaled position and velocity. We should also reset the values of _delta here to zero.

private void applyLiquidConstraints()
{
    // Prepare simulation
    for (int i = 0; i < _numActiveParticles; i++)
    {
        int index = _activeParticles[i];
        Particle particle = _liquid[index];

        // Scale positions and velocities
        _scaledPositions[index] = particle.position * MULTIPLIER;
        _scaledVelocities[index] = particle.velocity * MULTIPLIER;

        // Reset deltas
        _delta[index] = Vector2.Zero;
    }
}

Now for the important stuff. We need to loop through every particle and:

  1. Find the neighboring particles
  2. Calculate pressure
  3. Apply forces

We haven’t implemented the dynamic grid that will be used for neighbor searches, so we’ll just skip over that for now. However, we do need to give Particle a few new member variables in preparation for the neighbor search code. I will be storing the neighbors’ indices in an array in Particle. We’re going to be calculating distances between all of the neighbors, so we can store those in Particle as well. Give Particle the following member variables:

public float[] distances;
public int[] neighbors;
public int neighborCount;

Now we’re ready to port the code that calculates a particle’s pressure. This is the code we’re going to focus on:

// Particle pressure calculated by particle proximity
// Pressures = 0 iff all particles within range are idealRad distance away
float[] vlen = new float[neighbors.size()];
float p = 0.0f;
float pnear = 0.0f;
for(int a = 0; a < neighbors.size(); a++) {
    Integer n = (Integer)neighbors.get(a);
    int j = n.intValue();
    float vx = xs[j]-xs[i];//liquid[j].m_sweep.c.x - liquid[i].m_sweep.c.x;
    float vy = ys[j]-ys[i];//liquid[j].m_sweep.c.y - liquid[i].m_sweep.c.y;
    
    //early exit check
    if(vx > -idealRad && vx < idealRad && vy > -idealRad && vy < idealRad) {
        float vlensqr = (vx * vx + vy * vy);
        //within idealRad check
        if(vlensqr < idealRad*idealRad) {
            vlen[a] = (float)Math.sqrt(vlensqr);
            if (vlen[a] < Settings.EPSILON) vlen[a] = idealRad-.01f;
            float oneminusq = 1.0f-(vlen[a] / idealRad);
            p = (p + oneminusq*oneminusq);
            pnear = (pnear + oneminusq*oneminusq*oneminusq);
        } else {
            vlen[a] = Float.MAX_VALUE;
        }
    }
}

A couple notes:

Add the following to FluidSimulation:

public const float IDEAL_RADIUS_SQ = IDEAL_RADIUS * IDEAL_RADIUS;

This is what the ported code should look like:

// Calculate pressure
float p = 0.0f;
float pnear = 0.0f;
for (int a = 0; a < particle.neighborCount; a++)
{
    Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index];
    float distanceSq = relativePosition.LengthSquared();

    //within idealRad check
    if (distanceSq < IDEAL_RADIUS_SQ)
    {
        particle.distances[a] = (float)Math.Sqrt(distanceSq);
        //if (particle.distances[a] < Settings.EPSILON) particle.distances[a] = IDEAL_RADIUS - .01f;
        float oneminusq = 1.0f - (particle.distances[a] / IDEAL_RADIUS);
        p = (p + oneminusq * oneminusq);
        pnear = (pnear + oneminusq * oneminusq * oneminusq);
    }
    else
    {
        particle.distances[a] = float.MaxValue;
    }
}

Now for the code that applies the forces. Here is the code we’re focusing on:

// Now actually apply the forces
float pressure = (p - 5F) / 2.0F; //normal pressure term
float presnear = pnear / 2.0F; //near particles term
float changex = 0.0F;
float changey = 0.0F;
for(int a = 0; a < neighbors.size(); a++) {
    Integer n = (Integer)neighbors.get(a);
    int j = n.intValue();
    float vx = xs[j]-xs[i];//liquid[j].m_sweep.c.x - liquid[i].m_sweep.c.x;
    float vy = ys[j]-ys[i];//liquid[j].m_sweep.c.y - liquid[i].m_sweep.c.y;
    if(vx > -idealRad && vx < idealRad && vy > -idealRad && vy < idealRad) {
        if(vlen[a] < idealRad) {
            float q = vlen[a] / idealRad;
            float oneminusq = 1.0f-q;
            float factor = oneminusq * (pressure + presnear * oneminusq) / (2.0F*vlen[a]);
            float dx = vx * factor;
            float dy = vy * factor;
            float relvx = vxs[j] - vxs[i];
            float relvy = vys[j] - vys[i];
            factor = visc * oneminusq * deltaT;
            dx -= relvx * factor;
            dy -= relvy * factor;
            xchange[j] += dx;
            ychange[j] += dy;
            changex -= dx;
            changey -= dy;
        }
    }
}
xchange[i] += changex;
ychange[i] += changey;

Once again, I’m going to remove the early exit check. Here’s what the ported code looks like:

// Apply forces
float pressure = (p - 5f) / 2.0f; //normal pressure term
float presnear = pnear / 2.0f; //near particles term
Vector2 change = Vector2.Zero;
for (int a = 0; a < particle.neighborCount; a++)
{
    Vector2 relativePosition = _scaledPositions[particle.neighbors[a]] - _scaledPositions[index];

    if (particle.distances[a] < IDEAL_RADIUS)
    {
        float q = particle.distances[a] / IDEAL_RADIUS;
        float oneminusq = 1.0f - q;
        float factor = oneminusq * (pressure + presnear * oneminusq) / (2.0F * particle.distances[a]);
        Vector2 d = relativePosition * factor;
        Vector2 relativeVelocity = _scaledVelocities[particle.neighbors[a]] - _scaledVelocities[index];

        factor = VISCOSITY * oneminusq * DT;
        d -= relativeVelocity * factor;
        _delta[particle.neighbors[a]] += d;
        change -= d;
    }
}
_delta[index] += change;

Now we have to move the particles by updating their positions and velocities:

// Move particles
for (int i = 0; i < _numActiveParticles; i++)
{
    int index = _activeParticles[i];
    Particle particle = _liquid[index];

    particle.position += _delta[index] / MULTIPLIER;
    particle.velocity += _delta[index] / (MULTIPLIER * DT);
}

Spatial Partitioning

We need to implement the grid that will allow us to quickly find neighboring particles by breaking space up into small cells.

Here is a graphic that should help visualize the idea behind spatial partitioning. The blue dots are the fluid particles. The green cell represents a cell that contains a particle whose neighbors need to be found. The orange cells are the neighboring cells.

Instead of trying to define a 2d array of grid cells large enough to contain all the liquid (and trying to keep liquid inside that region of space), I’m going to use a dynamic meshing approach that I stumbled across in Dynamic Meshing for Material Point Method Computations

This will allow us to only have grid cells where particles exist (the cells are the blue rectangles):

This type of partitioning is important, because it means liquid can be used anywhere. The creation/removal of cells adds a little overhead, but it’s worth it.

The dynamic grid approach uses dictionaries instead of arrays to hold the cells. All a grid cell consists of is an X and Y coordinate, and a list of the particles in that cell. Lets add the grid structure to FluidSimulation as a member variable:

private Dictionary<int, Dictionary<int, List<int>>> _grid;

We also should initialize it in the FluidSimulation constructor:

_grid = new Dictionary<int, Dictionary<int, List<int>>>();

The grid cells need a size, so lets add one to FluidSimulation (we’re using a small number to match the scale of the Box2D world):

public const float CELL_SIZE = 0.5f;

We need to write the methods that will let us convert a position into a grid coordinate. To find the grid coordinates, we just divide the position by the cell size, and drop the fractional portion. Add the following two methods to FluidSimulation:

private int getGridX(float x) { return (int)Math.Floor(x / CELL_SIZE); }
private int getGridY(float y) { return (int)Math.Floor(y / CELL_SIZE); }

Since the particles can move into and out of cells, we need to update the grid for every particle that was moved. Add the following variables to Particle to keep track of the grid coordinates:

public int ci;
public int cj;

Add the following code just after the particle’s position and velocity gets updated in applyLiquidConstraints():

// Update particle cell
int x = getGridX(particle.position.X);
int y = getGridY(particle.position.Y);

if (particle.ci == x && particle.cj == y)
    continue;
else
{
    _grid[particle.ci][particle.cj].Remove(index);

    if (_grid[particle.ci][particle.cj].Count == 0)
    {
        _grid[particle.ci].Remove(particle.cj);

        if (_grid[particle.ci].Count == 0)
        {
            _grid.Remove(particle.ci);
        }
    }

    if (!_grid.ContainsKey(x))
        _grid[x] = new Dictionary<int, List<int>>();
    if (!_grid[x].ContainsKey(y))
        _grid[x][y] = new List<int>(20);

    _grid[x][y].Add(index);
    particle.ci = x;
    particle.cj = y;
}

All this code does is update the grid cells by comparing the particle’s previous cell coordinates (ci and cj) to the new coordinates (x and y). If the coordinates are the same, no update needs to happen. If the coordinates are different, the particle is removed from the list of particles in the previous cell, and added to the list of particles in the new cell. If a cell is empty after removing a particle, the cell is removed. If a cell doesn’t exist when trying to add a particle, it is created.

Now that the grid structure is set up, we need to go back and implement the neighbor search.

Neighbor Search

First, we need to define how many neighbors we should worry about. In my experience, 75 seems to be enough. So lets set that as a constant in FluidSimulation:

public const int MAX_NEIGHBORS = 75;

Now that we know how many neighbors there can be, lets initialize the neighbors and distances arrays in the Particle constructor:

distances = new float[FluidSimulation.MAX_NEIGHBORS];
neighbors = new int[FluidSimulation.MAX_NEIGHBORS];

We need to write a new method in FluidSimulation called findNeighbors().

private void findNeighbors(Particle particle)
{
    particle.neighborCount = 0;
    Dictionary<int, List<int>> gridX;
    List<int> gridY;

    for (int nx = -1; nx < 2; nx++)
    {
        for (int ny = -1; ny < 2; ny++)
        {
            int x = particle.ci + nx;
            int y = particle.cj + ny;
            if (_grid.TryGetValue(x, out gridX) && gridX.TryGetValue(y, out gridY))
            {

                for (int a = 0; a < gridY.Count; a++)
                {
                    if (gridY[a] != particle.index)
                    {
                        particle.neighbors[particle.neighborCount] = gridY[a];
                        particle.neighborCount++;

                        if (particle.neighborCount >= MAX_NEIGHBORS)
                            return;
                    }
                }
            }
        }
    }
}

This code loops through both the current particle’s cell and the 8 surrounding cells. The particles in those cells have their indices added to the current particle’s neighbors array. If the neighborCount is greater than or equal to the MAX_NEIGHBORS, it stops finding neighbors.

We need to place a call to this method in applyLiquidConstraints(), just before the pressure is calculated:

# ...
int index = _activeParticles[i];
Particle particle = _liquid[index];

// Find neighbors
findNeighbors(particle);

// Calculate pressure
float p = 0.0f;
float pnear = 0.0f;
# ...

Simple Rendering

I will describe how to render the fluid in another post. For now I just want to render the particles as dots. I’m going to add a member variable to FluidSimulation:

private Texture2D _pixel;

and create the texture in the FluidSimulation constructor:

_pixel = new Texture2D(_spriteBatch.GraphicsDevice, 1, 1);
_pixel.SetData<Color>(new[] { Color.White });

Now I just need to add a draw() method that loops through all the active particles and draws them:

public void draw()
{
    _spriteBatch.Begin();

    for (int i = 0; i < _numActiveParticles; i++)
    {
        Particle particle = _liquid[_activeParticles[i]];

        _spriteBatch.Draw(_pixel, particle.position * _scale, new Rectangle(0, 0, 2, 2), Color.LightBlue, 0f, new Vector2(1, 1), 1f, SpriteEffects.None, 0f);
    }

    _spriteBatch.End();
}

Don’t forget to call _fluidSimulation.draw() in the Game1 Draw() method.

Running the Simulation

Now that everything is set up to run the simulation, we just have to give FluidSimulation an update() method and handle particle creation.

The update() method is pretty straightforward:

public void update()
{
    MouseState mouseState = Mouse.GetState();

    _mouse = new Vector2(mouseState.X, mouseState.Y) / _scale;
    if (mouseState.LeftButton == ButtonState.Pressed)
        createParticle();

    applyLiquidConstraints();
}

The createParticle() method:

public void createParticle(int numParticlesToSpawn = 4)
{
    IEnumerable<Particle> inactiveParticles = from particle in _liquid
                                                where particle.alive == false
                                                select particle;
    inactiveParticles = inactiveParticles.Take(numParticlesToSpawn);

    foreach (Particle particle in inactiveParticles)
    {
        if (_numActiveParticles < MAX_PARTICLES)
        {
            Vector2 jitter = new Vector2((float)(_random.NextDouble() * 2 - 1), (float)(_random.NextDouble()) - 0.5f);

            particle.position = _mouse + jitter;
            particle.velocity = Vector2.Zero;
            particle.alive = true;
            particle.ci = getGridX(particle.position.X);
            particle.cj = getGridY(particle.position.Y);

            // Create grid cell if necessary
            if (!_grid.ContainsKey(particle.ci))
                _grid[particle.ci] = new Dictionary<int, List<int>>();
            if (!_grid[particle.ci].ContainsKey(particle.cj))
                _grid[particle.ci][particle.cj] = new List<int>();
            _grid[particle.ci][particle.cj].Add(particle.index);

            _activeParticles.Add(particle.index);
            _numActiveParticles++;
        }
    }
}

Basically all createParticle does is look for dead particles. If it finds some, it reactivates numParticlesToSpawn-many of them and adds them to the list of active particles.

It also has to create grid cells if those cells don’t exist at the position of the particle.

You should now be able to create particles with your mouse. Currently, there’s no boundaries. Which means that if you applied a gravitational force to the liquid, it would just fall off the screen. You can tell the simulation is working by how the particles form droplets/pools after being created.

In the next post, I’ll describe how to add some simple boundaries and apply a gravitational force use the Task Parallel Library to run the simulation asynchronously on multiple threads.

Thanks for reading! Please leave a reply if you have any questions or comments.