Octree and Barnes Hut Implementation (#20)

* made progress on octree class; working on building the octree class

* finished implementing octree class and members/functions; need to fix segmentation fault in testing

* update testing and previous implementation; debug octree implementation and makefile configuration

* add inclusion of <memory> to octree files

* add tests for body implementation; slightly edit body.cpp implementation

* fixed implementation of octree

* fixed octree insertion implementation; added testing for octree insert, but need more coverage

* fixed basic pointers vs shared pointers in octree implementation; added an envOctree member to the Grav Environment class; added some tests

* fixed issue with types; need to add more coverage and implement barnes-hut algo walking down the tree

* change every double to a float; implemented barnes-hut force algo; need to test

* implemented barnes-hut force calculation; need to implement testing

* Add functional header

* added testing for barnes-hut force calculation algorithm

* add an include header of math.h to test_environment

* fix type issue in test of barnes-hut algo

* added more testing for octree class; changed totalMass to a float from a pointer

* added line to workaround makefile issue on John's laptop; changed to floats from doubles; implemented work with Adam's config stuff

* change include to INC_DIR variable; fix warning in test_statistics.cpp

* small change to push and test workflow

* commented out test case that detects for errors in loading file

* add line to coverage to debug

* added debugging print to coverage workflow

* added dummy executable for the statistics file; changed += from Makefile; added download of bc in coverage workflow

* adjut Makefile for John's local machine; testing automation on github

* Fix small typo in Makefile to fix automation (LDLFLAGS --> LDFLAGS)

* change initilization of coordinate looping in environment to at first, position of first pointer

---------

Co-authored-by: Adam Boesky <[email protected]>

.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,7 @@
 CXX = g++
-CXXFLAGS = -std=c++17 -Wall --coverage
+CXXFLAGS = -g -std=c++17 -Wall --coverage
+# LINE BELOW REQUIRED FOR JOHN'S LOCAL CONFIGURATIONS #
+# LDFLAGS = -L/opt/homebrew/Cellar/yaml-cpp/0.8.0/lib -lyaml-cpp
 LDFLAGS = -lyaml-cpp
 SRC_DIR = src
 INC_DIR = include
@@ -19,6 +21,8 @@ TEST_SRCS = $(wildcard $(TEST_DIR)/*.cpp)
 TEST_OBJS = $(patsubst $(TEST_DIR)/%.cpp,$(TEST_OBJ_DIR)/%.o,$(TEST_SRCS))
 
 # Include directories
+# LINE BELOW REQUIRED FOR JOHN'S LOCAL CONFIGURATIONS #
+# INC_DIRS = -I $(INC_DIR) -I /opt/homebrew/Cellar/yaml-cpp/0.8.0/include
 INC_DIRS = -I $(INC_DIR)
 
 # Linking step for src files
@@ -62,7 +66,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) {};