# Non-real time simulation of overlapping repulsive balls

I want to make a non-real time simulation of overlapping repulsive balls, initially in 2D and later in 3D. With "non-real" I mean that rendering is not needed after each time step.

First of all, consider a simple closed domain for simplicity. In practice, domains will be complex and irregular but always closed.

The balls on the boundaries are fixed and could be overlapping. A fixed ball duplicates itself to produce a free ball of the same size whenever no other ball overlaps it. Both fixed and free balls repel each other but fixed balls cannot move. Note that, duplicant ball should be sufficiently tilted to start repulsion. In elastic colliding balls case, after two balls collide they change direction with some velocity but in this case the balls can stop quickly once they stop overlapping. Free balls move until there is no motion or let's say we solve motion problem until convergence. Then each fixed ball produce a free ball again and this process goes on until no fixed ball can duplicate due to being overlapped by any other ball(s).

I think GPU (CUDA) would be faster to solve this problem but initially I am thinking to write on CPU. However, before proceeding to coding I would like to know "feasibility" of this work. That is, considering a million of balls, approximately how long it would take to simulate this or similar kind of problems in non-real time. For a million of balls, if solution time is in order of minutes, I will dive into the problem. So, my question to experienced programmers is that approximately how long it takes to simulate this kind of problem for a millions of variables? Is it in order of seconds, minutes or hours?

• Can you clarify what the game development application is here? What's the player experience this serves? It looks from the constraints like you might be using it for procedural generation of tile-able blocks with relatively uniformly spaced features - if that's the case, there may be more straightforward means to this end than the physics approach you've described. Or are these particular physics vital to the application you intend? – DMGregory Mar 25 '16 at 5:29
• This is not for game development but asked to this community because I thought game developers have more experience on this kind of problems. The aim is to obtain uniform distribution of points in fast and more importantly parallel manner on GPU. Other techniques (Delaunay, frontal) in my field (mesh generation) are sequential. – Shibli Mar 25 '16 at 9:05
• Gotcha. The reason I asked is because it had the whiff of an X/Y Problem. Generating uniform-ish distributions in parallel is definitely in gamedev's wheelhouse, but as a heads-up we often solve problems like this with a "close enough" hack rather than the provable optimality you may expect from other fields. For us, it usually just has to look plausible, so we use things like randomly displaced grids, approximate poisson disc distributions, or (aperiodic) tiles of pre-generated distributions... Is that okay? – DMGregory Mar 25 '16 at 15:27

It would be on the order of seconds for CPU-based code, and much faster for GPU-based code. I have two sources of experience for this.

The first, writing gravitational simulations (force-based repulsion, exact intersection resolution. As you can see I got much farther with force-based repulsion!). Because this is a sparse problem, for short-ranged collision detection, I use a hash table. Particle positions are mapped to a pair of integers (a,b), and the hash table looks up lists of particles from pairs of integers. There are a million particles in one of the linked videos, and doing the force calculation on a single thread (including the much more expensive gravitational calculations!) took three to four seconds, if I recall correctly.

The second, writing statistical simulations of rigid disk phase transitions. See image below for what I'm talking about. In this case, the disks are densely packed and so a hash table makes no sense. You might as well just create a grid to store the particles in. Every grid would point to an array of pointers to particles (so that if multiple disks lie in a single grid cell, they can be properly handled). Again looking up adjacent disks is easy: just look at the adjacent array indices. The avoidance of mutual intersection of any two disks can be assured much more easily than in what you want to do: You START in a valid configuration where no two disks intersect. Then you choose one disk, and move it randomly in the x,y directions. If the new position intersects any disks, you move it back to where it was. If it doesn't, you leave it. Repeating this millions of times (which can be done very quickly because of the quick collision detection) gives the random images you see below. I would expect single threaded CPU performance for this kind of code for a million particles to be around a second, maybe less (for proposing and checking a new random move for each of a million particles). For reference, this approach is called a "Markov Monte-Carlo approach".

(see Werner Krauth's faculty webpage, his book, or his Coursera course for more. The idea is that as the density of nonintersecting disks increases, there is a phase transition from short-range order [liquid] to infinite-range order [solid] at a critical temperature [rather, critical density].)

To detect collisions, just check if the distance between two adjacent circles is less than the sum of their radii. (if d

To handle collisions between two movable particles, if your particles are at positions x1 and x2, the algorithm I used for the "exact intersection resolution" video above was like this:

Vector x1, x2; //assuming we have a vector class all written.
double d=(x1-x2).lengths();//calculate sqrt((x1.x-x2.x)^2+(x1.y-x2.y)^2) and store it in d.
if(d<2*r){ //if the distance is less than the sum of the radii, there is a collision
Vector direction=(x1-x2).normalized(); //assuming "v.normalized" returns the unit vector in the direction of v
Vector center=(x1+x2)*0.5; //average position of the two circles
x1=center+direction*r; //move it to +r along the axis of collision
x2=center-direction*r; //move it to -r along the axis of collision
}


The problem is that in resolving the collision between TWO particles, we may have introduced collisions of either one of the particles with some third particle. You can see this problem in the "exact intersection resolution" video! That is why I prefer force-based methods.

If your goal is just to generate pictures/distributions like Krauth's above, use the Monte Carlo approach. If you MUST use the algorithm in your OP, I would really really suggest a physics based approach (ie, storing a position and velocity and applying a strong spring repulsion force, then updating for many timesteps dt). Going from a valid configuration of particles to a valid configuration of particles is easy (Monte Carlo). Going from an invalid configuration of intersecting particles to a valid configuration is hard, impossible in some cases, and is best solved by using a physics force-based approach.

• in the first link it is mentioned that render time was 8 days not 3 to 4 seconds. – Shibli Jan 16 '17 at 19:49
• @Shibli 3 to 4 seconds per timestep, and this included the much more computationally intensive long-range 1/r^2 forces. – user10968 Jan 16 '17 at 19:58