Fluid Simulation Using Implicit Particles
Advanced Game Programming
Dan Englesson
Joakim Kilby
Joel Ek
December 20, 2011
Abstract
This report covers the implementation of a fluid simulation using a hybrid method.
The method used is heavily based on the contents of the book Fluid Simulation for
Computer Graphics by Robert Bridson [1] and is commonly referred to as the PIC/FLIP
method. In the method, the fluid surface is tracked using particles and its mass conserved
by enforcing zero divergence on the deforming velocity field. We store the velocity field
on a staggered grid as described in [2], which greatly helps to fulfill the mass-conserving
criteria. We also reconstruct the actual fluid surface by evaluating the Improved Blob-
bies signed distance function introduced in [3] on a regular grid followed by a standard
implementation of Marching Cubes.
We have tested our implementation using five different simulation cases and obtain
great results, all of which are shown in the report together with benchmark data.
1 Introduction
Fire, water and air, once thought of as
three of the four fundamental elements
making up all known substances, seem to
have more in common than once thought.
All three are different types of fluids, or
compounds which are able to flow and de-
form under the influence of commonly oc-
curring forces such as gravity.
The complex behavior of fluids has al-
ways been fascinating to man and it is only
in the latest decades that we have come
to fully understand the complex equations
governing their dynamic behavior. Fluid
simulation is greatly useful since it can aid
in the design of buildings and construc-
tions but also since it allows one to simu-
late an event which never has to happen.
For instance, the flooding of a city, a great
explosion or anything else that would fit in
a modern feature film.
2 Background
In order to understand how the differ-
ent steps of the method work, some back-
ground information is needed. This sec-
tion contains an explanation of the funda-
mental equations, a brief repetition of rel-
evant vector calculus and a short descrip-
tion of the infamous Stable Fluids method.
It also explains the difference between the
Eulerian and the Lagrangian viewpoints
as well as the benefits of using a staggered
grid.
2.1 Governing equations
The two equations governing fluid flow are
called the Navier-Stokes equations. The
fundamental equation is given in equation
1 which describes how a velocity field V
evolves over time. Naturally, the time
derivative of this velocity field must be
subject to external forces. This is modeled
by the F term in equation 1 and is com-
monly representing gravity, wind or user
interaction. The ��V2V term is the diffu-
sion term, which models the viscosity of
the fluid.
Viscosity is a commonly used term in
fluid simulation which simply means the
tendency of the fluid to resist flow and is
observed as its thickness. For example,
water flows easily and therefore has low
viscosity. Honey is thick and resists flow
which is why its viscosity is high. The
amount of viscosity is directly controlled
by the scalar ��. It should be noted that
this term is commonly discarded when
simulating water. The amount of numeri-
cal dissipation introduced by interpolation
errors is often a good replacement for this
term, which otherwise would be solved by
a Poisson equation.
∂V
∂t
= F + ��V2V - V · VV -
Vp
��
(1)
V · V = 0
(2)
The V · VV term is the self-advection
term which models how the velocity field
flows within itself. This can seem a bit
strange but it is only natural as the ve-
locity field represents the movements of
the mass of a fluid. The final term,
∇p
��
,
is the subtraction of the pressure gradi-
ent Vp over the density ��. When com-
bining the entire equation with the con-
straint in equation 2, this term is replaced
by a subtraction of a pseudo-pressure gra-
dient. This as a result of any gradient
of a scalar being curl-free and therefore
removed when enforcing a divergence-free
velocity field due to the Helmholtz-Hodge
decomposition as stated in [4].
2.2 Representing fluids
There are two common ways of represent-
ing fluids which are based on completely
different viewpoints. One way is to think
of the fluid as a collection of tiny parti-
cles, or atoms, which together make up
the volume of the fluid. This represen-
tation is a Lagrangian representation and
has given rise to a whole family of purely
1
particle-based fluid simulations of which
smoothed particle hydrodynamics, or SPH
is the most common. These methods can
be beneficial as their definition is intuitive
and understandable but have a great flaw.
The calculations are based on the assump-
tion that the particle density always is suf-
ficient. This is not the case as particles are
able to spread out, allowing for regions
with low density. In these regions, the
calculations become less accurate as there
simply is not enough information stored in
the nearby particles.
Another way of representing fluids is
from the Eulerian point of view. For this
representation, the fluid is observed from
a domain of interest which is divided into
a number of cells. This forms a grid,
which is why the Eulerian representation
also is known as the grid-based represen-
tation. Grid-based representations do not
suffer from the same problems with low-
density regions as particle-based represen-
tations do. Instead, other problematic is-
sues arise. In order for the fluid surface to
be detailed enough, the grid must consist
of a multitude of cells. It is not uncommon
for purely grid-based fluid simulations to
use millions of cells. This effect is due
to the Nyquist criterion which states that
the sampling frequency of a signal must be
twice that of the highest frequency com-
ponent. In other words, the simulation
is only guaranteed to capture surface de-
tails larger than the magnitude of two cell
widths. This is why the number of cells
must be extremely large for purely grid-
based methods.
A fluid simulation can also use a hy-
brid method, which combines the sur-
face capturing possibilities of a particle-
based method and couples it with an
auxiliary grid to enforce accurate mass-
conservation. Such a method is explained
in detail in the method section.
2.3 Operators and vector fields
In order to understand how the veloc-
ity field is evolved and how the mass-
conserving constraint is maintained, the
following three operators need to be un-
derstood. Equation 3 shows the gradient
operator. It operates on a scalar field and
yields a vector field with each of the par-
tial derivatives as its components. When
the operator is applied to the cells of a
grid structure, it should be interpreted as
the local change of the scalar field q. In
other words, the exchange of the property
described by the scalar q.
Vq =
( ∂q
∂x
,
∂q
∂y
,
∂q
∂z
)
(3)
Equation 4 shows the divergence oper-
ator. Divergence is related to the gradi-
ent operator and involves the same partial
derivatives but instead of measuring each
derivative separately, it measures the total
local change. As shown below, it operates
on a vector field and yields a scalar field.
The operator should be interpreted as the
total increase or decrease of any property
being transported by the vector field u.
V · u =
∂ u
∂x
+
∂ u
∂y
+
∂ u
∂z
(4)
The Laplace operator is the most com-
plicated of the three operators discussed
here. It can be seen in equation 5 and op-
erates on a scalar field q. Even though the
expression involves partial derivatives of
the second order which might seem com-
plex, it is actually the divergence operator
applied to the gradient of a scalar field. In
other words, a concise way of stating a se-
quential application of the two operators
previously defined. When applied to the
cells of a grid structure, it should be inter-
preted as the total interchange of material
between all adjacent cells.
V2q =
∂2q
∂x2
+
∂2q
∂y2
+
∂2q
∂z2
(5)
2.4 The staggered grid
In order to apply the continuous operators
to a discrete grid structure, the operators
2
first need to be differentiated. This re-
quires an introduction to a special type of
grid, the staggered marker and cell (MAC)
grid, introduced by Foster and Metaxas in
[2]. The principal difference compared to
a common regular grid is that the com-
ponents of the velocity field are separated
and offset. This is shown in figure 1 which
shows that the sample points of the veloc-
ity field coincide with the cell face bound-
aries. Please note that all scalars which
are to be transported by the velocity field
have their sample points located at the
center of a cell. This is important for the
definition of the differentiated operators.
Figure 1: The staggered grid
The benefit of using this structure is
that the evaluation of the discrete counter-
parts of the operators previously defined
becomes greatly simplified. The structure
introduces the curious half index for which
a positive half index is the sample value
stored at the face shared between a cell
and its consequent cell.
Consider the discrete gradient shown in
equation 6 where s is the size of a cell. The
gradient is applied to the scalar q, which is
stored in the center of the cells. The eval-
uation position must therefore be exactly
between two adjacent cells. Please note
that this evaluation is made component-
wise and results in three different scalars,
each with different evaluation positions.
Also note that the evaluation positions co-
incide with the sample locations of the ve-
locity field components.
(Vq)i+1/2,j,k =
(qi+1,j,k−qi,j,k
s
)
(Vq)i,j+1/2,k =
(qi,j+1,k−qi,j,k
s
)
(Vq)i,j,k+1/2 =
(qi,j,k+1−qi,j,k
s
)
(6)
The definition of the discrete divergence
operator is shown in equation 7. As with
the discrete gradient, s denotes the cell
size. Divergence operates on component
pairs of the velocity field and its evalua-
tion position must therefore be the in the
center of a cell.
(V · u)i,j,k
=
ui+1/2,j,k− ui−1/2,j,k
s
+
ui,j+1/2,k− ui,j−1/2,k
s
+
ui,j,k+1/2− ui,j,k−1/2
s
(7)
The discrete Laplace operator is as
stated a combination of the gradient and
the divergence operator. In order to un-
derstand how it is defined on the stag-
gered grid, first consider the application
of the gradient operator. As previously
stated, the evaluation position of the gra-
dient is between two adjacent cells. If the
divergence operator is applied to those six
scalars, we obtain equation 8. As two op-
erators are applied in effect, the scaling
factor of s2 is natural.
(V2q)i,j,k
=
qi+1,j,k+qi−1,j,k
s2
+
qi,j+1,k+qi,j−1,k
s2
+
qi,j,k+1+qi,j,k−1
s2
-
6qi,j,k
s2
(8)
Please note that both the discrete
Laplace operator and the discrete diver-
gence operator are evaluated at the cell
center. Since Poisson equations have the
form shown in equation 9, using a stag-
gered grid is clearly beneficial for solving
these equations. Especially since it allows
for the evaluation of central differences us-
ing a width of only a single cell size.
V · u = V2q
(9)
3
2.5 Stable Fluids
The Stable Fluids method was introduced
by Jos Stam in 1999 [4]. It was the first
unconditionally stable method for fluid
simulation and introduced the concept of
semi-Lagrangian advection. The Stable
Fluids method calculates an approximate
solution to the Navier-Stokes equations;
equation 1 and 2.
A solution to equation 1 is obtained by
sequentially calculating the contribution
from each part of the equation, where for
each step the input is given by the output
of the previous step.
In order to enforce the constraint given
by equation 2, the final velocity field ob-
tained in the solution method is projected
onto its divergence-free part.
The entire process is illustrated in fig-
ure 2 where w0 is the velocity field from
the previous time step and w4 is the final,
divergence-free, velocity field.
Figure 2: The Stable Fluids method
3 Method
3.1 Mesh conversion
In order to represent polygon mesh ob-
jects as fluids or solid objects for simula-
tion purposes, a conversion must be made
from polygons into a voxel representation.
To achieve this, a polygon object is read
from an .obj file and as each triangle of the
object is processed, its vertex coordinates
are transformed into discrete grid coordi-
nates. The voxel corresponding to those
grid coordinates, the voxel in which the
vertex lies, is then set as the cell-type the
object is to represent, e.g. fluid or solid.
To avoid problems with triangles that
are significantly larger than the grid-
resolution, a number of points are se-
lected, in a random fashion, on the trian-
gle surface and transformed into grid co-
ordinates as well.
This method creates a voxelized shell,
or surface, resembling the original poly-
gon object. However for simulation, a sur-
face is insufficient and the voxelized model
must be filled to create a volume.
The volume is created by filling out the
gaps between voxels in the representation.
This is achieved by selecting a voxel which
is marked with the correct cell-type and
then stepping along the x, y and z direc-
tions in turn until another voxel is found
which has the same cell-type. All voxels
between the two are then set as the correct
cell-type.
This method works under the assump-
tion that all polygon objects are Mani-
fold surfaces; there can be no holes in the
mesh, no self-intersection of polygons and
no interior shells within the mesh.
3.2 Hybrid particle methods
There exist two main approaches for sim-
ulating fluid effects in computer graph-
ics as mentioned earlier, namely the use
of Lagrangian particles or Eulerian grids.
The two approaches have their strengths
and weaknesses, the Lagrangian particles
are very good with the advection part but
have problem with the pressure and in-
compressibility constraint. Fortunately,
the Eulerian grid approach is excellent in
solving the pressure and incompressibil-
ity constraint but due to interpolation er-
rors from the semi-Lagrangian advection,
the method has problem with the advec-
tion part. One can see that where the
Lagrangian approach has its difficulties,
the Eulerian approach is very good and
by combining these two methods by let-
ting Lagrangian particles handle the ad-
vection part and the Eulerian grid han-
dle the pressure and incompressibility con-
straint, a better simulation of water effects
can be achieved.
Many different hybrid approaches exists
such as particle level sets, The Particle in
Cell (PIC) and the Fluid Implicit Parti-
4
cle (FLIP) method, where the two latter
methods will be further discussed in this
paper. The particle level set method was
introduced by Foster and Fedkiw in [5]
and the PIC method was introduced as
early as 1963 by Harlow in [6] and was
later improved by Brackbill and Ruppel
[7] with the FLIP method in 1986. The
FLIP method where then introduced to
incompressible flow by Zhu and Bridson
in 2005 [3] and is today a state-of-the-art
technique for fluid simulation.
3.2.1 The PIC/FLIP method
Even though the PIC and FLIP methods
are basically the same method but with
a slight difference in the velocity update,
the fluid characteristics differ quite a bit.
A fluid that only uses PIC is more viscous
than a FLIP fluid, which is due to numeri-
cal dissipation because of the double inter-
polation of the velocity, but more on that
in a bit. The FLIP method however, has
little viscosity and is therefore very well-
suited for water effects but unwanted visi-
ble noise on the surface is present. By lin-
early combining the PIC and FLIP meth-
ods, one obtain a fluid which has little vis-
cosity and is free of surface noise.
The severe numerical dissipation of the
PIC method is caused by interpolation.
The particle velocities are transferred to
the grid through interpolation which in-
troduces some smoothing and then the
smoothed velocity values are being used to
solve the new velocities with the Navier-
Stokes equations and are then interpo-
lated back to the particles replacing the
old velocities. As mentioned, due to this
excessive double interpolation, a PIC fluid
will appear more viscous than a FLIP
fluid. Brackbill and Ruppel solved this
problem by instead of interpolating the
newly calculated velocities to the particles
from the grid and replacing the velocities,
they interpolated the change in velocity
and added it to the already existing parti-
cle velocities. By doing so, only smoothing
from one interpolation is performed and
thus the smoothing is not accumulated as
in the PIC method and makes the FLIP
method almost free of numerical dissipa-
tion. Below follows a stepwise method-of-
solution for a PIC/FLIP fluid simulation
and a more in-depth view of the different
steps.
1. Initialize the grid and the particle po-
sitions and velocities.
2. Transfer particle velocities to a stag-
gered grid.
3. FLIP: Save a copy of the grid veloci-
ties.
4. Calculate and apply external forces.
5. Enforce the Dirichlet boundary con-
dition.
6. Classify all voxels as fluid, solid or air.
7. Calculate the pseudo-pressure gradi-
ent using a Preconditioned Conjugate
Gradient method.
8. FLIP: Update particle velocities by
subtracting the new grid velocities
from the saved grid velocities, inter-
polate and add the difference to the
particle velocities.
9. PIC: Update particle velocities by in-
terpolating and replacing the old ve-
locities with the new grid velocities to
the particles.
10. PIC/FLIP: Update the particle veloc-
ities by taking the interpolated FLIP
and PIC velocities from 8. and 9. and
for each particle weigh the PIC and
FLIP velocities with a factor �� be-
tween zero and one.
11. Update the particle positions with an
ODE solver with the newly created
velocities.
5
3.2.2 Initialize particles
Every voxel that is classified as fluid will
have seeded particles inside. The parti-
cles are seeded in a jittered grid, much like
super-sampling in a renderer, where eight
particles in each fluid voxel is jittered by
randomly placing the particles in a 2x2x2
sub-grid. A higher number of particles can
be used in each voxel but does not neces-
sarily give significantly better results but
might instead slow down the calculations.
Bridson stated in [3] that eight particles
per fluid cell was enough to yield good re-
sults.
3.2.3 Transfer particles to grid
Since the particles are randomly placed
in space and the calculations of the pres-
sure and incompressibility are performed
on the grid, some kind of interpolation is
needed to transfer the nearby particle ve-
locities to the grid. The grid velocities are
updated from the nearby particles through
trilinear interpolation of the weighted av-
erage particle velocities that lies in a cube
twice the grid cell width where the center
is at the grid-velocity component.
3.2.4 Calculate external forces
The external forces such as gravity and
forces from any user interaction are added
to the velocities by simple Euler inte-
gration as seen in equation 10, where V
stands for velocity, F is the external forces
and ∆t is the time step.
Vnew = Vold + F∆t
(10)
A simple Euler integration is enough
since this task is unconditionally stable for
reasonably small ∆t. This is because there
exists no feedback loop between the exter-
nal forces and the velocity field.
3.2.5 Enforce Dirichlet boundary
condition
The Dirichlet boundary condition states
that there should be no flow into or out
of solid cells to which has the normal n,
see equation 11 where V is the velocity in
the fluid cell and n is the normal to the
neighboring solid.
V · n = 0
(11)
If a fluid cell has a solid neighboring cell,
the velocity components are checked and
if any of the velocity components point to-
wards a neighboring solid cell, the velocity
is projected to go along with the surface of
the solid cell. Since the simulation is per-
formed on a staggerd MAC-grid and the
velocities are divided into separate com-
ponents, the projection of the velocity is
very simple. One just has to set it to zero
if the component points into a neighboring
solid cell.
3.2.6 Classify voxels
Firstly, the type of every non-solid voxel
is cleared and set to air. Then voxels con-
taining one or more particles are set to
fluid. Solid voxels are set to be solid from
the beginning of the simulation and and
do not change their type.
3.2.7 Conserving mass
A mass-conserving velocity field is synony-
mous to a divergence-free velocity field.
Conserving the mass of a fluid is there-
fore equal to enforcing zero divergence on
the velocity field transporting mass.
This part of a simulation step is the
most computationally heavy one as it in-
volves solving a Poisson equation. Re-
call the previously stated Laplace opera-
tor and how it was defined as a summation
of the partial derivatives of the second or-
der. It should be emphasized that there
is a dependency between every pair of two
adjacent cells and that the solution to a
Poisson equation in effect means solving a
massive system of linear equations. Before
this task can be undertaken, a derivation
of how the Poisson equation is obtained is
required.
6
The Helmholtz-Hodge decomposition
states that any vector field can be ex-
pressed as a divergence-free part Vdf and
a curl-free part Vcf as shown in equation
12.
V = Vdf + Vcf
(12)
Furthermore, vector calculus states that
the gradient of any scalar fields is, by defi-
nition, curl-free. We can therefore replace
the curl-free part of our velocity field with
the gradient of an unknown scalar field,
Vq. This is shown in equation 13.
V = Vdf + Vq
(13)
As with any equation, applying an oper-
ator to both sides of the equality sign does
not change equality. Thus, it is allowed
to apply the divergence operator to both
sides. Since divergence is a linear opera-
tor, V·(Vdf +Vq) is equal to V·Vdf +V·Vq
and we get equation 14.
V · V = V · Vdf + V·Vq
(14)
Equation 14 can be simplified as we in
the decomposition step defined Vdf as be-
ing divergence-free. The divergence oper-
ator applied to Vdf must therefore be equal
to zero. Also, since we defined the Laplace
operator as the consequent application of
a gradient operator and a divergence op-
erator, we get equation 15.
V · V = V2q
(15)
Equation 15 is a Poisson equation in
which q is the unknown scalar field which
solves the equation. By rearranging equa-
tion 13 we obtain equation 16 which
clearly shows that we can enforce mass-
conservation on the velocity field by find-
ing the unknown scalar field q and sub-
tracting its gradient, Vq.
Vdf = V - Vq
(16)
Equation 15 is the equation which needs
to be solved in order to find the scalar field
q, often referred to as the pseudo-pressure.
The left-hand side of the equation can be
calculated simply by evaluating the diver-
gence of the velocity field. The right-hand
side contains the unknown scalar field and
the Laplace operator. The definition of
the operator is shown in equation 8 and
states that adjacent values of the unknown
scalar field should be used in the calcula-
tions. Naturally, this is not possible as the
scalar field is unknown. However, the co-
efficients are fully known as they only de-
pend on the local configuration of a cell.
That is, if there are solid cells directly ad-
jacent or if the fluid is free to exchange
material through every cell face.
The definition of the Laplace operator
becomes somewhat different for boundary
fluid cells as the Dirichlet boundary con-
dition states that there should be no flow
through solid boundaries. The result on
the Laplace operator is that the coefficient
for the solid cell (when considering the ad-
jacent fluid cell) is set to zero and that the
central coefficient is increased by one.
As the coefficients are calculated and
stored, a massive system of linear equa-
tions is obtained. The system is on the
form Ax = b where A is the coefficient ma-
trix and b is the divergence of every cell.
The number of equations in this system is
directly related to the number of cells in
the simulation grid. If there are w · h · d
cells, the size of the coefficient matrix is
(w · h · d)2.
Since there can be at most six other cells
adjacent to every cell, this system is very
sparse and can be solved efficiently by it-
erative methods taking advantage of this
property. The memory requirement of the
sparse coefficient matrix is w·h·d·4. This
can be compared with the virtual size of
the matrix previously stated.
In order to solve the Poisson equation,
Bridson [1] recommended the precondi-
tioned conjugate gradient method with a
preconditioner of the modified incomplete
Cholesky factorization type. This method
has fast convergence and is able to operate
on a sparse system. An implementation
7
of the PCG method was made with the
pseudo-code in Bridson��s book as a guide.
Though it should be noted that any iter-
ative method that does not require an ex-
plicit representation of the A matrix could
be used, such as the Jacobi iterative tech-
nique.
When the Poisson equation is solved,
the gradient of the resulting scalar field
is subtracted from the velocity field and
mass-conservation is guaranteed.
3.2.8 Update particle velocities
The particles need to be updated with
the newly calculated velocities which are
stored on the MAC-grid. This is done by
trilinearly interpolating the velocities of
the eight neighboring grid-velocities to the
particle and, for PIC, update the velocity
with the new velocity or for FLIP, inter-
polate the change in velocity and add it
to the existing particle velocity. A linear
combination of both PIC and FLIP can
be used to get a low viscosity, water-like,
fluid with no surface noise. This can be
seen in equation 17 where unew
p
is the new
particle velocity and �� is the PIC blend-
ing factor. The factor determines how
much PIC, viscosity or numerical dissipa-
tion there should be where 1.0 is pure PIC
and 0.0 is pure FLIP. The lerp() functions
represents the trilinear interpolation func-
tion.
unew
p
= �� · lerp( unew
grid, xp)
+ (1 - ��)[ uold
p
+ lerp(∆ ugrid, xp)]
(17)
3.2.9 The CFL condition
The CFL condition states that the a parti-
cle should always move less than one grid-
cell in each sub step. It is done by taking
the cell-width and dividing it by the maxi-
mum velocity in the grid to get a stabledt.
The stabledt is then compared to the ac-
tual time step dt and if it is larger than dt,
the stabledt is set to dt. The particles are
then advected in six sub steps until it has
reached dt. It is common to have around
five sub steps.
3.2.10 Update particle positions
The particles positions are advected with
a Runge Kutta 2 ODE solver, which is
stable as opposed to a simple Euler ODE
solver. However, the particles can occa-
sionally penetrate solid boundaries due to
errors in the RK2 solver. This can cause
the particle that has penetrated a solid
voxel to become stuck. To fix this prob-
lem, the particles that have penetrated a
solid voxel are moved back in the normal
direction to half the cell-width outside of
the solid voxel.
3.3 Intermediate storage
The fluid solver is divided into two differ-
ent parts; a simulation program which is
responsible for all calculations regarding
the movement of the fluid and a visualiza-
tion program which takes simulation data
as input and produces .obj files, meshes,
from the data. The visualization program
also has the capability to directly visualize
particles or surfaces.
As the simulation program progresses,
it sequentially writes the positions of the
particles for the current time step into a
binary file. The reason for writing only
particle positions comes from the fact that
the number of particles are very large and
that the positions are all that is necessary
to reconstruct surfaces as they determine
the position of the fluid.
The position of a particle is repre-
sented by three floating point numbers,
the Cartesian coordinates x, y and z. As
each floating point number occupies four
bytes of memory, the position of a particle
occupies 12 bytes.
Table 2, in appendix A, shows the size
of the simulation data for a few different
amounts of particles, 250 frames are saved
in each of the files.
8
3.4 Surface reconstruction
When a simulation has been performed
and stored on disk, the surface must be
reconstructed in order to be rendered in
a standard rendering pipeline. As most
rendering pipelines are optimized to han-
dle triangles, this is often reduced to tri-
angulating the simulation data. This is
no easy task as there are no trivial solu-
tions to the problem. Surface reconstruc-
tion is commonly achieved in two distinct
passes. First, a scalar field is initialized
from the underlying data on a regular grid
using a function. The only requirement
on this function is that it maps the data
set onto a real-valued scalar. Note that
real-valued implies that the function can
assume both positive and negative values
and that scalar means single value. In the
second pass, a certain value (an iso-level)
of the calculated scalar field is chosen to
be visualized. The reconstructed surface
is created where the function assumes this
value.
3.4.1 The signed distance function
When dealing with point clouds, the map-
ping function commonly measures dis-
tance to the nearest point or average dis-
tance to some of the nearest points. Zhu
and Bridson introduced a method called
Improved Blobbies in [3]. The method
calculates, for every grid corner, the av-
erage position and radius of the nearby
points and determines whether the grid
corner is located within the average radius
of the average position. This gives neg-
ative distances for corners that are close
to many points and positive distances for
those who have few points nearby. In ef-
fect, this forms a signed distance function
for which the definition can be seen in
equation 18. However, when there are no
points present within the fixed search ra-
dius, the function is undefined and this
can cause problems.
The function is defined at every grid
corner, which is why it makes sense to
evaluate it for every grid corner. How-
ever, this is highly inefficient as the time
complexity for such an implementation is
O(m · n) where n is the number of points
and m the number of grid corners. For
optimization, we set the search radius R
equal to three times the the grid spacing,
the particle radius ri to half the grid spac-
ing and iterate through the particles. This
makes the influence of every point easy to
determine as it only can affect the corners
located within its proximity. The time
complexity is thus reduced from O(m · n)
to O(n).
��( xg) = len( xg -��
i
wixi)-��
i
wiri (18)
The points are weighted as shown in
equation 19 using a well-shaped kernel
function, shown in equation 20. The ker-
nel function decays as the distance be-
tween the point and the grid corner grows,
being exactly zero at a distance equal to
the search radius, R. Every weight is nor-
malized with respect to all contributions
from nearby points. Implementing this
weighting functionality is therefore most
effectively done by accumulating contribu-
tions and their kernel values followed by a
final normalizing step.
wi =
k1(len( xg - xi)/R)
��
j k1(len( xg - xj)/R)
(19)
k1(s) = max(0,(1 - s2)3)
(20)
Since the kernel squares its input vari-
able and the input always is the norm of
a vector, we replace this kernel with an
optimization, avoiding square roots. The
improved kernel is shown in equation 22
and the modified weight function in equa-
tion 21. Note that the distance between xg
and xi is calculated as the square distance,
or dot product, which does not require a
square root operation.
wi =
k2(( xg - xi) · ( xg - xi)/R2)
��
j k2(( xg - xj) · ( xg - xj)/R2)
(21)
9
k2(s) = max(0,(1 - s)3)
(22)
3.4.2 Marching Cubes
A simple and elegant surface recon-
struction algorithm was introduced by
Lorensen and Cline in [8]. The new
method went under the name Marching
Cubes and has since become the indus-
try standard for iso-surface generation.
The method is based on evaluating a real-
valued function into a discrete scalar field
defined on a regular grid. It should be
noted that this grid is different from the
one used during the simulation pass. They
are not required to be of the same resolu-
tion, yet it is implied that this grid should
be of a resolution determined by the den-
sity of the tracker particles in the simula-
tion. If eight particles are seeded per fluid
cell, the optimal resolution for the surface
reconstruction grid will be that of the sim-
ulation grid according to the Nyquist cri-
terion.
The first step of the method is to eval-
uate the signed distance function for ev-
ery cell corner in the grid. Since our
signed distance function is not defined ev-
erywhere in the simulation domain, spe-
cial consideration has to be made about
where to evaluate the function. Cell cor-
ners where the function is undefined is
flagged as outside for the following steps
of the algorithm.
Following the first step is to, given a
certain iso-value, classify cell corners as
either inside or outside of the surface.
Since we are using a signed distance func-
tion which represents the signed distance
to the surface, the iso-value of interest
is zero. Determining which cell corners
that are outside thus corresponds to test-
ing whether the value is greater than zero
and vice versa.
Depending on how the corners are
classed as either inside or outside, a config-
uration is formed for every cell. This con-
figuration needs to be enumerated. Since
there are eight corners of a cell, there
Figure 3: The corner and edge numbering
can be exactly 256 different combinations,
each forming a unique case. The configu-
ration is perfect for storage in a single un-
signed byte which can be used as an index
into a case table, detailing which trian-
gles that form the surface. The notation
shown in figure 3 is used to encode a one
for every corner that is flagged as inside.
Figure 4: A unique case
For cells where only corners 0 and 1
are inside, the corresponding index will be
00000011 where the least significant bit is
located to the right. This is equal to 3
and thus case 3 is to be used from the
triangle case table. From the standard
implementation of Marching Cubes, this
corresponds to placing two triangles. One
between vertices 3, 1 and 8 and one be-
tween vertices 8, 1 and 9. This case is vis-
ible in figure 4 which clearly shows that a
surface is constructed which encapsulates
corners 0 and 1. The actual positions of
10
the vertices are determined by interpolat-
ing the scalar values located at the cor-
ners as shown in equation 23 where �� is
the blending factor between vertex i and
vertex j.
�� =
0 - ��i
��j - ��i
(23)
The original Marching Cubes algorithm
also details how the vertex normals can
be calculated from the values in the scalar
field. The gradient is simply evaluated
at every grid corner and interpolated to
the vertex positions. However, this is also
problematic as the function can be unde-
fined outside of the surface. Instead sur-
face normals are summed at the each ver-
tex and finally normalized. As the vertex
normals are calculated from the topology
of the surface rather than from the gradi-
ent of the scalar field, the resulting quality
of the vertex normals is lower. The two
calculation methods are not identical yet
both produce vertex normals of sufficient
quality for our purposes.
It should be noted that the original
Marching Cubes algorithm can create am-
biguous cases for which the generated sur-
face will have holes. The method used in
[9] details how these cases can be found
and corrected properly.
3.5 Exporting meshes
During the surface reconstruction pass, a
vertex list, a vertex normal list and a
triangle list is created. Vertices pend-
ing addition to the vertex list are com-
pared to the existing vertices. If they
are located within a small threshold value
(10−5), the old vertex is reused and tri-
angles are therefore able to share vertices.
This reduces the number of vertices in the
model and can also be used to disable the
creation of the unwanted tiny stretched
triangles, simply by raising the threshold
value. It should be noted that our im-
plementation uses a brute force method
which has a time complexity of O(n) in
which every new vertex is compared to
the already existing ones. An efficient im-
plementation should make use of a data
structure to reduce the time complexity
of this operation.
From these lists, the surface must be
stored in a highly accessible and standard-
ized format. We chose the object file for-
mat (.obj) from Alias Wavefront since it
is text-based and easy to implement. The
format is also compliant with the indexed
face set data structure, which consists of
the lists previously described.
4 Results
Our implementation was tested using five
different test cases. Test cases one and
two were designed to show how the viscos-
ity of the simulated fluid can be changed
simply by altering the PIC influence fac-
tor. All simulated cases are shown in ta-
ble 1, where additional simulation data are
shown in 2 in appending A.
Case Description
1
Dam break (FLIP)
2
Dam break (PIC)
3
Dam break with dragon
4
Falling cube
5
Splashing dragon
Table 1: List of simulation cases
All images shown below in appendix B,
were rendered in Autodesk 3DSMAX us-
ing simulation data from our implementa-
tion. For each case, a set of 250 frames
were simulated, reconstructed and ren-
dered.
5 Discussion
We have found that the use of hybrid
methods as opposed to pure Eulerian
methods for fluid simulation allows for
the capture of highly detailed behavior in
the fluid. Our method is not limited by
the grid resolution since the fundamen-
tal representation of the fluid is the set
of particles. This allows us to perform
11
simulations at a relatively low grid reso-
lution while still producing detailed sur-
faces, saving both time and memory. Our
method also captures the detail in splashes
with greater accuracy than Eulerian meth-
ods.
There is however a drawback to our
method; the surface reconstruction from
a set of particles, e.g. a point cloud, is
harder to perform than that of a level
set and the complexity of the generated
surfaces create a larger number of poly-
gons. The complexity of the surface can
be somewhat simplified by setting a limit
on the number of particles that must be
present in a voxel to generate a surface.
This will however reduce the detail in fine
features such as splashes and a trade-off
must be made between detail and com-
plexity.
Our implementation is able to perform
simulations in a relatively short time even
as the number of particles increase. This is
mainly due to clever data structures and a
mindfulness of how the cache memory and
caching works. As table 2 clearly shows,
the bulk of the time it takes to produce
a frame is spent on rendering and sur-
face reconstruction. We have noted that
as the complexity of our simulation data
increases, e.g. when large splashes occur,
the reconstruction phase takes longer and
produces more complex geometry with a
larger amount of polygons. This is some-
thing which we plan to remedy and a pos-
sible solution is mentioned in the next sec-
tion.
6 Improvements & future
work
A major area of improvement in our im-
plementation is the surface generation. A
first step would be to change the way in
which a particle may influence the grid
points. In the current implementation,
a particle influences grid points within a
sphere of a certain radius from its location.
This sphere-of-influence could be changed
into a more elliptical shape with an ori-
entation and size based on the local par-
ticle density. The most crucial improve-
ment this would yield is the ability to de-
fine sharper features in the surface. As
it stands today, extremely sharp features
are smoothed out because of the spherical
shape of influence from a particle. This is
described in more detail in [10].
We would also like to introduce post-
processing operations on the generated
surfaces, mainly to reduce the number of
unnecessary polygons and by doing so de-
creasing the rendering time. Amongst the
operations we would like to implement are;
Mesh-smoothing to reduce noise-like be-
havior which may be introduced by low
local particle density.
Decimation simply to reduce the num-
ber of polygons. We believe that mesh
decimation based on polygon area and
curvature would be very effective in reduc-
ing the amount of polygons whilst preserv-
ing the fine detail in the mesh.
As for the simulation part of our im-
plementation, we would like to introduce
seeding of foam particles at voxels where
the magnitude of the curl is above a cer-
tain threshold to emphasize the wild be-
havior of the fluid. We would also like to
introduce a two-way coupling in our sim-
ulator, meaning that the fluid may inter-
act with solid objects, or other fluids, in
a more complex way, e.g. having floating
solid objects or blending of oil and water.
We would also like to improve the grid
structure so that we can represent solid
objects which have a slope in a correct
way. Using the existing structure to do
this would result in a staircase-like behav-
ior as a result of the voxel representation.
Bridson has suggested a way of achieving
this in [1].
Finally we have considered implement-
ing an adaptive grid structure. In this
structure, the grid would only be defined
in the close vicinity to where the fluid ac-
tually is and it would follow the fluid as it
moves. We believe that this would reduce
12
the simulation time significantly.
References
[1] R. Bridson, Fluid Simulation for
Computer Graphics.
A K Pe-
ters/CRC Press, Sept. 2008.
[2] N. Foster and D. Metaxas, ��Realis-
tic animation of liquids,�� in Graph-
ical Models and Image Processing,
pp. 23–30, 1995.
[3] Y. Zhu and R. Bridson, ��Animating
sand as a fluid,�� in ACM SIGGRAPH
2005 Papers, SIGGRAPH ��05, (New
York, NY, USA), pp. 965–972, ACM,
2005.
[4] J. Stam, ��Stable fluids,�� in Pro-
ceedings of the 26th annual confer-
ence on Computer graphics and in-
teractive techniques, SIGGRAPH ��99,
(New York, NY, USA), pp. 121–128,
ACM Press/Addison-Wesley Pub-
lishing Co., 1999.
[5] N. Foster and R. Fedkiw, ��Practi-
cal animation of liquids,�� in Pro-
ceedings of the 28th annual confer-
ence on Computer graphics and in-
teractive techniques, SIGGRAPH ��01,
(New York, NY, USA), pp. 23–30,
ACM, 2001.
[6] M. Evans and F. Harlow, The
particle-in-cell method for hydrody-
namics calculations. LA-2139, 1957.
[7] J. Brackbill and H. Ruppel, ��FLIP:
A method for adaptively zoned,
particle-in-cell calculations of fluid
flows in two dimensions,�� Journal
of Computational Physics, vol. 65,
pp. 314–343, Aug. 1986.
[8] W. E. Lorensen and H. E. Cline,
��Marching cubes: A high resolution
3d surface construction algorithm,��
SIGGRAPH
Comput.
Graph.,
vol. 21, pp. 163–169, August 1987.
[9] T. S. Newman and H. Yi, ��A
survey of the marching cubes al-
13
gorithm,�� Computers & Graphics,
vol. 30, pp. 854–879, Oct. 2006.
[10] J. Yu and G. Turk, ��Reconstructing
surfaces of particle-based fluids
using
anisotropic kernels,��
in
Proceedings of the 2010 ACM SIG-
GRAPH/Eurographics Symposium
on Computer Animation, SCA ��10,
(Aire-la-Ville, Switzerland, Switzer-
land), pp. 217–225, Eurographics
Association, 2010.
A Benchmarks
1
2
3
4
5
Simulation grid
128x64x64 128x64x64 128x64x64 128x64x128 128x64x64
Particle density
23
23
23
23
43
Particles
917 600
917 600
595 200
438 976
402 880
PIC influence
5%
45%
5%
5%
5%
Simulation
34 min
18 min
25 min
8 min
7 min
Reconstruction
82 min
64 min
71 min
155 min
36 min
Rendering
483 min
396 min
310 min
926 min
366 min
Simulation data
2.56 GB
2.56GB
1.66 GB
1.22 GB
1.12 GB
Surface data
820 MB
709 MB
829 MB
1.34 GB
572 MB
Table 2: Benchmark data for all simulation cases
B Images
Figure 5: Dam break (FLIP)
Figure 6: Dam break (PIC)
Figure 7: Dam break with dragon
Figure 8: Falling cube
Figure 9: Splashing dragon