While working out some PixelFlow Demos, i needed some simple (CPU-based for the moment) particle systems which is able to solve collisions and constraints. What started as a basic particle simulation experiment turned out to be quite a powerful tool.

Velocity/Position update

To move a particle, information about its position and velocity is required.

There are several ways to update a particles position, I only mention the 2 most common (probably).

Euler Integration

  v += a * dt
  c += v * dt
  a = 0

Verlet Integration

  v = c - p
  p = c
  c += v + a * 0.5 * dt * dt
  a = 0

v ... velocity
c ... current position
p ... previous position
a ... acceleration
dt ... timestep

I was used to use Euler Integration, but it kind of failed when using it for softbody dynamics. So Verlet Integration it is.

For a more detailed overview i can recommend the following articles



Collision Detection

To solve collisions, information about a particles size (radius) is required. A collision can be handled (in its simplest way), by moving two particles, that are colliding (the distance between them is lower than the sum of their radii) apart by half the collision distance (can be weighted by their masses, etc…).

In a Particle System, one particle most likely has to handle many collisions at once, which usually takes multiple subsequent iterations of the collision update step (relaxation). Obviously this task is very expensive and grows with the number of particles (each particle needs to check all other particles).

For reducing the number of collision checks some space-partitioning techniques can be used, e.g. Quadtree/Octree, BVH, etc.

PixelFlow uses a regular collision Grid at the moment using PPLL (per pixel linked lists) to efficiently store particles in grid-cells. I chose to use the grid because it is very easy to implement, still fast, and can also be implemented on the GPU using GLSL-Shaders (http://thomasdiewald.com/blog/?p=2099).

  dxyz      = pb.c - pa.c
  dmin      = pb.r + pa.r
  dd_cur_sq = squaredLength(dxyz)
  dd_min_sq = dmin * dmin
  if (dd_cur_sq < dd_min_sq) { 
    force = (dd_min_sq / (dd_cur_sq + dd_min_sq) - 0.5f)
    pa.collision -= dxyz * force

Since the particles are stored in a grid and the current position is used to lookup-gridcells, the position cant be updated immediately which is why a temporary “collision” variable is used.



Additionally, all kind of constraints can be added to generate more complex simulations.

For example, a very simple constraint can be: keep two particles apart by maintaining a certain distance (rest-length) between them.

Just this simple rule is all it takes to generate stunning cloth/softbody-simulations. The code is almost identical to the one for handling collisions.

  dxyz      = pb.c - pa.c
  dd_cur_sq = squaredLength(dxyz)
  force     = (dd_rest_sq / (dd_cur_sq + dd_rest_sq) - 0.5f)
  pa.c     -= dxyz * force
  pb.c     += dxyz * force

In contrary to the collision solver, the position can be updated immediately. However, this adds some noticeable error to the result. Since particles previously computed may be changed due to the current results (read/write race), which in turn adds an error to the previous computation. Its noticable in certain situations. Iteratively executing the solver (I use 8 iterations by default) diminishes the error. If more accurary is required, the position-offset must be added/subtracted to a temporary buffer, same as its done by the collision handler.


Cloth Simulation