Analysis.h

Go to the documentation of this file.

#pragma once
 
 
#include "objects/finders/NeighbourFinder.h"
#include "objects/wrappers/Expected.h"
#include "objects/wrappers/ExtendedEnum.h"
#include "objects/wrappers/Function.h"
#include "system/Settings.h"
 
NAMESPACE_SPH_BEGIN
 
class Storage;
class Path;
class ILogger;
struct PlotPoint;
enum class QuantityId;
namespace Post {
enum class HistogramId;
}
 
SPH_EXTEND_ENUM(QuantityId, Post::HistogramId);
 
namespace Post {
 
Array<Size> findNeighbourCounts(const Storage& storage, const Float particleRadius);
 
 
enum class ComponentFlag {
    OVERLAP = 0,
 
    SEPARATE_BY_FLAG = 1 << 0,
 
    ESCAPE_VELOCITY = 1 << 1,
 
    SORT_BY_MASS = 1 << 2,
};
 
Size findComponents(const Storage& storage,
    const Float particleRadius,
    const Flags<ComponentFlag> flags,
    Array<Size>& indices);
 
Array<Size> findLargestComponent(const Storage& storage,
    const Float particleRadius,
    const Flags<ComponentFlag> flags);
 
 
struct Tumbler {
    Size index;
 
    Float beta;
};
 
Array<Tumbler> findTumblers(const Storage& storage, const Float limit);
 
enum class MoonEnum {
    LARGEST_FRAGMENT, 
    RUNAWAY,          
    MOON,             
    IMPACTOR,         
    UNOBSERVABLE,     
};
 
Array<MoonEnum> findMoons(const Storage& storage, const Float radius = 1._f, const Float limit = 0._f);
 
Size findMoonCount(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    ArrayView<const Vector> v,
    const Size i,
    const Float radius = 1._f,
    const Float limit = 0._f);
 
Vector getCenterOfMass(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    ArrayView<const Size> idxs = nullptr);
 
SymmetricTensor getInertiaTensor(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    const Vector& r0,
    ArrayView<const Size> idxs = nullptr);
 
SymmetricTensor getInertiaTensor(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    ArrayView<const Size> idxs = nullptr);
 
Vector getAngularFrequency(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    ArrayView<const Vector> v,
    const Vector& r0,
    const Vector& v0,
    ArrayView<const Size> idxs = nullptr);
 
Vector getAngularFrequency(ArrayView<const Float> m,
    ArrayView<const Vector> r,
    ArrayView<const Vector> v,
    ArrayView<const Size> idxs = nullptr);
 
Float getSphericity(IScheduler& scheduler, const Storage& strorage, const Float resolution = 0.05_f);
 
enum class HistogramId {
    RADII = -1,
 
    EQUIVALENT_MASS_RADII = -2,
 
    VELOCITIES = -3,
 
    ROTATIONAL_FREQUENCY = -4,
 
    ROTATIONAL_PERIOD = -5,
 
    ROTATIONAL_AXIS = -6,
};
 
using ExtHistogramId = ExtendedEnum<HistogramId>;
 
enum class HistogramSource {
    COMPONENTS,
 
    PARTICLES,
};
 
 
struct HistogramParams {
 
    Interval range;
 
    Size binCnt = 0;
 
    Float referenceDensity = 2700._f;
 
    Float massCutoff = 0._f;
 
    Float velocityCutoff = INFTY;
 
    bool centerBins = true;
 
    struct ComponentParams {
 
        Float radius = 2._f;
 
        Flags<Post::ComponentFlag> flags = Post::ComponentFlag::OVERLAP;
 
    } components;
 
    Function<bool(Size index)> validator = [](Size UNUSED(index)) { return true; };
};
 
struct HistPoint {
    Float value;
 
    Size count;
 
    bool operator==(const HistPoint& other) const;
};
 
Array<HistPoint> getDifferentialHistogram(ArrayView<const Float> values, const HistogramParams& params);
 
Array<HistPoint> getDifferentialHistogram(const Storage& storage,
    const ExtHistogramId id,
    const HistogramSource source,
    const HistogramParams& params);
 
Array<HistPoint> getCumulativeHistogram(const Storage& storage,
    const ExtHistogramId id,
    const HistogramSource source,
    const HistogramParams& params);
 
 
class LinearFunction {
private:
    Float a, b;
 
public:
    LinearFunction(const Float slope, const Float offset)
        : a(slope)
        , b(offset) {}
 
    INLINE Float operator()(const Float x) const {
        return a * x + b;
    }
 
    Float slope() const {
        return a;
    }
 
    Float offset() const {
        return b;
    }
 
    Float solve(const Float y) const {
        SPH_ASSERT(a != 0._f);
        return (y - b) / a;
    }
};
 
LinearFunction getLinearFit(ArrayView<const PlotPoint> points);
 
 
class QuadraticFunction {
private:
    Float a, b, c;
 
public:
    QuadraticFunction(const Float a, const Float b, const Float c)
        : a(a)
        , b(b)
        , c(c) {}
 
    INLINE Float operator()(const Float x) const {
        return (a * x + b) * x + c;
    }
 
    Float quadratic() const {
        return a;
    }
    Float linear() const {
        return b;
    }
    Float constant() const {
        return c;
    }
 
    StaticArray<Float, 2> solve(const Float y) const {
        SPH_ASSERT(a != 0);
        const Float disc = sqr(b) - 4._f * a * (c - y);
        if (disc < 0._f) {
            return EMPTY_ARRAY;
        } else if (disc == 0._f) {
            return { -b / (2._f * a) };
        } else {
            const Float sqrtDisc = sqrt(disc);
            Float x1 = (-b - sqrtDisc) / (2._f * a);
            Float x2 = (-b + sqrtDisc) / (2._f * a);
            if (x1 > x2) {
                std::swap(x1, x2);
            }
            return { x1, x2 };
        }
    }
};
 
QuadraticFunction getQuadraticFit(ArrayView<const PlotPoint> points);
 
} // namespace Post
 
NAMESPACE_SPH_END