Octree and Barnes Hut Implementation

.github/workflows/coverage.yml

@@ -16,9 +16,13 @@ jobs:
       run: |
         sudo apt-get update -q
         sudo apt-get install libyaml-cpp-dev -y
+        sudo apt-get install bc -y
 
     - name: Evaluate coverage
       run: |
+        make coverage > coverage_output.txt
+        cat coverage_output.txt
+
         COVERAGE=$(make coverage | grep 'TOTAL COVERAGE' | grep -Eo '[0-9]+')
 
         if [ "$COVERAGE" -gt 89 ]; then

Makefile

@@ -1,5 +1,5 @@
 CXX = g++
-CXXFLAGS = -std=c++17 -Wall --coverage
+CXXFLAGS = -g -std=c++17 -Wall --coverage
 LDFLAGS = -lyaml-cpp
 SRC_DIR = src
 INC_DIR = include
@@ -19,7 +19,7 @@ TEST_SRCS = $(wildcard $(TEST_DIR)/*.cpp)
 TEST_OBJS = $(patsubst $(TEST_DIR)/%.cpp,$(TEST_OBJ_DIR)/%.o,$(TEST_SRCS))
 
 # Include directories
-INC_DIRS = -I $(INC_DIR)
+INC_DIRS = -I $(INC_DIR) #-I /opt/homebrew/Cellar/yaml-cpp/0.8.0/include
 
 # Linking step for src files
 $(BIN_DIR)/$(TARGET): $(OBJS)
@@ -62,7 +62,7 @@ coverage:
 	for filename in $(filter-out $(SRC_DIR)/simulation.cpp, $(SRCS)); do \
 		obj_file=$(OBJ_DIR)/$${filename%.cpp}.o; \
 		result=$$(gcov --object-directory=$(OBJ_DIR) $${filename} -n | grep -v ".*simulation.*" | grep -v ".*\.h" | grep -A 1 "src"); \
-		results+="\n\n$$result"; \
+		results="$$results\n\n$$result"; \
 		lines=$$(echo $$result | grep -Eo "[0-9]+$$"); \
 		line_pct=$$(echo $$result | grep -Eo '([0-9]+\.[0-9]+)\%' | sed 's/\%//'); \
 		total_lines=$$((total_lines + lines)); \

include/body.h

@@ -11,7 +11,7 @@ class Body : public Particle {
     public:
 
         // Child class constructor
-        explicit Body(const std::array<double, 3>* initial_position, const std::array<double, 3>* initial_velocity, double mass, float radius=0);
+        explicit Body(const std::array<float, 3>* initial_position, const std::array<float, 3>* initial_velocity, float mass, float radius=0);
 
         // New physical member
         float radius;

include/environment.h

@@ -2,39 +2,50 @@
 #include <vector>
 #include <array>
 #include <map>
-#include <memory> // Include for smart pointers
-#include "particle.h"
+#include <memory>
+#include <functional>
 
+#include "./particle.h"
+#include "./octree.h"
 
 template <typename T>
-class GravitationalEnvironment {
-public:
-    // Constructors
-    GravitationalEnvironment(const std::vector<std::shared_ptr<T>>& particlePtrs, const bool log, std::string logFilePrefix = "run");
-    GravitationalEnvironment(const std::string configFileName, const bool log, std::string logFilePrefix = "run");
+class GravitationalEnvironment{
+
+    public:
+        // Constructors
+        GravitationalEnvironment(const std::vector<std::shared_ptr<T>>& particlePtrs, const bool log, std::string logFilePrefix="run", std::string forceAlgorithm="pair-wise");
+        GravitationalEnvironment(const std::string configFileName, const bool log, std::string logFilePrefix = "run", std::string forceAlgorithm="pair-wise");
 
-    // Define member functions
-    void loadParticlesFromConfig(std::string configFileName);
-    std::vector<std::array<double, 3>> getForces(const double timestep);
-    void updateAll(const std::vector<std::array<double, 3>>& forces, const double timestep);
-    void step(const double timestep);
-    void simulate(const double duration, const double timestep);
-    std::string getStepLog() const;
-    std::string getLogHeader() const;
-    void reset();
+        // Callable member that we will set to pair-wise or Barnes-Hut force algorithm
+        std::function<std::vector<std::array<float, 3>>(float)> getForces;
 
-    // Instantiation of the physical members
-    std::vector<std::shared_ptr<T>> particlePtrs; // Changed to std::shared_ptr
-    bool log;
-    double time;
-    int nParticles; // This should be calculated in the constructor.
-    std::string logFileName;
+        // Define member functions for force algorithms
+        std::vector<std::array<float, 3>> getForcesPairWise(const float timestep);
+        std::vector<std::array<float, 3>> getForcesBarnesHut(const float timestep);
+
+        void loadParticlesFromConfig(std::string configFileName);
+        std::array<float, 3> calculateForceBarnesHut(std::shared_ptr<T> objPtr, std::shared_ptr<Octree<T>> currPtr, std::array<float, 3> netForce, float theta);
+        void updateAll(const std::vector<std::array<float, 3>>& forces, const float timestep);
+        void step(const float timestep);
+        void simulate(const float duration, const float timestep);
+        std::string getStepLog() const;
+        std::string getLogHeader() const;
+        void reset();
+
+        // Instantiation of the physical members
+        std::vector<std::shared_ptr<T>> particlePtrs;
+        bool log;
+        float time;
+        int nParticles;
+        std::string logFileName;
+        Octree<T> envOctree;
 
-private:
-    // Instantiation of the physical members
-    std::string logFilePrefix;
+    private:
+        // Instantiation of the physical members
+        std::string logFilePrefix;
 };
 
 // Helper functions
 int getLargestLabelNumber(const std::vector<std::string>& filenames, const std::string logFilePrefix);
+float getEuclidianDistance(std::array<float, 3> coords1, std::array<float, 3> coords2);
 std::map<std::string, std::map<std::string, std::string>> loadConfig(const std::string& fileName);

include/octree.h

@@ -10,31 +10,34 @@ class Octree {
     public:
 
         // Octree constructor
-        Octree(std::array<float, 2>& xCoords, std::array<float, 2>& yCoords, std::array<float, 2>& zCoords);
+        Octree(std::array<float, 2>& xCoords, std::array<float, 2>& yCoords, std::array<float, 2>& zCoords, bool internal);
 
         // Member functions
-        void clear();
+        void clearOctree();
+        void updateCoords(std::array<float, 2>& newXCoords, std::array<float, 2>& newYCoords, std::array<float, 2>& newZCoords);
         void insert(std::shared_ptr<T> objPtr);
         void build(std::vector<std::shared_ptr<T>>& objPtrs);
 
         // Members
         std::vector<std::shared_ptr<T>> objPtrs;
         std::array<float, 3> centerOfMass;
-        float* totalMass;
+        float totalMass;
+        bool internal;
 
         // Dimensions of the current octant
         std::array<float, 2> xCoords;
         std::array<float, 2> yCoords;
         std::array<float, 2> zCoords;
 
         // Octree children --> 0-7 based on 2D convention in postive z, and then 2D convention in negative z, observing from above        
-        Octree<T>* child0;
-        Octree<T>* child1;
-        Octree<T>* child2;
-        Octree<T>* child3;
-        Octree<T>* child4;
-        Octree<T>* child5;
-        Octree<T>* child6;
-        Octree<T>* child7;
+        std::shared_ptr<Octree<T>> child0;
+        std::shared_ptr<Octree<T>> child1;
+        std::shared_ptr<Octree<T>> child2;
+        std::shared_ptr<Octree<T>> child3;
+        std::shared_ptr<Octree<T>> child4;
+        std::shared_ptr<Octree<T>> child5;
+        std::shared_ptr<Octree<T>> child6;
+        std::shared_ptr<Octree<T>> child7;
+
 
 };

include/particle.h

@@ -7,13 +7,13 @@ class Particle {
 
     public:
         // Constructors
-        Particle(const std::array<double, 3>* initial_position, const std::array<double, 3>* initial_velocity, double mass);
+        Particle(const std::array<float, 3>* initial_position, const std::array<float, 3>* initial_velocity, float mass);
 
         // Define member functions
-        void update (const std::array<double, 3>* force, const double timestep);
+        void update (const std::array<float, 3>* force, const float timestep);
 
         // Instantiation of the physical members
-        std::array<double, 3> position;
-        std::array<double, 3> velocity;
-        double mass;
+        std::array<float, 3> position;
+        std::array<float, 3> velocity;
+        float mass;
 };

include/statistics.h

@@ -8,11 +8,14 @@ extern std::mt19937 GENERATOR;
 
 // Distribution template used to sample from random distributions
 template <typename Distribution>
-std::vector<double> sampleFromDistribution(size_t n, Distribution& distribution) {
+std::vector<float> sampleFromDistribution(size_t n, Distribution& distribution) {
     // Generate samples
-    std::vector<double> samples(n);
+    std::vector<float> samples(n);
     for (size_t i = 0; i < n; ++i) {
         samples[i] = distribution(GENERATOR);
     }
     return samples;
 }
+
+// Dummy executable to work with testing framework
+int dummyExecutable();

src/body.cpp

@@ -11,5 +11,5 @@
 float DENSITY = 4E3;
 
 // Constructor definition
-Body::Body(const std::array<double, 3>* initial_position, const std::array<double, 3>* initial_velocity, double mass, float radius)
+Body::Body(const std::array<float, 3>* initial_position, const std::array<float, 3>* initial_velocity, float mass, float radius)
     : Particle(initial_position, initial_velocity, mass), radius((radius == 0) ? pow(3 * mass / (4 * M_PI * DENSITY), 1./3.) : radius) {};