Compute on all points in geometry nodes

Question

I try (for a while now) to make a gravitation simulation with geometry nodes. The computation is not complex, and I get something that is almost OK.

The problem is that it seems that it's not possible to make a computation from the entire points. We only can get the "nearest" point position.

So, the simulation works only partially, the gravitation is applied for all points with its neighbor, but it should be computed with all the others points.

I provide the example blend file.

https://drive.google.com/file/d/19iRfLL1_WBIOQFTdYoSAAh9RkzFrNNlT/view?usp=sharing

Answer 1

(Using Blender 3.6.8)

Let $N$ be the number of bodies. Gravity acceleration for body $i$ due to bodies $j$ is a function of $$ \forall i, 1 \le i \le N, \ \vec{a_i} = \sum_{j \neq i}^N \frac{\vec{M_i M_j}}{\|\vec{M_i M_j}\|^3} \simeq \sum_{j}^N \frac{\vec{M_i M_j}}{\|\vec{M_i M_j}\|^3 + \varepsilon}$$ where $M_k$ is the position of body $k$ and $\varepsilon$ a small positive not zero value.
As a conclusion, two loops over $i$ and $j$ are required.

---

The following approach is using a Simulation Zone for the $i$ loop, and an Attribute Statistic node for the $j$ loop. From Blender 4, a Repeat Zone should be preferred over the Simulation Zone.

---

The following approach is using a Cartesian grid to mimic the $\left[\vec{M_i M_j}\right]_{N \times N}$ matrix. The doublet $(i,j)$ is encoded in the grid points index $n$ as $n=i+j \times N$. Furthermore, $i$ specifies those points ID used by the Accumulate Field node, which is then computing the summation over $j$ for each $i$. The double-loop over $(i,j)$ is performed through the Set Position node.

---

Resources:

Answer 2

(if you're using this as a template for your gravity sim, notice I don't use a gravitational constant)

Repeat Zones are a bit slow, not sure if the Stef's answer would be slower than mine, worth to check… Meanwhile here's an answer using the "quadratic geometry explosion" (just made that term up) pattern, that I also recently used here. And, like there, I'm optimizing distance checks by applying negative $-100 \%$ scale, so that the point I compare to is at origin, and position vector length equals distance to that point. Not sure if this is actually a useful optimization here, though. I also capture radius on the parent-point so it transfers to a spawned instance, and upon realizing instances transfers to the children-points. Now I can easily access the radius of a current point, and a parent point (captured). I use radii to store point masses (could use a separate attribute, but that would increase cache size)

The Attribute Statistic node returns a single value, but an Accumulate Field node also sums up values and can be separated by group. The positioning on index to sample nearest at that index could be optimized away as well as the last Capture Attribute node could be removed…

Notice this relies on the input geometry. I prepared two examples of input geometry for testing:

Test #1: Two-body problem

Test #2: Unstable Constellation

This is why triple-star systems don't exist: there's no stabilizing mechanism that puts bodies back at equal distances when two of them get closer than the other due to some external chaos (here the chaos is just floating point datatype inaccuracy).