BarnesHut.h

Go to the documentation of this file.

#pragma once
 
 
#include "common/ForwardDecl.h"
#include "gravity/IGravity.h"
#include "objects/containers/List.h"
#include "objects/finders/KdTree.h"
#include "objects/geometry/Multipole.h"
#include "physics/Constants.h"
#include "sph/kernel/GravityKernel.h"
#include <atomic>
 
NAMESPACE_SPH_BEGIN
 
enum class MultipoleOrder;
 
struct BarnesHutNode : public KdNode {
    Vector com;
 
    MultipoleExpansion<3> moments;
 
    Float r_open;
 
    BarnesHutNode(const Type& type)
        : KdNode(type) {
#ifdef SPH_DEBUG
        com = Vector(NAN);
        r_open = NAN;
#endif
    }
};
 
 
class BarnesHut : public IGravity {
protected:
    ArrayView<const Vector> r;
 
    Array<Float> m;
 
    KdTree<BarnesHutNode> kdTree;
 
    GravityLutKernel kernel;
 
    Float thetaInv;
 
    MultipoleOrder order;
 
    Size maxDepth;
 
    Float gravityConstant = Constants::gravity;
 
public:
    BarnesHut(const Float theta,
        const MultipoleOrder order,
        const Size leafSize = 25,
        const Size maxDepth = 50,
        const Float gravityConstant = Constants::gravity);
 
    BarnesHut(const Float theta,
        const MultipoleOrder order,
        GravityLutKernel&& kernel,
        const Size leafSize = 25,
        const Size maxDepth = 50,
        const Float gravityConstant = Constants::gravity);
 
    virtual void build(IScheduler& pool, const Storage& storage) override;
 
    virtual void evalAll(IScheduler& pool, ArrayView<Vector> dv, Statistics& stats) const override;
 
    virtual Vector eval(const Vector& r0) const override;
 
    virtual Float evalEnergy(IScheduler& scheduler, Statistics& stats) const override;
 
    virtual RawPtr<const IBasicFinder> getFinder() const override;
 
    MultipoleExpansion<3> getMoments() const;
 
protected:
    Vector evalImpl(const Vector& r0, const Size idx) const;
 
    class NodeTask;
 
    struct TreeWalkState {
        List<Size> checkList;
 
        Array<Size> particleList;
 
        Array<Size> nodeList;
 
        Size depth = 0;
 
        TreeWalkState clone() const;
    };
 
    struct TreeWalkResult {
        std::atomic_int approximatedNodes = { 0 };
 
        std::atomic_int exactNodes = { 0 };
    };
 
    void evalNode(IScheduler& scheduler,
        ArrayView<Vector> dv,
        const Size nodeIdx,
        TreeWalkState data,
        TreeWalkResult& result) const;
 
    void evalParticleList(const LeafNode<BarnesHutNode>& leaf,
        ArrayView<Size> particleList,
        ArrayView<Vector> dv) const;
 
    void evalNodeList(const LeafNode<BarnesHutNode>& leaf,
        ArrayView<Size> nodeList,
        ArrayView<Vector> dv) const;
 
    Vector evalExact(const LeafNode<BarnesHutNode>& node, const Vector& r0, const Size idx) const;
 
    void buildLeaf(BarnesHutNode& node);
 
    void buildInner(BarnesHutNode& node, BarnesHutNode& left, BarnesHutNode& right);
};
 
NAMESPACE_SPH_END