diff --git a/lib/node.cpp b/lib/node.cpp index 7fea2fa0fa96869dc482f447b8c30761f4feb55e..8c3969aa23460da90cbf8bdca3838a8718c0616a 100644 --- a/lib/node.cpp +++ b/lib/node.cpp @@ -19,7 +19,7 @@ Node::Node() { */ Node::Node(Node* root, Node* parent, double* masses, double* particles, unsigned nParticles, unsigned depth, - double size, double center[3], std::vector<unsigned> localParticles, unsigned nLocalParticles) { + double size, double center[3], std::vector<unsigned>& localParticles, unsigned nLocalParticles) { parent_ = parent; particles_ = particles; nParticles_ = nParticles; @@ -37,9 +37,10 @@ Node::Node(Node* root, Node* parent, double* masses, double* particles, unsigned mass_ += masses_[localParticles_[i]]; } + // TODO // allocate children // check if child contains no particles, in that case set child to nullptr - // TODO + // check if node is leaf } Node::~Node() { @@ -63,12 +64,18 @@ Node& Node::operator=(const Node& other) { return *this; } -void Node::update() { - // update this node +void Node::update(std::vector<unsigned>& localParticles) { // TODO + // update center of mass + // maybe subdivide? + // update this node // update children + isLeaf_ = true; for (unsigned i=0; i<8; ++i) { - if (children_[i] != nullptr) children_[i]->update(); + if (children_[i] != nullptr) { + // children_[i]->update(); + isLeaf_ = false; + } } } @@ -81,10 +88,18 @@ double* Node::calculateForce(unsigned particle) { + (center_[1] - py)*(center_[1] - py) + (center_[2] - pz)*(center_[2] - pz)); double theta = size_/lambda; - // if angle too large, call this function for every child. If a child is a nullptr, calculate force on particle from particles in this node - // else calculate force on particle from particles in this node + + // TODO + // If angle too large, call this function for every child which is not a nullptr. + // If the node is a leaf, calculate force on particle from particles in this node. + // If the node is not a leaf, call this function for every child which is not a nullptr. + // If angle small enough, calculate force on particle from this node. if (theta > theta0_) { - // TODO + if (isLeaf_) { + // TODO + } else { + // TODO + } } else { // TODO } diff --git a/lib/node.hpp b/lib/node.hpp index ebb859a7e05b804eaf004327cad6f35ea986f86b..cd56829f06e1ced15c138184e5b83790316ea964 100644 --- a/lib/node.hpp +++ b/lib/node.hpp @@ -10,7 +10,7 @@ public: // whoever creates the node is responsible for deciding which particles are in which octant and for allocating the localParticles array // node will write to localParticles array Node(Node* root, Node* parent, double* masses, double* particles, unsigned nParticles, unsigned depth, - double size, double center[3], std::vector<unsigned> localParticles, unsigned nLocalParticles); + double size, double center[3], std::vector<unsigned>& localParticles, unsigned nLocalParticles); // destuctor ~Node(); // copy constructor @@ -18,7 +18,7 @@ public: // copy assignment Node& operator=(const Node&); // update tree: visit each node, check if it needs to be split or merged (for example if a particle has left the region) - void update(); + void update(std::vector<unsigned>& allParticles); // calculate force on a particle recursively double* calculateForce(unsigned particle); private: @@ -35,6 +35,7 @@ private: unsigned nLocalParticles_; double mass_; const double theta0_ = 0.4; + bool isLeaf_; }; #endif /* NODE_HPP */ \ No newline at end of file diff --git a/lib/tree.cpp b/lib/tree.cpp index 08a382405957f107198f5cf6387879c6edb44193..7ff2bce087b6ed3664cfd366aa16d3273013a51e 100644 --- a/lib/tree.cpp +++ b/lib/tree.cpp @@ -17,6 +17,9 @@ Tree::Tree(double* masses, double* particles, double* forces, unsigned nParticle root_ = new Node(); *root_ = Node(root_, root_, masses_, particles_, nParticles_, 0, size_, center_, localParticles_, nParticles_); + + allParticles_ = std::vector<unsigned>(nParticles_); + for (unsigned i=0; i < nParticles_; ++i) allParticles_[i] = i; } Tree::~Tree() { @@ -29,5 +32,5 @@ void Tree::calculateForce(unsigned particle) { } void Tree::update() { - root_->update(); + root_->update(allParticles_); } \ No newline at end of file diff --git a/lib/tree.hpp b/lib/tree.hpp index 5b55b41e8823d7fdbd5b40efd593ebe609471d61..40b3682dad3e5688d9d1348e7641ead7bb4511dd 100644 --- a/lib/tree.hpp +++ b/lib/tree.hpp @@ -3,6 +3,8 @@ #include "node.hpp" +#include <vector> + class Tree { public: // constructor @@ -22,6 +24,7 @@ private: unsigned nParticles_; double size_; double center_[3]; + std::vector<unsigned> allParticles_; }; #endif /* TREE_HPP */ \ No newline at end of file