conserving mass and momentum when absorbing particles

core/math/MathUtils.h

@@ -399,6 +399,11 @@ INLINE constexpr Float sphereSurfaceArea(const Float radius) {
     return 4._f * PI * pow<2>(radius);
 }
 
+/// Inverse to \ref sphereVolume.
+INLINE Float volumeEquivalentRadius(const Float volume) {
+    return cbrt(3 * volume / (4 * PI));
+}
+
 /// Checks if two values are equal to some given accuracy.
 /// \note We use <= rather than < on purpose as EPS for integral types is zero.
 INLINE bool almostEqual(const Float& f1, const Float& f2, const Float& eps = EPS) {

core/quantities/Attractor.cpp

@@ -52,12 +52,29 @@ void Attractor::interact(IScheduler& scheduler, Storage& storage, const Float dt
     case ParticleInteractionEnum::ABSORB: {
         /// \todo parallelize
         ArrayView<const Vector> r = storage.getValue<Vector>(QuantityId::POSITION);
+        ArrayView<const Vector> v = storage.getDt<Vector>(QuantityId::POSITION);
+        ArrayView<const Float> m = storage.getValue<Float>(QuantityId::MASS);
+        ArrayView<const Float> rho;
+        if (storage.has(QuantityId::DENSITY)) {
+            rho = storage.getValue<Float>(QuantityId::DENSITY);
+        }
         Array<Size> toRemove;
+        Float absorbed_mass = 0._f;
+        Float absorbed_volume = 0._f;
+        Vector absorbed_momentum = Vector(0);
         for (Size i = 0; i < r.size(); ++i) {
             if (getSqrLength(position - r[i]) < sqr(radius)) {
                 toRemove.push(i);
+                absorbed_mass += m[i];
+                absorbed_volume += !rho.empty() ? m[i] / rho[i] : 0;
+                absorbed_momentum += m[i] * (v[i] - velocity);
             }
         }
+        if (absorbed_mass > 0) {
+            velocity += absorbed_momentum / mass;
+            mass += absorbed_mass;
+            radius = volumeEquivalentRadius(sphereVolume(radius) + absorbed_volume);
+        }
         storage.remove(toRemove, Storage::IndicesFlag::INDICES_SORTED | Storage::IndicesFlag::PROPAGATE);
         break;
     }