# 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

• an integer representing the maximum amount of particles the simulation can handle
• an integer representing the number of active particles
• an array holding all of the particles
• a list of integers representing the indices of active particles
• a float representing the particles’ radius
• a float representing the fluid’s viscosity
• a float representing the change in time

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;

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
float vlensqr = (vx * vx + vy * vy);
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:

• In my experience, removing the early exit check didn’t have a noticeable effect on the simulation or its performance, so I’m just going to remove it.
• The value of `idealRad * idealRad` is never going to change, so we can move that to a constant.

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
{
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;
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>();

```            _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.