DTreeForce.h

Go to the documentation of this file.

 
#pragma once
 
#include <ogdf/energybased/dtree/DTreeForceTypes.h>
#include <ogdf/energybased/dtree/DTreeWSPD.h>
 
namespace ogdf {
namespace energybased {
namespace dtree {
 
 
template<int Dim>
class DTreeForce {
public:
    template<int _Dim, typename ForceFunc, bool UseForcePrime>
    friend class DTreeWSPDCallback;
 
    using WSPD = DTreeWSPD<Dim>;
    using Tree = typename WSPD::Tree;
 
    explicit DTreeForce(int numPoints);
 
    virtual ~DTreeForce();
 
    double position(int i, int d) const;
 
    void setPosition(int i, int d, double c);
 
    double mass(int i) const;
 
    void setMass(int i, double m);
 
    double force(int i, int d) const;
 
    double force_prime(int i) const;
 
    template<typename ForceFunc, bool UseForcePrime>
    void computeForces(ForceFunc forceFunc);
 
    int numPoints() const { return m_numPoints; };
 
    const WSPD& wspd() const;
 
    const Tree& tree() const;
 
private:
    struct PointData {
        double mass;
 
        double force[Dim];
 
        double force_prime;
    };
 
    struct NodeData {
        double mass;
 
        double centerOfMass[Dim];
 
        double force[Dim];
 
        double force_prime;
    };
 
    WSPD& wspd();
 
    void allocate();
 
    void deallocate();
 
    void bottomUpPhase(int curr);
 
    void topDownPhase(int curr);
 
    void resetPointForces();
 
    void resetTreeNodes();
 
    PointData* m_pointData;
 
    NodeData* m_nodeData;
 
    int m_numPoints;
 
    WSPD* m_pWSPD;
};
 
template<int Dim, typename ForceFunc, bool UseForcePrime>
class DTreeWSPDCallback : public IWSPD {
public:
    DTreeWSPDCallback(DTreeForce<Dim>& treeForce, ForceFunc forceFunc)
        : m_treeForce(treeForce), m_forceFunc(forceFunc) { }
 
    virtual void onWellSeparatedPair(int a, int b) {
        // force vector
        double delta[Dim];
        double force;
        double force_prime;
 
        double dist = computeDeltaAndDistance<Dim>(m_treeForce.m_nodeData[a].centerOfMass,
                m_treeForce.m_nodeData[b].centerOfMass, delta);
        m_forceFunc(dist, force, force_prime);
#if 0
        m_forceFunc(m_treeForce.m_nodeData[a].centerOfMass, m_treeForce.m_nodeData[b].centerOfMass, force, force_prime);
#endif
 
        // compute the force vector for each dim
        for (int d = 0; d < Dim; d++) {
            m_treeForce.m_nodeData[a].force[d] +=
                    force * delta[d] / dist * m_treeForce.m_nodeData[b].mass;
            m_treeForce.m_nodeData[b].force[d] -=
                    force * delta[d] / dist * m_treeForce.m_nodeData[a].mass;
        };
 
        if (UseForcePrime) {
            m_treeForce.m_nodeData[a].force_prime += force_prime * m_treeForce.m_nodeData[b].mass;
            m_treeForce.m_nodeData[b].force_prime += force_prime * m_treeForce.m_nodeData[a].mass;
        }
    }
 
    DTreeForce<Dim>& m_treeForce;
    ForceFunc m_forceFunc;
};
 
template<int Dim>
DTreeForce<Dim>::DTreeForce(int numPoints) : m_numPoints(numPoints) {
    // get the memory
    allocate();
 
    m_pWSPD->setSeparationFactor(1.0);
}
 
template<int Dim>
DTreeForce<Dim>::~DTreeForce() {
    // free the mem
    deallocate();
}
 
template<int Dim>
void DTreeForce<Dim>::allocate() {
    m_pWSPD = new WSPD(m_numPoints);
    m_nodeData = new NodeData[tree().maxNumNodes()];
    m_pointData = new PointData[m_numPoints];
 
    // initial mass setting for a point
    for (int i = 0; i < m_numPoints; i++) {
        m_pointData[i].mass = 1.0;
    }
}
 
template<int Dim>
void DTreeForce<Dim>::deallocate() {
    delete m_pWSPD;
    delete[] m_nodeData;
    delete[] m_pointData;
}
 
template<int Dim>
const typename DTreeForce<Dim>::WSPD& DTreeForce<Dim>::wspd() const {
    return *m_pWSPD;
}
 
template<int Dim>
typename DTreeForce<Dim>::WSPD& DTreeForce<Dim>::wspd() {
    return *m_pWSPD;
}
 
template<int Dim>
const typename DTreeForce<Dim>::Tree& DTreeForce<Dim>::tree() const {
    return wspd().tree();
}
 
template<int Dim>
double DTreeForce<Dim>::position(int i, int d) const {
    return wspd().point(i).x[d];
}
 
template<int Dim>
void DTreeForce<Dim>::setPosition(int i, int d, double c) {
    wspd().setPoint(i, d, c);
}
 
template<int Dim>
double DTreeForce<Dim>::mass(int i) const {
    return m_pointData[i].mass;
}
 
template<int Dim>
void DTreeForce<Dim>::setMass(int i, double m) {
    m_pointData[i].mass = m;
}
 
template<int Dim>
double DTreeForce<Dim>::force(int i, int d) const {
    return m_pointData[i].force[d];
}
 
template<int Dim>
double DTreeForce<Dim>::force_prime(int i) const {
    return m_pointData[i].force_prime;
}
 
template<int Dim>
template<typename ForceFunc, bool UseForcePrime>
void DTreeForce<Dim>::computeForces(ForceFunc forceFunc) {
    // reset all point forces
    resetPointForces();
 
    if (numPoints() <= 1) {
        return;
    }
 
    // first let the WSPD instance update the quadtre
    wspd().update();
 
    // now do the bottom up phase
    bottomUpPhase(tree().rootIndex());
 
    DTreeWSPDCallback<Dim, ForceFunc, UseForcePrime> wspdCallBack(*this, forceFunc);
 
    // run the WSPD cell cell interaction
    wspd().computeWSPD(&wspdCallBack);
 
    // finally, push down the forces
    topDownPhase(tree().rootIndex());
}
 
template<int Dim>
void DTreeForce<Dim>::bottomUpPhase(int curr) {
    // reset the force and the center of mass
    for (int d = 0; d < Dim; d++) {
        m_nodeData[curr].force[d] = 0.0;
        m_nodeData[curr].force_prime = 0.0;
        m_nodeData[curr].centerOfMass[d] = 0.0;
    }
 
    // reset the mass
    m_nodeData[curr].mass = 0.0;
 
    // if this is an inner node
    if (tree().numChilds(curr)) {
        // iterate over the children
        for (int i = 0; i < tree().numChilds(curr); i++) {
            // the child index
            const int child = tree().child(curr, i);
 
            // compute the size of the subtree
            bottomUpPhase(child);
 
            // node curr
            for (int d = 0; d < Dim; d++) {
                // sum up the center of mass coordinates
                m_nodeData[curr].centerOfMass[d] +=
                        m_nodeData[child].centerOfMass[d] * m_nodeData[child].mass;
            }
 
            // add the mass
            m_nodeData[curr].mass += m_nodeData[child].mass;
        };
    } // else it is a leaf
    else {
        // and the remaining points
        for (int i = 0; i < tree().numPoints(curr); ++i) {
            int pointIndex = tree().point(curr, i);
            // node curr
            for (int d = 0; d < Dim; d++) {
                m_nodeData[curr].centerOfMass[d] += position(pointIndex, d) * mass(pointIndex);
            }
 
            m_nodeData[curr].mass += mass(pointIndex);
        }
    }
    // node curr
    for (int d = 0; d < Dim; d++) {
        m_nodeData[curr].centerOfMass[d] /= m_nodeData[curr].mass;
    }
}
 
template<int Dim>
void DTreeForce<Dim>::topDownPhase(int curr) {
    // if this is an inner node
    if (tree().numChilds(curr)) {
        // iterate over the children
        for (int i = 0; i < tree().numChilds(curr); i++) {
            // the child index
            const int child = tree().child(curr, i);
 
            // node curr
            for (int d = 0; d < Dim; d++) {
                // push the force down to the child
                m_nodeData[child].force[d] += m_nodeData[curr].force[d];
            }
 
            m_nodeData[child].force_prime += m_nodeData[curr].force_prime;
 
            // compute the size of the subtree
            topDownPhase(child);
        };
    } // else it is a leaf
    else {
        // and the remaining points
        for (int i = 0; i < tree().numPoints(curr); ++i) {
            // node curr
            int pointIndex = tree().point(curr, i);
 
            // node curr
            for (int d = 0; d < Dim; d++) {
                // push the force down to the child
                m_pointData[pointIndex].force[d] = m_nodeData[curr].force[d] * mass(pointIndex);
            }
 
            m_pointData[pointIndex].force_prime = m_nodeData[curr].force_prime * mass(pointIndex);
        }
    }
}
 
template<int Dim>
void DTreeForce<Dim>::resetPointForces() {
    for (int i = 0; i < m_numPoints; i++) {
        for (int d = 0; d < Dim; d++) {
            m_pointData[i].force[d] = 0.0;
        }
        m_pointData[i].force_prime = 0.0;
    }
}
 
}
}
}