LAR (Localized Attraction-Repulsion): generative matter

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.

1. 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
1101300.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
1101300.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.

2. 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
1101300.00110090
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
201300.00115040
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
20180.00014080
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:

Figure 11. A hacky collision model

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));
}

Yeah also 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
20180.000140800.50.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
10.9225180.000140800.50
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
10.9220180.000118130.50
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
10.820300.000113200.50
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
10.832300.000113200.50
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
10.232300.000113200.50
Figure 18. A slower growth of a RD pattern

We have some weirdos running around. Honestly I don't know where they come from, I will just ignore them because those patterns are too great. Well as a matter a fact I kinda know where they come from: particules which were randomly placed on top of each other at initialization. One solution could be to only compute the forces between 2 particles if their distance is higher than 0.5.

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. There are optimizations that could be made, such as storing the particles in a grid and only compare particles from adjacent cells. And other less naive approaches. 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

3. 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 Shader 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:

References

  1. Coulomb's law, Wikipédia
  2. Compute Shaders, Khronos

Leave a Reply

Your email address will not be published. Required fields are marked *