Localized Attraction-Repulsion (LAR)

Published on Jul 29, 2020

In this article I will describe a particle simulation that can produce different types of matter based on different input rules. This simulation is based on physical properties and physical rules but the rules were altered in a way that allows for the emergence of intriguing natural-like patterns. The implementation of part 1 and 2 will be done in Processing and the source code will be available at the end. The following illustrations demonstrates a few possible usages for this simulation.

Attraction Repulsion between particles

This simulation takes place in a 2d space, although it could also be done in 3D. It is based on attraction / repulsion between a large amount of particles. Each particle of our simulation will be attracted / repelled by the others. We will first cover the mathematics to simulate such phenomena by taking 2 particles $P\_1$ and $P\_2$ .

Both particles have a position in the 2D space, respectively $(x\_{p1}; y\_{p1})$ and $(x\_{p2}; y\_{p2})$ . The following illustration shows the attraction and repulsion forces generated by $P\_2$ over $P\_1$ .

Figure 1. (in green) the attraction, (in red) the repulsion

We will be using Coulomb’s law [1], an experimental law of physics that quantifies the amount of force between two stationary, electrically charged particles of charge $q\_1$ and $q\_2$ . The value of $q\_1$ and $q\_2$ doesn’t really matter in our case, you will see why later. The law of Coulombs states that the force of attraction $\overrightarrow{F\_a}$ between $P\_2$ and $P\_1$ is:

F\_a = k \cdot \frac{q\_1 \cdot q\_2}{r^2}

where $k$ is the Coulomb’s constant
and $r$ is the distance between the particles

In our case, we want to control the attraction strength between the particules, and it would be great to do so with a single value. We will consider that all of our particles have the same electric charge $q$ . Thus the formula can be simplified:

F\_a = k \cdot \frac{q^2}{r^2}

Because $k$ and $q$ are both constants, we can define them with a single value $S\_a$ , our the strength of the force of attraction between the particles:

S\_a = k \cdot q^2

Therefore our formula can be simplified to the following:

F\_a = \frac{S\_a}{r^2}

The exact same concept can be applied with the force of repulsion, let $S\_r$ be the repulsion strength.

To compute the actual vector of the force of attraction, the operation is simple. We can first compute the direction $\overrightarrow{D}$ between $P\_1$ and $P\_2$ :

\overrightarrow{D} = \begin{bmatrix}x\_{p2} - x\_{p1} \\ y\_{p2} - y\_{p1}\end{bmatrix}

The scalar distance between $P\_1$ and $P\_2$ will be the magnitude of $\overrightarrow{D}$ :

r = ||D||

Computing $\overrightarrow{F\_a}$ and $\overrightarrow{F\_r}$ is now a trivial task. It can be done by taking the normalized vector of the direction, and multiply it by the scalar value of the force:

\overrightarrow{F\_a} = \frac{\overrightarrow{D}}{||D||} \cdot |F\_a|

\overrightarrow{F\_r} = - \frac{\overrightarrow{D}}{||D||} \cdot |F\_r|

Which, in code, translates into:

// compute the force of p2 over p1
PVector D = p2.position.copy().sub(p1.position);
float r = D.mag();
D = D.normalize();
PVector Fa = D.copy().mult(attrStrength).div(r*r);
PVector Fr = D.copy().mult(-repStrength).div(r*r);

Figure 2. Graphical output of the code below (vectors were scaled up so that we can see them)

Now that we know how to compute these forces, we can already apply this model to our 2 particles. To update a particle position, we can, on each iteration, compute the sum of the forces applied to it. This sum will be the acceleration. The acceleration can then be added to the velocity, and the velocity to the position. At the end of the iteration, the particle acceleration is set back to 0, and to simulate the friction we multiply the velocity by a scalar factor, the friction decay (0.98 for instance). I already covered this in a previous article, if you are struggling with this step you can check the source code at the end of the article.

For now, we will work in pixel coordinates. Which means that we should take large values for the attraction and repulsion strength for any motion to happen: $S\_a = 180$ and $S\_r = 200$ . It is important to note that keeping a repulsion strength higher than the attraction strength will ensure that the particles keep a correct distanciation between them.

Figure 3. Attraction / Repulsion between 2 particles

Time to throw some more particles in !

Figure 4. With more particles, same rules applied

We already observe an interesting effect, the particles, over time, tend to get distributed evenly in the space. Before jumping to the next part, let’s introduce a last component to our system. The particles end up being to far away from the center over time. A nice solution is to add some attraction to the center of the world. This is a possible implementation, where the attraction to the center gets higher the further we get from such center:

// we add attraction to the center
p1.acceleration.add(p1.position.copy().sub(center).mult(-0.01));

Let $C\_a$ be the attraction strength towards the center. This is the outcome of this system with the followign values:

$S\_a$	$S\_r$	$C\_a$
110	130	0.01

Figure 5. Addition of some attraction to the center

The particles seems to end up in a sort of vortex over time. This is an interesting behavior, already great ideas could be built from there. This behavior is logical though, think about how gravitation works 🤔

By selecting values carefully for our parameters, we can already get a system that stabilizes himelf:

$S\_a$	$S\_r$	$C\_a$
110	130	0.001

Figure 6. The system stabilizes with the correct set of parameters

Figure 7. 500 particles, same parameters

That’s it, these are the very basis for an attraction repulsion system. In the following part, I will propose an idea that allows for more interesting behaviors to emerge.

Attraction range and repulsion range

Currently, every particle is affected by every other particle. This is a realistic behavior, however we are not bounded by physics here. We are going to introduce 2 new parameters: the attraction range $R\_a$ and the repulsion range $R\_r$ . Particles will only be attracted / repelled to each other if the scalar distance between them is under those threshold. Implementation is trivial.

$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$
110	130	0.001	100	90

Figure 8. Birth of local vortex

😱 we can already see an interesting behavior with that vortex in the middle only. But let’s explore with some different values…

$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$
20	130	0.001	150	40

Figure 9. A better vortex ?

😱😱 that is even a better vortex than the previous one ! Let’s dig more:

$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$
20	18	0.0001	40	80

Figure 10. I wish there was a better vortex

Hum… there seem to be an issue with those values. Indeed, when particles gets attracted too close, the force generated becomes so high that they end up being sent to the moon. We can fix this issue by adding some collisions to our system.

I will implement a really dumb collision model; it does the job but by no mean is physically correct:

Which translates into the following code:

if (r < p1.radius + p2.radius) {
  PVector mv = D.copy().mult(-((p1.radius + p2.radius) - r));
  p1.velocity.add(mv);
  p2.velocity.add(mv.mult(-0.2));
}

I added a radius to the particles that is only used for the collisions and for drawing the particles. It was set to 10 in the following illustration.

Figure 12. Introduction of our beautiful collision model

😱 BOOM ! Up to our state-of-the-art collision model, the vortex is back !

Please note that 2 more parameters can be added to the list-of-things-we-can-tweak to get funny results: the collision strength $Col$ (a factor that indicates by how much a particle should be moved out of another one it overlaps) and the collision response $Col\_r$ (a factor that indicates how much energy is transferred upon collision, it was set to 0.2 in the code above). In the next illustration, I lowered the radius of the particles to 5.

$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
20	18	0.0001	40	80	0.5	0.1

Figure 13. Emergence of localized behaviors

Let’s try to emphasize this localized behavior. One thing that could be done is to increase the friction. The decay of friction $d\_f$ is currently set to 0.98. By setting it to a smaller value, velocity will decrease faster over time (resulting in increased friction). I will also decrease the maximum velocity $max\_s$ from 3 to 1.

$max\_s$	$d\_f$	$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
1	0.92	25	18	0.0001	40	80	0.5	0

Figure 14. Particles group into clusters

With this particular set of settings, we can observe the emergence of clusters of particles. Beautiful.

$max\_s$	$d\_f$	$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
1	0.92	20	18	0.0001	18	13	0.5	0

Figure 15. A big cluster of particles

With these settings we get a big cluster of particles. As a quick note, when I work on a new algorithm, this is usually the step where I look back and try to see if pushing it further could be interesting. In that particular case, I think the wide variety of results and the variety of parameters are both indicators that better results could come down the line. So let’s dig more.

The next illustration is where this algorithm totally blew my mind. Let’s try to find a configuration where particles tend to stick together. For that to happen, I think increasing the friction while keeping a great attraction at a lower than the repulsion could to the trick.

$max\_s$	$d\_f$	$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
1	0.8	20	30	0.0001	13	20	0.5	0

Figure 16. Reaction-Diffusion patterns

I was shocked when I saw this. In the middle, some reaction diffusion patterns seem to emerge, however they end up being lost due to the overall repulsion being too high. By increasing the attraction strength just above the repulsion strength, we can ensure that the particles with keep bonding.

$max\_s$	$d\_f$	$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
1	0.8	32	30	0.0001	13	20	0.5	0

Figure 17. Consistent Reaction-Diffusion patterns

Let’s try a slower growth:

$max\_s$	$d\_f$	$S\_a$	$S\_r$	$C\_a$	$R\_a$	$R\_r$	$Col$	$Col\_r$
1	0.2	32	30	0.0001	13	20	0.5	0

Figure 18. A slower growth of a RD pattern

We have some weirdos running around: overlapping particules at initialization are strongly repelled. One solution could be to only compute the forces between 2 particles if their distance is higher than soem threshold.

And let’s throw some more particles in, 800. I also increased to friction decay to 0.95 to get some motion back.

Figure 19. With 800 particles

At this point, the next step on the line is to scale this system up. 800 particles is fun, but going to higher scales would destroy the CPU. Various optimizations could be made, such as space partitioning algorithms. OR the computations could be moved to the GPU, by using compute shaders.

The source code of the implementation on Processing, up to this point, is available on Github:

https://github.com/bcrespy/LocalizedAttractionRepulsion

Scaling this system up

Compute Shaders [2] have the ability to store small chunks of data into the GPU memory, that can be shared across the different cores of the Graphics Card. I won’t go into details in this article because it is out of the scope, you should know that usually, when it comes to particles interacting with each other, Compute Shaders can increase performances by a lot.

The following illustrations were generated using the same algorithm presented above. The only difference is the scale of the space in which the particles evolve, it is not in pixel coordinates but is within [-2; 2]. This is why the overlay values are much smaller than the ones presented before.

The particles are colored based on the magnitude of their velocity from white (no velocity) to red (high velocity). The number of particles is 12 544, and each particle interacts with each other particle.

Figure 20. RD patterns with clusters of particles

Figure 21. Small patterns followed by small clusters

Figure 22. Clean RD patterns with rotating clusters

Figure 23. A weird kind of matter

Because this simulation relies on particles, we can interact with the particles easily. In the following illustrations, a vector field is generated from a point that moves in the space. Basically, particles are repelled by this point. I used it to show you that this simulation mimics the behavior of liquids, which is pretty neat.

Figure 24. Interacting with fast reconstructing RD patterns

Figure 25. Stricky fluid behavior

I used this technique when I did some visuals for Inkie. I layered many other techniques to get this particular look / interactions with the lasers, but the idea remains the same:

I have also explored this system furthermore on my instagram:

Follow me on instagram for more content

That will be it for this article. Thanks for reading through, I hope you will find some inspiration in that technique and build your own stuff with it 👽

Link to the repo of the Processing project:

https://github.com/bcrespy/LocalizedAttractionRepulsion

References

Coulomb’s law, Wikipédia
Compute Shaders, Khronos