#include <limits> #include <cmath> #include <iostream> #include <algorithm> #include "Node.hpp" namespace nbody { Node::Node(Tree* _tree):tree(_tree) {} Box Node::getBB() const { return bb; } void Node::setBB(const Box& bb_) { bb = bb_; } //check if node needs to be splitted during tree build bool Node::isSplitable() const { bool result = true; if (bodies.size() <= tree->maxLeafBodies) { result = false; } //this is to prevent errors with collocated particles if (bb.volume() <= std::numeric_limits<float>::epsilon()) { result = false; } return result; } void Node::extendBBforBodies() { bb.extendForBodies(bodies); } void Node::extendBBtoCube() { bb.extendToCube(); } std::vector<Body> Node::getBodies() const { return bodies; } void Node::insertBefore(Node* node) { node->next = this; node->prev = this->prev; this->prev->next = node; this->prev = node; } void Node::insertAfter(Node* node) { this->next->insertBefore(node); } void Node::unlink() { this->next->prev = this->prev; this->prev->next = this->next; this->next = this; this->prev = this; } bool Node::isCorrect() const { if (afterSubtree == nullptr) { std::cerr << "after subtree null\n"; return false; } if (!bb.isCorrectBox()) { std::cerr << "bb wrong\n"; return false; } for (std::size_t i = 0; i < 3; i++) { if (bb.min[i] > bb.max[i]) { std::cerr << "bb " << i << " min " << bb.min[i] << " max " << bb.max[i] << '\n'; return false; } } if (std::any_of(std::begin(bodies), std::end(bodies), [&](const Body& b) {return !isContained(b, bb); })) { std::cerr << "bb out of bounds\n"; return false; } if (!leaf) { Node* current = next; std::size_t children = 0; while (current != nullptr && current != afterSubtree) { current = current->afterSubtree; children++; } if (current == nullptr) { std::cerr << "afterSubtree null\n"; return false; } if (children != tree->numberOfChildren()) { std::cerr << "wrong number of children " << children << '\n'; return false; } current = next; for (std::size_t i = 0; i < tree->numberOfChildren(); i++) { current = current->afterSubtree; } if (current != afterSubtree) { std::cerr << "last sibling afterSubtree inconsistent\n"; return false; } } if (!leaf && !bodies.empty()) { std::cerr << "non-empty inner node\n"; return false; } if (leaf && nextSibling != nullptr && next != nextSibling) { std::cerr << "wrong next sibling\n"; return false; } return true; } //update overall node information void Node::update() { double position[3] = {0.0, 0.0, 0.0}; double mass = 0.0; if (leaf) { for (const auto& body : bodies) { mass += body.mass; for (std::size_t i = 0; i < 3; i++) { position[i] += body.position[i] * body.mass; } } } else { for (Node* node = next; node != nullptr; node = node->nextSibling) { mass += node->representative.mass; } } for (std::size_t i = 0; i < 3; i++) { representative.position[i] = position[i] / mass; } representative.mass = mass; } //get criterion to check if node is sufficient for force evaluation double Node::getL() const { return bb.maxSidelength(); } void Node::print(std::size_t parallelId) const { bb.printBB(parallelId); for (const auto& body : bodies) { std::cout << " "; body.print(parallelId); } } //check if node is sufficient for force evaluation bool Node::sufficientForBody(const Body& body) const { double distance = 0.0; for (std::size_t i = 0; i < 3; i++) { distance += (representative.position[i] - body.position[i]) * (representative.position[i] - body.position[i]); } return distance > std::pow(getL(), 2.0); } //check if node is sufficient for force evaluation for all bodies in box bool Node::sufficientForBox(const Box& box) const { return bb.distanceToBox(box) > getL(); } void Node::setBodies(const std::vector<Body>& bodies_) { bodies = bodies_; } void Node::setBodies(std::vector<Body>&& bodies_) { bodies = std::move(bodies_); } //get local bodies void Node::extractLocalBodiesTo(std::vector<Body>& result) { std::copy_if(std::begin(bodies), std::end(bodies), std::back_inserter(result), [](const Body& b) {return !b.refinement; }); bodies.clear(); } } // namespace nbody