From ad9bd9c696a01778d95c3945da0c89b101f5dac3 Mon Sep 17 00:00:00 2001
From: "armindamon.riess" <armindamon.riess@uzh.ch>
Date: Wed, 14 Dec 2022 13:46:37 +0100
Subject: [PATCH] function to save tree, fixed 2 errors in node.cpp

---
 lib/nBodySim.cpp |  8 +++++++-
 lib/nBodySim.hpp |  2 ++
 lib/node.cpp     | 19 ++++++++++++++-----
 lib/node.hpp     |  3 +++
 lib/tree.cpp     |  4 ++++
 lib/tree.hpp     |  3 +++
 main.cpp         |  3 +++
 7 files changed, 36 insertions(+), 6 deletions(-)

diff --git a/lib/nBodySim.cpp b/lib/nBodySim.cpp
index 1701e28..179a270 100644
--- a/lib/nBodySim.cpp
+++ b/lib/nBodySim.cpp
@@ -81,7 +81,7 @@ void nBodySim::runSimulation(double dt, unsigned nSteps) {
 
         std::cout << std::right << "Step "    << std::setw(5)  << i 
                                 // << ": MFM = " << std::setw(12) << calculateMeanForceMagnitude() 
-                                // << ", DEV = " << std::setw(12) << verifyForces() 
+                                << ", DEV = " << std::setw(12) << verifyForces() 
                                 // << ", MCD = " << std::setw(8)  << calculateMeanCenterDistance()
                                 << ", MID = " << std::setw(8)  << calculateMeanInterparticleDistance() 
                                 << ", COM = " << "(" << COM[0] << ", " << COM[1] << ", " << COM[2] << ")"
@@ -342,4 +342,10 @@ void nBodySim::saveState2file(unsigned step, std::ofstream& file) const {
     for (unsigned i=0; i < nParticles_; ++i) {
         file << potential_[i] << std::endl;
     }
+}
+
+
+void nBodySim::saveTree2file(std::string filename) const {
+    std::ofstream file(filename);
+    tree_->save2file(file);
 }
\ No newline at end of file
diff --git a/lib/nBodySim.hpp b/lib/nBodySim.hpp
index daef913..4b6be7d 100644
--- a/lib/nBodySim.hpp
+++ b/lib/nBodySim.hpp
@@ -43,6 +43,8 @@ public:
     void write2file(const double* array, std::ofstream& file, unsigned dim, unsigned skip) const;
     // Write current state of simulation to file
     void saveState2file(unsigned step, std::ofstream& file) const;
+    // save location and size of all nodes in tree
+    void saveTree2file(std::string filename) const;
 
     // Functions to get private variables
     double getMass() { return mass_; }
diff --git a/lib/node.cpp b/lib/node.cpp
index 22ebc3c..5f75e6e 100644
--- a/lib/node.cpp
+++ b/lib/node.cpp
@@ -2,7 +2,7 @@
 
 #include <vector>
 #include <cmath>
-#include <iostream>
+#include <cassert>
 
 Node::Node(double mass, double* positions, unsigned nParticles, double softening, unsigned depth, 
         double size, Vector center, List localParticles, unsigned nLocalParticles) {
@@ -50,8 +50,8 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
     for (unsigned i=0; i<8; ++i) {
         childParticles_[i] = List{};
     }
-    for (unsigned i=0; i<nLocalParticles; ++i) {
-        childParticles_[getOctant(localParticles[i])].push_back(i);
+    for (unsigned i : localParticles) {
+        childParticles_[getOctant(i)].push_back(i);
     }
 
     updateMassAndCOM();
@@ -80,13 +80,13 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
             }
         } else {
             // if child wouldn't contain enough particles, don't create a new child
-            if (childParticles_[i].size() < minNParticles_) continue;
+            if (childParticles_[i].size() <= minNParticles_) continue;
             
             // if child is a nullptr and particleDivision[i] contains enough particles, create a new child
             double childSize = size_/2;
             Vector childCenter(3);
             for (unsigned j=0; j<3; ++j) {
-                childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
+                childCenter[j] = center_[j] + (i & (1 << (2-j)) ? childSize : -childSize);
             }
             children_[i] = new Node(mass_, positions_, nParticles_, softening_, depth_+1, 
                     childSize, childCenter, childParticles_[i], childParticles_[i].size());
@@ -258,4 +258,13 @@ unsigned Node::getOctant(unsigned particle) {
 
 Vector Node::getCOM() const {
     return centerOfMass_;
+}
+
+void Node::save2file(std::ofstream& file) const {
+    file << center_[0] << " " << center_[1] << " " << center_[2] << " " << size_ << " " << nLocalParticles_ << std::endl;
+    for (unsigned i=0; i<8; ++i) {
+        if (children_[i] != nullptr) {
+            children_[i]->save2file(file);
+        }
+    }
 }
\ No newline at end of file
diff --git a/lib/node.hpp b/lib/node.hpp
index b1ab8c2..1b4e732 100644
--- a/lib/node.hpp
+++ b/lib/node.hpp
@@ -2,6 +2,7 @@
 #define NODE_HPP
 
 #include <vector>
+#include <fstream>
 
 using Vector = std::vector<double>;
 using Matrix = std::vector<Vector>;
@@ -34,6 +35,8 @@ public:
     unsigned getOctant(unsigned particle);
     // get center of mass
     Vector getCOM() const;
+    // save node to file
+    void save2file(std::ofstream& file) const;
 
 private:
     // Node* root_;
diff --git a/lib/tree.cpp b/lib/tree.cpp
index 04bd078..c9b84e0 100644
--- a/lib/tree.cpp
+++ b/lib/tree.cpp
@@ -35,4 +35,8 @@ void Tree::update() {
 
 Vector Tree::getCOM() const {
     return root_->getCOM();
+}
+
+void Tree::save2file(std::ofstream& file) const {
+    root_->save2file(file);
 }
\ No newline at end of file
diff --git a/lib/tree.hpp b/lib/tree.hpp
index 90cf1a8..dc034c0 100644
--- a/lib/tree.hpp
+++ b/lib/tree.hpp
@@ -4,6 +4,7 @@
 #include "node.hpp"
 
 #include <vector>
+#include <fstream>
 
 class Tree {
 public:
@@ -18,6 +19,8 @@ public:
     void update();
     // get center of mass
     Vector getCOM() const;
+    // save tree to file
+    void save2file(std::ofstream& file) const;
 
 private:
     Node* root_;
diff --git a/main.cpp b/main.cpp
index 0df677a..0ab77d5 100644
--- a/main.cpp
+++ b/main.cpp
@@ -28,8 +28,11 @@ int main(int argc, char** argv) {
     // double* positions = sim.getPositions();
     unsigned N = sim.getNParticles();
     std::cout << "N = " << N << std::endl;
+    
+    sim.saveTree2file("../out/tree.txt");
 
     try {
+        std::cout << "Running simulation..." << std::endl;
         auto start = std::chrono::high_resolution_clock::now();
         sim.runSimulation(dt, nSteps);
         auto stop = std::chrono::high_resolution_clock::now();
-- 
GitLab