diff --git a/core/math/MathUtils.h b/core/math/MathUtils.h index 68b31277..85bc2e8a 100644 --- a/core/math/MathUtils.h +++ b/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) { diff --git a/core/quantities/Attractor.cpp b/core/quantities/Attractor.cpp index 3c2446a0..d8286d36 100644 --- a/core/quantities/Attractor.cpp +++ b/core/quantities/Attractor.cpp @@ -52,12 +52,29 @@ void Attractor::interact(IScheduler& scheduler, Storage& storage, const Float dt case ParticleInteractionEnum::ABSORB: { /// \todo parallelize ArrayView r = storage.getValue(QuantityId::POSITION); + ArrayView v = storage.getDt(QuantityId::POSITION); + ArrayView m = storage.getValue(QuantityId::MASS); + ArrayView rho; + if (storage.has(QuantityId::DENSITY)) { + rho = storage.getValue(QuantityId::DENSITY); + } Array 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; }