Functions.h

Go to the documentation of this file.

#pragma once
 
 
#include "math/rng/Rng.h"
#include "objects/geometry/SymmetricTensor.h"
#include "objects/wrappers/Flags.h"
#include "physics/Constants.h"
 
NAMESPACE_SPH_BEGIN
 
namespace Analytic {
 
class StaticSphere {
private:
    Float r0;
 
    Float rho;
 
public:
    StaticSphere(const Float r0, const Float rho)
        : r0(r0)
        , rho(rho) {}
 
    INLINE Float getPressure(const Float r) const {
        if (r > r0) {
            return 0._f;
        }
        return 2._f / 3._f * PI * Constants::gravity * sqr(rho) * (sqr(r0) - sqr(r));
    }
 
    INLINE Vector getAcceleration(const Vector& r) const {
        const Float l = getLength(r);
        const Float l0 = min(r0, l);
        return -Constants::gravity * rho * sphereVolume(l0) * r / pow<3>(l);
    }
 
    INLINE Float getEnergy() const {
        return -16._f / 15._f * sqr(PI) * sqr(rho) * Constants::gravity * pow<5>(r0);
    }
};
 
} // namespace Analytic
 
namespace Rigid {
 
INLINE SymmetricTensor sphereInertia(const Float m, const Float r) {
    return SymmetricTensor::identity() * (0.4_f * m * sqr(r));
}
 
INLINE SymmetricTensor parallelAxisTheorem(const SymmetricTensor& I, const Float m, const Vector& a) {
    return I + m * (SymmetricTensor::identity() * getSqrLength(a) - symmetricOuter(a, a));
}
 
} // namespace Rigid
 
INLINE Float getCriticalFrequency(const Float rho) {
    return sqrt(sphereVolume(1._f) * rho * Constants::gravity);
}
 
Float evalBenzAsphaugScalingLaw(const Float D, const Float rho);
 
Float getImpactEnergy(const Float R, const Float r, const Float v);
 
Float getEffectiveImpactArea(const Float R, const Float r, const Float phi);
 
Float getImpactorRadius(const Float R_pb, const Float v_imp, const Float QOverQ_D, const Float rho);
 
Float getImpactorRadius(const Float R_pb,
    const Float v_imp,
    const Float phi,
    const Float Q_effOverQ_D,
    const Float rho);
 
 
class ImpactCone {
private:
    AffineMatrix frame;
    Float v_c;
 
    UniformRng rng;
 
public:
    explicit ImpactCone(const AffineMatrix& frame, const Float cutoffVelocity)
        : frame(frame)
        , v_c(cutoffVelocity) {}
 
    void getFragments(const Float m_tot, const Size N, Array<Vector>& r, Array<Vector>& v, Array<Float>& m) {
        constexpr Float theta = PI / 4._f;
        const Float m_frag = m_tot / N;
        for (Size i = 0; i < N; ++i) {
            const Float phi = 2._f * PI * rng();
            v.push(frame * sphericalToCartesian(v_c, theta, phi));
            r.push(frame.translation());
            m.push(m_frag);
        }
    }
};
 
class CollisionMC {
private:
    UniformRng rng;
 
public:
    Array<Float> operator()(const Float QoverQ_D, const Float M_tot, const Float m_min) {
        const Float largest = max(M_LR(QoverQ_D, M_tot), M_LF(QoverQ_D, M_tot));
        const Float exponent = q(QoverQ_D) + 1._f;
 
        Array<Float> fragments;
        fragments.push(largest);
        Float m_partial = largest;
        while (M_tot - m_partial > m_min) {
            const Float m = pow(rng(), 1._f / exponent) - m_min; 
            if (m + m_partial > M_tot) {
                continue;
            }
            fragments.push(m);
            m_partial += m;
        }
        return fragments;
    }
 
private:
    Float M_LR(const Float QoverQ_D, const Float M_tot) {
        if (QoverQ_D < 1._f) {
            return (-0.5_f * (QoverQ_D - 1._f) + 0.5_f) * M_tot;
        } else {
            return (-0.35_f * (QoverQ_D - 1._f) + 0.5_f) * M_tot;
        }
    }
 
    Float M_LF(const Float QoverQ_D, const Float M_tot) {
        return 8.e-3_f * (QoverQ_D * exp(-sqr(0.25_f * QoverQ_D))) * M_tot;
    }
 
    Float q(const Float QoverQ_D) {
        return -10._f + 7._f * pow(QoverQ_D, 0.4_f) * exp(-QoverQ_D / 7._f);
    }
};
 
NAMESPACE_SPH_END