Q is now again a private variable of the nodes

lib/node.cpp

@@ -55,6 +55,7 @@ void Node::update(List localParticles, unsigned nLocalParticles) {
     }
 
     updateMassAndCOM();
+    calculateQ();
 
     isLeaf_ = true;
     for (unsigned i=0; i<8; ++i) {
@@ -145,14 +146,8 @@ Vector Node::calculateForce(unsigned particle) {
     return force;
 }
 
-
-Vector Node::multipole(Vector y, double ynorm) const {
-    Vector result = {0, 0, 0};
-    for (unsigned i=0; i<3; ++i) {
-        result[i] += totalMass_*mass_/(ynorm*ynorm) * y[i]/ynorm;
-    }
-    // calculate matrix Q
-    Matrix Q = {{0,0,0},{0,0,0},{0,0,0}};
+void Node::calculateQ() {
+    Q_ = {{0,0,0},{0,0,0},{0,0,0}};
     for (unsigned i=0; i<3; ++i) {
         for (unsigned j=0; j<3; ++j) {
             // calculate the entry in Qij
@@ -171,7 +166,14 @@ Vector Node::multipole(Vector y, double ynorm) const {
             }
         }
     }
+}
 
+Vector Node::multipole(Vector y, double ynorm) const {
+    Vector result = {0, 0, 0};
+    for (unsigned i=0; i<3; ++i) {
+        result[i] += totalMass_*mass_/(ynorm*ynorm) * y[i]/ynorm;
+    }
+    /*
     // y^T * Q * y
     double yTQy = 0;
     for (unsigned i=0; i<3; ++i) {
@@ -191,7 +193,7 @@ Vector Node::multipole(Vector y, double ynorm) const {
     for (unsigned i=0; i<3; ++i) {
         result[i] += 0.5*(-5 * yTQy / ynorm + QQTy[i] * y[i])/(ynorm * ynorm*ynorm*ynorm*ynorm);
     }
-
+    */
     return result;
 }
 
@@ -211,6 +213,10 @@ void Node::updateMassAndCOM() {
     }
 }
 
+Vector Node::getCOM() const {
+    return centerOfMass_;
+}
+
 unsigned Node::getOctant(unsigned particle) {
     double px = positions_[3*particle + 0];
     double py = positions_[3*particle + 1];
@@ -246,10 +252,6 @@ 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) {

lib/node.hpp

@@ -27,15 +27,17 @@ public:
     void update(List particles, unsigned nParticles);
     // Calculate force on a particle recursively
     Vector calculateForce(unsigned particle);
-    // Calculate the matrix Q
+    // Calculate matrix Q for multipole expansion
+    void calculateQ();
+    // Calculate multipole expansion of a node
     Vector multipole(Vector y, double ynorm) const;
     // Update center of mass and total mass in the node
     void updateMassAndCOM();
+    // get center of mass
+    Vector getCOM() const;
     // Function that takes a particle index in the global particles_ array and returns the octant it is in relative to the center fo the current node 
     // this function doesn't check if the particle is actually in the node.
     unsigned getOctant(unsigned particle);
-    // get center of mass
-    Vector getCOM() const;
     // save node to file
     void save2file(std::ofstream& file) const;
 
@@ -58,8 +60,10 @@ private:
     List childParticles_[8];
     unsigned nLocalParticles_;
 
+    Matrix Q_;
+
     double totalMass_;
-    const double theta0_ = 0.0;
+    const double theta0_ = 0.4;
     const double minNParticles_ = 8;
     bool isLeaf_;
 };