implemented constructor and update()

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], List& localParticles, unsigned nLocalParticles) {
     parent_ = parent;
     particles_ = particles;
     nParticles_ = nParticles;
@@ -40,7 +40,36 @@ Node::Node(Node* root, Node* parent, double* masses, double* particles, unsigned
     // TODO
     // allocate children
     // check if child contains no particles, in that case set child to nullptr
-    // check if node is leaf
+    // check if node is leaf (it contains a single particle), in that case all children are nullptr
+    if (nLocalParticles <= 1) {
+        isLeaf_ = true;
+        for (unsigned i=0; i<8; ++i)
+            children_[i] = nullptr;
+    } else {
+        isLeaf_ = false;
+        
+        // for each local particle, determine which octant it is in
+        List* particleDivision = new List[nLocalParticles];
+        for (unsigned i=0; i<nLocalParticles; ++i) {
+            particleDivision[getOctant(localParticles_[i])].push_back(i);
+        }
+
+        // for each octant, create a child node
+        for (unsigned i=0; i<8; ++i) {
+            if (particleDivision[i].size() == 0) {
+                children_[i] = nullptr;
+            } else {
+                double childSize = size_/2;
+                double childCenter[3];
+                for (unsigned j=0; j<3; ++j) {
+                    childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
+                }
+                children_[i] = new Node(root, this, masses, particles, nParticles, depth+1, 
+                        childSize, childCenter, particleDivision[i], particleDivision[i].size());
+            }
+        }
+        delete[] particleDivision;
+    }
 }
 
 Node::~Node() {
@@ -64,29 +93,55 @@ Node& Node::operator=(const Node& other) {
     return *this;
 }
 
-void Node::update(std::vector<unsigned>& localParticles) {
+void Node::update(List& localParticles) {
     // TODO
     // update center of mass
     // maybe subdivide? 
     // update this node
     // update children
+
+    // for each local particle, determine which octant it is in
+    List* particleDivision = new List[nLocalParticles_];
+    for (unsigned i=0; i<nLocalParticles_; ++i) {
+        particleDivision[getOctant(localParticles_[i])].push_back(i);
+    }
+
     isLeaf_ = true;
     for (unsigned i=0; i<8; ++i) {
         if (children_[i] != nullptr) {
-            // children_[i]->update();
-            isLeaf_ = false;
+            // if child is not a nullptr and particleDivision[i] is empty, delete the child
+            // otherwise update the child
+            if (particleDivision[i].size() == 0) {
+                delete children_[i];
+                children_[i] = nullptr;
+            } else {
+                children_[i]->update(particleDivision[i]);
+                isLeaf_ = false;
+            }
+        } else {
+            // if child is a nullptr and particleDivision[i] is not empty, create a new child
+            if (particleDivision[i].size() != 0) {
+                double childSize = size_/2;
+                double childCenter[3];
+                for (unsigned j=0; j<3; ++j) {
+                    childCenter[j] = center_[j] + (i & (1 << j) ? childSize : -childSize);
+                }
+                children_[i] = new Node(root_, this, masses_, particles_, nParticles_, depth_+1, 
+                        childSize, childCenter, particleDivision[i], particleDivision[i].size());
+                isLeaf_ = false;
+            }
         }
     }
+
+    delete[] particleDivision;
 }
 
 double* Node::calculateForce(unsigned particle) {
     // calculate theta
-    double px = particles_[3*particle + 0];
-    double py = particles_[3*particle + 1];
-    double pz = particles_[3*particle + 2];
-    double lambda = std::sqrt((center_[0] - px)*(center_[0] - px)
-                            + (center_[1] - py)*(center_[1] - py)
-                            + (center_[2] - pz)*(center_[2] - pz));
+    double dx = center_[0] - particles_[3*particle + 0];
+    double dy = center_[1] - particles_[3*particle + 1];
+    double dz = center_[2] - particles_[3*particle + 2];
+    double lambda = std::sqrt(dx*dx + dy*dy + dz*dz);
     double theta = size_/lambda;
 
     // TODO
@@ -103,6 +158,7 @@ double* Node::calculateForce(unsigned particle) {
     } else {
         // TODO
     }
+
     return nullptr;
 }
 
@@ -113,4 +169,39 @@ void Node::calculateCenterOfMass() {
     for (unsigned i : localParticles_)
         for (unsigned j=0; j<3; ++j)
             centerOfMass_[j] += particles_[3*i + j] / nLocalParticles_;
+}
+
+unsigned Node::getOctant(unsigned particle) {
+    double px = particles_[localParticles_[0]];
+    double py = particles_[localParticles_[1]];
+    double pz = particles_[localParticles_[2]];
+    if (px < center_[0]) {
+        if (py < center_[1]) {
+            if (pz < center_[2]) {
+                return 0;
+            } else {
+                return 1;
+            }
+        } else {
+            if (pz < center_[2]) {
+                return 2;
+            } else {
+                return 3;
+            }
+        }
+    } else {
+        if (py < center_[1]) {
+            if (pz < center_[2]) {
+                return 4;
+            } else {
+                return 5;
+            }
+        } else {
+            if (pz < center_[2]) {
+                return 6;
+            } else {
+                return 7;
+            }
+        }
+    }
 }
\ No newline at end of file