Collision.h

Go to the documentation of this file.

#pragma once
 
 
#include "objects/containers/FlatSet.h"
#include "physics/Functions.h"
#include "quantities/Storage.h"
#include "system/Settings.h"
 
NAMESPACE_SPH_BEGIN
 
enum class CollisionResult {
    NONE,          
    BOUNCE,        
    FRAGMENTATION, 
    MERGER,        
    EVAPORATION,   
};
 
class ICollisionHandler : public Polymorphic {
public:
    virtual void initialize(Storage& storage) = 0;
 
    virtual CollisionResult collide(const Size i, const Size j, FlatSet<Size>& toRemove) = 0;
};
 
class IOverlapHandler : public Polymorphic {
public:
    virtual void initialize(Storage& storage) = 0;
 
    virtual bool overlaps(const Size i, const Size j) const = 0;
 
    virtual void handle(const Size i, const Size j, FlatSet<Size>& toRemove) = 0;
};
 
template <typename T>
Tuple<T, T> minMax(const T& t1, const T& t2) {
    if (t1 < t2) {
        return { t1, t2 };
    } else {
        return { t2, t1 };
    }
}
 
template <typename T>
INLINE T weightedAverage(const T& v1, const Float w1, const T& v2, const Float w2) {
    SPH_ASSERT(w1 + w2 > 0._f, w1, w2);
    return (v1 * w1 + v2 * w2) / (w1 + w2);
}
 
INLINE bool areParticlesBound(const Float m_sum, const Float h_sum, const Vector& dv, const Float limit) {
    const Float vEscSqr = 2._f * Constants::gravity * m_sum / h_sum;
    const Float vRelSqr = getSqrLength(dv);
    return vRelSqr * limit < vEscSqr;
}
 
// ----------------------------------------------------------------------------------------------------------
// Collision Handlers
// ----------------------------------------------------------------------------------------------------------
 
class NullCollisionHandler : public ICollisionHandler {
public:
    virtual void initialize(Storage& UNUSED(storage)) override {}
 
    virtual CollisionResult collide(const Size UNUSED(i),
        const Size UNUSED(j),
        FlatSet<Size>& UNUSED(toRemove)) override {
        return CollisionResult::NONE;
    }
};
 
class MergingCollisionHandler : public ICollisionHandler {
private:
    ArrayView<Vector> r, v;
    ArrayView<Float> m;
    ArrayView<Vector> L;
    ArrayView<Vector> omega;
    ArrayView<SymmetricTensor> I;
    ArrayView<Tensor> E;
 
 
    Float bounceLimit;
    Float rotationLimit;
 
public:
    explicit MergingCollisionHandler(const RunSettings& settings) {
        bounceLimit = settings.get<Float>(RunSettingsId::COLLISION_BOUNCE_MERGE_LIMIT);
        rotationLimit = settings.get<Float>(RunSettingsId::COLLISION_ROTATION_MERGE_LIMIT);
    }
 
    explicit MergingCollisionHandler(const Float bounceLimit, const Float rotationLimit)
        : bounceLimit(bounceLimit)
        , rotationLimit(rotationLimit) {}
 
    virtual void initialize(Storage& storage) override {
        ArrayView<Vector> dv;
        tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
        m = storage.getValue<Float>(QuantityId::MASS);
        omega = storage.getValue<Vector>(QuantityId::ANGULAR_FREQUENCY);
 
        if (storage.has(QuantityId::MOMENT_OF_INERTIA)) {
            I = storage.getValue<SymmetricTensor>(QuantityId::MOMENT_OF_INERTIA);
            E = storage.getValue<Tensor>(QuantityId::LOCAL_FRAME);
            L = storage.getValue<Vector>(QuantityId::ANGULAR_MOMENTUM);
        }
    }
 
    virtual CollisionResult collide(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        // set radius of the merger so that the volume is conserved
        const Float h_merger = root<3>(pow<3>(r[i][H]) + pow<3>(r[j][H]));
 
        // conserve total mass
        const Float m_merger = m[i] + m[j];
 
        // merge so that the center of mass remains unchanged
        const Vector r_merger = weightedAverage(r[i], m[i], r[j], m[j]);
 
        // converve momentum
        const Vector v_merger = weightedAverage(v[i], m[i], v[j], m[j]);
 
        Vector omega_merger;
        SymmetricTensor I_merger;
        Vector L_merger;
        Tensor E_merger = Tensor::identity();
 
        // Never modify particle values below UNTIL we know the collision is accepted; save the preliminary
        // values to _merger variables!
 
        if (I) {
            // So far this is just an experimental branch, not intended to be used for "serious" simulations.
            // This assert is just a check that it does not get enabled by accident.
            SPH_ASSERT(Assert::isTest);
 
            // compute inertia tensors in inertial frame
            const SymmetricTensor I1 = transform(I[i], convert<AffineMatrix>(E[i]));
            const SymmetricTensor I2 = transform(I[j], convert<AffineMatrix>(E[j]));
 
            // sum up the inertia tensors, but first move them to new origin
            I_merger = Rigid::parallelAxisTheorem(I1, m[i], r_merger - r[i]) +
                       Rigid::parallelAxisTheorem(I2, m[j], r_merger - r[j]);
 
            // compute the total angular momentum - has to be conserved
            L_merger = m[i] * cross(r[i] - r_merger, v[i] - v_merger) + //
                       m[j] * cross(r[j] - r_merger, v[j] - v_merger) + //
                       L[i] + L[j];
            // L = I*omega  =>  omega = I^-1 * L
            omega_merger = I_merger.inverse() * L_merger;
 
            // compute the new local frame of the merger and inertia tensor in this frame
            Eigen eigen = eigenDecomposition(I_merger);
            I_merger = SymmetricTensor(eigen.values, Vector(0._f));
            SPH_ASSERT(isReal(I_merger));
 
            E_merger = convert<Tensor>(eigen.vectors);
            SPH_ASSERT(isReal(E_merger));
 
            SPH_ASSERT(isReal(L_merger) && isReal(omega_merger), L_merger, omega_merger);
            SPH_ASSERT(almostEqual(getSqrLength(E_merger.row(0)), 1._f, 1.e-6_f));
            SPH_ASSERT(almostEqual(getSqrLength(E_merger.row(1)), 1._f, 1.e-6_f));
            SPH_ASSERT(almostEqual(getSqrLength(E_merger.row(2)), 1._f, 1.e-6_f));
 
        } else {
            L_merger = m[i] * cross(r[i] - r_merger, v[i] - v_merger) + //
                       m[j] * cross(r[j] - r_merger, v[j] - v_merger) + //
                       Rigid::sphereInertia(m[i], r[i][H]) * omega[i] + //
                       Rigid::sphereInertia(m[j], r[j][H]) * omega[j];
            omega_merger = Rigid::sphereInertia(m_merger, h_merger).inverse() * L_merger;
        }
 
        if (!this->acceptMerge(i, j, h_merger, omega_merger)) {
            return CollisionResult::NONE;
        }
 
        // NOW we can save the values
 
        m[i] = m_merger;
        r[i] = r_merger;
        r[i][H] = h_merger;
        v[i] = v_merger;
        v[i][H] = 0._f;
        omega[i] = omega_merger;
 
        if (I) {
            I[i] = I_merger;
            L[i] = L_merger;
            E[i] = E_merger;
        }
 
        SPH_ASSERT(isReal(v[i]) && isReal(r[i]));
        toRemove.insert(j);
        return CollisionResult::MERGER;
    }
 
private:
    INLINE bool acceptMerge(const Size i, const Size j, const Float h, const Vector& omega) const {
        if (!areParticlesBound(m[i] + m[j], r[i][H] + r[j][H], v[i] - v[j], bounceLimit)) {
            // moving too fast, reject the merge
            return false;
        }
        const Float omegaCritSqr = Constants::gravity * (m[i] + m[j]) / pow<3>(h);
        const Float omegaSqr = getSqrLength(omega);
        if (omegaSqr * rotationLimit > omegaCritSqr) {
            // rotates too fast, reject the merge
            return false;
        }
        return true;
    }
};
 
class ElasticBounceHandler : public ICollisionHandler {
protected:
    ArrayView<Vector> r, v;
    ArrayView<Float> m;
 
    struct {
        Float n;
 
        Float t;
 
    } restitution;
 
public:
    explicit ElasticBounceHandler(const RunSettings& settings) {
        restitution.n = settings.get<Float>(RunSettingsId::COLLISION_RESTITUTION_NORMAL);
        restitution.t = settings.get<Float>(RunSettingsId::COLLISION_RESTITUTION_TANGENT);
    }
 
    ElasticBounceHandler(const Float n, const Float t) {
        restitution.n = n;
        restitution.t = t;
    }
 
    virtual void initialize(Storage& storage) override {
        ArrayView<Vector> dv;
        tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
        m = storage.getValue<Float>(QuantityId::MASS);
    }
 
    virtual CollisionResult collide(const Size i, const Size j, FlatSet<Size>& UNUSED(toRemove)) override {
        const Vector dr = getNormalized(r[i] - r[j]);
        const Vector v_com = weightedAverage(v[i], m[i], v[j], m[j]);
        v[i] = this->reflect(v[i], v_com, -dr);
        v[j] = this->reflect(v[j], v_com, dr);
 
        // no change of radius
        v[i][H] = 0._f;
        v[j][H] = 0._f;
 
        SPH_ASSERT(isReal(v[i]) && isReal(v[j]));
        return CollisionResult::BOUNCE;
    }
 
private:
    INLINE Vector reflect(const Vector& v, const Vector& v_com, const Vector& dir) {
        SPH_ASSERT(almostEqual(getSqrLength(dir), 1._f), dir);
        const Vector v_rel = v - v_com;
        const Float proj = dot(v_rel, dir);
        const Vector v_t = v_rel - proj * dir;
        const Vector v_n = proj * dir;
 
        // flip the orientation of normal component (bounce) and apply coefficients of restitution
        return restitution.t * v_t - restitution.n * v_n + v_com;
    }
};
 
template <typename TPrimary, typename TFallback>
class FallbackHandler : public ICollisionHandler {
private:
    TPrimary primary;
    TFallback fallback;
 
public:
    FallbackHandler(const RunSettings& settings)
        : primary(settings)
        , fallback(settings) {}
 
    virtual void initialize(Storage& storage) override {
        primary.initialize(storage);
        fallback.initialize(storage);
    }
 
    virtual CollisionResult collide(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        CollisionResult result = primary.collide(i, j, toRemove);
        if (result == CollisionResult::NONE) {
            return fallback.collide(i, j, toRemove);
        } else {
            return result;
        }
    }
};
 
class FragmentationHandler : public ICollisionHandler {
public:
    // ParametricRelations, directionality of fragments
 
    virtual CollisionResult collide(const Size i, const Size j, FlatSet<Size>& UNUSED(toRemove)) override {
        (void)i;
        (void)j;
        NOT_IMPLEMENTED;
        return CollisionResult::FRAGMENTATION;
    }
};
 
 
// ----------------------------------------------------------------------------------------------------------
// Overlap Handlers
// ----------------------------------------------------------------------------------------------------------
 
class NullOverlapHandler : public IOverlapHandler {
public:
    virtual void initialize(Storage& UNUSED(storage)) override {}
 
    virtual bool overlaps(const Size UNUSED(i), const Size UNUSED(j)) const override {
        return false;
    }
 
    virtual void handle(const Size UNUSED(i),
        const Size UNUSED(j),
        FlatSet<Size>& UNUSED(toRemove)) override {}
};
 
class MergeOverlapHandler : public IOverlapHandler {
private:
    MergingCollisionHandler handler;
 
public:
    MergeOverlapHandler()
        : handler(0._f, 0._f) {}
 
    virtual void initialize(Storage& storage) override {
        handler.initialize(storage);
    }
 
    virtual bool overlaps(const Size UNUSED(i), const Size UNUSED(j)) const override {
        return true;
    }
 
    virtual void handle(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        handler.collide(i, j, toRemove);
    }
};
 
template <typename TFollowupHandler>
class RepelHandler : public IOverlapHandler {
private:
    TFollowupHandler handler;
 
    ArrayView<Vector> r, v;
    ArrayView<Float> m;
 
public:
    explicit RepelHandler(const RunSettings& settings)
        : handler(settings) {}
 
    virtual void initialize(Storage& storage) override {
        ArrayView<Vector> dv;
        tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
        m = storage.getValue<Float>(QuantityId::MASS);
 
        handler.initialize(storage);
    }
 
    virtual bool overlaps(const Size UNUSED(i), const Size UNUSED(j)) const override {
        // this function is called only if the spheres intersect, which is the only condition for this handler
        return true;
    }
 
    virtual void handle(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        Vector dir;
        Float dist;
        tieToTuple(dir, dist) = getNormalizedWithLength(r[i] - r[j]);
        dir[H] = 0._f; // don't mess up radii
        // can be only used for overlapping particles
        SPH_ASSERT(dist < r[i][H] + r[j][H], dist, r[i][H] + r[i][H]);
        const Float x1 = (r[i][H] + r[j][H] - dist) / (1._f + m[i] / m[j]);
        const Float x2 = m[i] / m[j] * x1;
        r[i] += dir * x1;
        r[j] -= dir * x2;
        SPH_ASSERT(almostEqual(getSqrLength(r[i] - r[j]), sqr(r[i][H] + r[j][H])),
            getSqrLength(r[i] - r[j]),
            sqr(r[i][H] + r[j][H]));
 
        SPH_ASSERT(isReal(v[i]) && isReal(v[j]));
        SPH_ASSERT(isReal(r[i]) && isReal(r[j]));
 
        // Now when the two particles are touching, handle the collision using the followup handler.
        handler.collide(i, j, toRemove);
    }
};
 
class InternalBounceHandler : public IOverlapHandler {
private:
    ElasticBounceHandler handler;
 
    ArrayView<Vector> r, v;
 
public:
    explicit InternalBounceHandler(const RunSettings& settings)
        : handler(settings) {
        // this handler allows overlaps of particles, so it should never be used with point particles, as we
        // could potentially get infinite accelerations
        SPH_ASSERT(settings.get<GravityKernelEnum>(RunSettingsId::GRAVITY_KERNEL) !=
                   GravityKernelEnum::POINT_PARTICLES);
    }
 
    virtual void initialize(Storage& storage) override {
        ArrayView<Vector> dv;
        tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
 
        handler.initialize(storage);
    }
 
    virtual bool overlaps(const Size i, const Size j) const override {
        // overlap needs to be handled only if the particles are moving towards each other
        const Vector dr = r[i] - r[j];
        const Vector dv = v[i] - v[j];
        return dot(dr, dv) < 0._f;
    }
 
    virtual void handle(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        handler.collide(i, j, toRemove);
    }
};
 
class MergeBoundHandler : public IOverlapHandler {
private:
    MergingCollisionHandler handler;
 
public:
    ArrayView<Vector> r, v, omega;
    ArrayView<Float> m;
 
    Float bounceLimit;
    Float rotationLimit;
 
public:
    explicit MergeBoundHandler(const RunSettings& settings)
        : handler(settings) {
        bounceLimit = settings.get<Float>(RunSettingsId::COLLISION_BOUNCE_MERGE_LIMIT);
        rotationLimit = settings.get<Float>(RunSettingsId::COLLISION_ROTATION_MERGE_LIMIT);
    }
 
    virtual void initialize(Storage& storage) override {
        ArrayView<Vector> dv;
        tie(r, v, dv) = storage.getAll<Vector>(QuantityId::POSITION);
        omega = storage.getValue<Vector>(QuantityId::ANGULAR_FREQUENCY);
        m = storage.getValue<Float>(QuantityId::MASS);
 
        handler.initialize(storage);
    }
 
    virtual bool overlaps(const Size i, const Size j) const override {
        const Float m_merger = m[i] + m[j];
        if (!areParticlesBound(m_merger, r[i][H] + r[j][H], v[i] - v[j], bounceLimit)) {
            // moving too fast, reject the merge
            return false;
        }
 
        const Float h_merger = root<3>(pow<3>(r[i][H]) + pow<3>(r[j][H]));
        const Float omegaCritSqr = Constants::gravity * (m[i] + m[j]) / pow<3>(h_merger);
 
        const Vector r_merger = weightedAverage(r[i], m[i], r[j], m[j]);
        const Vector v_merger = weightedAverage(v[i], m[i], v[j], m[j]);
        const Vector L_merger = m[i] * cross(r[i] - r_merger, v[i] - v_merger) + //
                                m[j] * cross(r[j] - r_merger, v[j] - v_merger) + //
                                Rigid::sphereInertia(m[i], r[i][H]) * omega[i] + //
                                Rigid::sphereInertia(m[j], r[j][H]) * omega[j];
        const Vector omega_merger = Rigid::sphereInertia(m_merger, h_merger).inverse() * L_merger;
        const Float omegaSqr = getSqrLength(omega_merger);
        if (omegaSqr * rotationLimit > omegaCritSqr) {
            // rotates too fast, reject the merge
            return false;
        }
 
        return true;
    }
 
    virtual void handle(const Size i, const Size j, FlatSet<Size>& toRemove) override {
        handler.collide(i, j, toRemove);
    }
};
 
NAMESPACE_SPH_END