diff --git a/CMakeLists.txt b/CMakeLists.txt index 09d24c8..fa64595 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -4,7 +4,7 @@ if(${CMAKE_VERSION} VERSION_LESS 3.12) cmake_policy(VERSION ${CMAKE_MAJOR_VERSION}.${CMAKE_MINOR_VERSION}) endif() -project(gravity VERSION 0.2.0 +project(gravity VERSION 0.4.0 DESCRIPTION "N-Body-Problem project of https://thoughts-on-cpp.com" LANGUAGES CXX) diff --git a/solver/CMakeLists.txt b/solver/CMakeLists.txt index 424950a..78c0069 100644 --- a/solver/CMakeLists.txt +++ b/solver/CMakeLists.txt @@ -1,5 +1,11 @@ add_library(solver - src/solver.cpp include/particle.h include/particleBuilder.h src/particleBuilder.cpp include/types.h) + src/solver.cpp + include/particle.h + include/particleBuilder.h + src/particleBuilder.cpp + include/types.h + src/types.cpp + src/particle.cpp) target_include_directories(solver PUBLIC diff --git a/solver/include/particle.h b/solver/include/particle.h index 5e46e7b..3421986 100644 --- a/solver/include/particle.h +++ b/solver/include/particle.h @@ -2,16 +2,38 @@ #define PARTICLE_H #include "types.h" +#include -struct Particle { +class Particle { +public: + Particle(); + Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position); + + bool operator==(const Particle &rhs) const; + bool operator!=(const Particle &rhs) const; + + double getMass() const; + + const Vector2D &getAcceleration() const; + + void setAcceleration(const Vector2D &acceleration); + + const Vector2D &getVelocity() const; + + void setVelocity(const Vector2D &velocity); + + const Vector2D &getPosition() const; + + void setPosition(const Vector2D &position); + +private: + static size_t IDCounter; + + size_t id; double mass; Vector2D acceleration; Vector2D velocity; Vector2D position; - - Particle() : mass(0) {} - Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position) - : mass(mass), acceleration(acceleration), velocity(velocity), position(position) {} }; #endif diff --git a/solver/include/particleBuilder.h b/solver/include/particleBuilder.h index 996a543..1cd9eb9 100644 --- a/solver/include/particleBuilder.h +++ b/solver/include/particleBuilder.h @@ -2,6 +2,7 @@ #define PARTICLEBUILDER_H #include "types.h" +#include class Particle; @@ -14,6 +15,7 @@ class ParticleBuilder { ParticleBuilder & acceleration(const Vector2D &acceleration); ParticleBuilder & mass(double mass); Particle build() const; + std::vector build(size_t numberOfParticles); private: double mMass; diff --git a/solver/include/solver.h b/solver/include/solver.h index ae04a36..65fa2ef 100644 --- a/solver/include/solver.h +++ b/solver/include/solver.h @@ -2,12 +2,26 @@ #define SOLVER_H #include +#include "types.h" class Particle; class Solver { public: - std::vector solve(const std::vector &particles, double epsilon) const; + explicit Solver(double mEpsilon); + + std::vector solve(const std::vector &particles) const; + +private: + std::vector calculateAcceleration(const std::vector &particles) const; + std::vector calculateVelocity(const std::vector &particles) const; + std::vector calculatePosition(const std::vector &particles) const; + static Particle AccumulateAcceleration(const std::vector &particles, const Particle &particle); + static double CalculateEquivalentMass(const Particle &particleA, const Particle &particleB); + + static const double G; + static const double EPSILON; + double mEpsilon; }; #endif \ No newline at end of file diff --git a/solver/include/types.h b/solver/include/types.h index e7cac2e..31e00f5 100644 --- a/solver/include/types.h +++ b/solver/include/types.h @@ -7,6 +7,30 @@ struct Vector2D { Vector2D() : x(0), y(0) {} Vector2D(double x, double y) : x(x), y(y) {} + + bool operator==(const Vector2D &rhs) const; + + bool operator!=(const Vector2D &rhs) const; + + double length() const; + + Vector2D& operator-=(const Vector2D& rhs); + + Vector2D& operator+=(const Vector2D& rhs); + + Vector2D& operator*=(const double& rhs); + + Vector2D& operator/=(const double& rhs); }; +Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs); + +Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs); + +Vector2D operator*(const Vector2D &lhs, const double &rhs); + +Vector2D operator*(const double &lhs, const Vector2D &rhs); + +Vector2D operator/(const Vector2D &lhs, const double &rhs); + #endif diff --git a/solver/src/particle.cpp b/solver/src/particle.cpp new file mode 100644 index 0000000..655dd20 --- /dev/null +++ b/solver/src/particle.cpp @@ -0,0 +1,54 @@ +#include "particle.h" +#include +#include +#include +#include + +size_t Particle::IDCounter = 0; + +Particle::Particle() + : mass(0) { + id = IDCounter++; +} + +Particle::Particle(double mass, const Vector2D &acceleration, const Vector2D &velocity, const Vector2D &position) + : mass(mass), acceleration(acceleration), velocity(velocity), position(position) { + assert(mass > 0); + id = IDCounter++; +} + +bool Particle::operator==(const Particle &rhs) const { + return id == rhs.id; +} + +bool Particle::operator!=(const Particle &rhs) const { + return !(rhs == *this); +} + +const Vector2D &Particle::getAcceleration() const { + return acceleration; +} + +void Particle::setAcceleration(const Vector2D &acceleration) { + Particle::acceleration = acceleration; +} + +const Vector2D &Particle::getVelocity() const { + return velocity; +} + +void Particle::setVelocity(const Vector2D &velocity) { + Particle::velocity = velocity; +} + +const Vector2D &Particle::getPosition() const { + return position; +} + +void Particle::setPosition(const Vector2D &position) { + Particle::position = position; +} + +double Particle::getMass() const { + return mass; +} diff --git a/solver/src/particleBuilder.cpp b/solver/src/particleBuilder.cpp index 237d7a7..9b995d2 100644 --- a/solver/src/particleBuilder.cpp +++ b/solver/src/particleBuilder.cpp @@ -1,5 +1,6 @@ #include "particleBuilder.h" #include "particle.h" +#include ParticleBuilder & ParticleBuilder::position(const Vector2D &position) { @@ -29,3 +30,21 @@ Particle ParticleBuilder::build() const { return {mMass, mAcceleration, mVelocity, mPosition}; } + +std::vector ParticleBuilder::build(size_t numberOfParticles) +{ + std::vector particle(numberOfParticles); + + std::mt19937 mt(std::random_device{}()); + std::uniform_real_distribution real_dist(1.0, static_cast(numberOfParticles)); + + for(Particle &p : particle) { + p = mass(real_dist(mt)) + .acceleration({ real_dist(mt), real_dist(mt) }) + .velocity({ real_dist(mt), real_dist(mt) }) + .position({ real_dist(mt), real_dist(mt) }) + .build(); + } + + return particle; +} diff --git a/solver/src/solver.cpp b/solver/src/solver.cpp index f5a6135..249c8ab 100644 --- a/solver/src/solver.cpp +++ b/solver/src/solver.cpp @@ -1,7 +1,85 @@ #include "solver.h" #include "particle.h" +#include +#include +#include +#include -std::vector Solver::solve(const std::vector &particles, double epsilon) const -{ - return std::vector(2) ; -} \ No newline at end of file + +const double Solver::G = 6.67408e-11; +const double Solver::EPSILON = 1e-3; + +Solver::Solver(double mEpsilon) : mEpsilon(mEpsilon) {} + +std::vector Solver::solve(const std::vector &particles) const { + std::vector solution = calculateAcceleration(particles); + solution = calculateVelocity(solution); + solution = calculatePosition(solution); + + return solution; +} + +std::vector Solver::calculateAcceleration(const std::vector &particles) const { + std::vector solution(particles.size()); + + std::transform(begin(particles), end(particles), begin(solution), [&particles](const Particle &particle) { + return AccumulateAcceleration(particles, particle); + }); + + return solution; +} + +std::vector Solver::calculateVelocity(const std::vector &particles) const { + std::vector solution(particles.size()); + + std::transform(begin(particles), end(particles), begin(solution), [epsilon = mEpsilon](Particle particle) { + const Vector2D v0 = particle.getVelocity(); + particle.setVelocity(v0 + particle.getAcceleration()*epsilon); + + return particle; + }); + + return solution; +} + +std::vector Solver::calculatePosition(const std::vector &particles) const { + std::vector solution(particles.size()); + + std::transform(begin(particles), end(particles), begin(solution), [epsilon = mEpsilon](Particle particle) { + const Vector2D v = particle.getVelocity(); + const Vector2D r0 = particle.getPosition(); + particle.setPosition(r0 + v*epsilon + particle.getAcceleration()*epsilon*epsilon); + + return particle; + }); + + return solution; +} + +Particle Solver::AccumulateAcceleration(const std::vector &particles, const Particle &particle) { + Particle particleA = particle; + const double e3 = EPSILON*EPSILON*EPSILON; + + std::for_each(begin(particles), end(particles), [&particleA, e3](const Particle &particleB) { + if(particleA != particleB) { + const double M = CalculateEquivalentMass(particleA, particleB); + const Vector2D r = particleB.getPosition() - particleA.getPosition(); + const double rLength = r.length(); + const double r3 = fabs(rLength*rLength*rLength); + + const Vector2D a0 = particleA.getAcceleration(); + particleA.setAcceleration(a0 + G*M*r/(r3 + e3)); + } + }); + + return particleA; +} + +double Solver::CalculateEquivalentMass(const Particle &particleA, const Particle &particleB) { + const double massA = particleA.getMass(); + assert(massA > 0); + + const double massB = particleB.getMass(); + + return massA*massB/massA; +} diff --git a/solver/src/types.cpp b/solver/src/types.cpp new file mode 100644 index 0000000..72b4fd2 --- /dev/null +++ b/solver/src/types.cpp @@ -0,0 +1,67 @@ +#include "types.h" + +#include +#include +#include +#include + +bool Vector2D::operator==(const Vector2D &rhs) const { + auto equalZero = std::numeric_limits::min(); + + return fabs(x - rhs.x) <= equalZero && + fabs(y - rhs.y) <= equalZero; +} + +bool Vector2D::operator!=(const Vector2D &rhs) const { + return !(rhs == *this); +} + +double Vector2D::length() const { + return sqrt(x*x + y*y); +} + +Vector2D& Vector2D::operator-=(const Vector2D& rhs) { + *this = *this - rhs; + return *this; +} + +Vector2D& Vector2D::operator*=(const double& rhs) { + *this = *this * rhs; + return *this; +} + +Vector2D& Vector2D::operator/=(const double& rhs) { + *this = *this / rhs; + return *this; +} + +Vector2D &Vector2D::operator+=(const Vector2D &rhs) { + *this = *this + rhs; + return *this; +} + +Vector2D operator-(const Vector2D &lhs, const Vector2D &rhs) { + return { lhs.x - rhs.x, + lhs.y - rhs.y }; +} + +Vector2D operator+(const Vector2D &lhs, const Vector2D &rhs) { + return { lhs.x + rhs.x, + lhs.y + rhs.y }; +} + +Vector2D operator*(const Vector2D &lhs, const double &rhs) { + return { lhs.x * rhs, + lhs.y * rhs }; +} + +Vector2D operator*(const double &lhs, const Vector2D &rhs) { + return rhs * lhs; +} + +Vector2D operator/(const Vector2D &lhs, const double &rhs) { + assert(fabs(rhs) > 0); + + return { lhs.x / rhs, + lhs.y / rhs }; +} \ No newline at end of file diff --git a/solverTest/CMakeLists.txt b/solverTest/CMakeLists.txt index d1ab1e9..e474563 100644 --- a/solverTest/CMakeLists.txt +++ b/solverTest/CMakeLists.txt @@ -1,9 +1,11 @@ add_executable(solverTest src/test-main.cpp - src/solverTest.cpp) + src/solverTest.cpp src/vectorTest.cpp) + target_link_libraries(solverTest Catch2 solver) enable_testing() -add_test(NAME solverTest COMMAND solverTest) \ No newline at end of file +add_test(NAME solverTest COMMAND solverTest) +add_test(NAME vectorTest COMMAND vectorTest) \ No newline at end of file diff --git a/solverTest/src/solverTest.cpp b/solverTest/src/solverTest.cpp index 56d0fa9..eb8f024 100644 --- a/solverTest/src/solverTest.cpp +++ b/solverTest/src/solverTest.cpp @@ -7,8 +7,7 @@ double Inverse(double value) { return -value; } -Particle GenerateStandardParticle(double xPosition, double yPosition) -{ +Particle GenerateStandardParticle(double xPosition, double yPosition) { ParticleBuilder particleBuilder; return particleBuilder @@ -17,9 +16,9 @@ Particle GenerateStandardParticle(double xPosition, double yPosition) .build(); } -TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) -{ - Solver solver; +TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) { + const double EPSILON = 0.1; + Solver solver(EPSILON); const std::vector particlesX = { GenerateStandardParticle(0.5, 0), GenerateStandardParticle(-0.5, 0)}; @@ -27,38 +26,87 @@ TEST_CASE( "Explicit euler algorithm with two point mass", "[euler]" ) const std::vector particlesY = { GenerateStandardParticle(0, 0.5), GenerateStandardParticle(0, -0.5)}; - const double EPSILON = 1e-3; - //Solution const double acceleration = -0.6674079993; //m/s^2 const double velocity = -0.06674079993; //m/s const double position = 0.48665184; //m SECTION( "Two still standing point mass are attracting each other in x-direction" ) { - std::vector result = solver.solve(particlesX, EPSILON); + std::vector result = solver.solve(particlesX); Particle &particle = result.front(); - REQUIRE(particle.acceleration.x == Approx(acceleration)); - REQUIRE(particle.velocity.x == Approx(velocity)); - REQUIRE(particle.position.x == Approx(position)); + CHECK(particle.getAcceleration().x == Approx(acceleration)); + CHECK(particle.getVelocity().x == Approx(velocity)); + CHECK(particle.getPosition().x == Approx(position)); particle = result.back(); - REQUIRE(particle.acceleration.x == Approx(Inverse(acceleration))); - REQUIRE(particle.velocity.x == Approx(Inverse(velocity))); - REQUIRE(particle.position.x == Approx(Inverse(position))); + CHECK(particle.getAcceleration().x == Approx(Inverse(acceleration))); + CHECK(particle.getVelocity().x == Approx(Inverse(velocity))); + REQUIRE(particle.getPosition().x == Approx(Inverse(position))); } SECTION( "Two still standing point mass are attracting each other in y-direction" ) { - std::vector result = solver.solve(particlesY, EPSILON); + std::vector result = solver.solve(particlesY); Particle &particle = result.front(); - REQUIRE(particle.acceleration.y == Approx(acceleration)); - REQUIRE(particle.velocity.y == Approx(velocity)); - REQUIRE(particle.position.y == Approx(position)); + CHECK(particle.getAcceleration().y == Approx(acceleration)); + CHECK(particle.getVelocity().y == Approx(velocity)); + CHECK(particle.getPosition().y == Approx(position)); particle = result.back(); - REQUIRE(particle.acceleration.y == Approx(Inverse(acceleration))); - REQUIRE(particle.velocity.y == Approx(Inverse(velocity))); - REQUIRE(particle.position.y == Approx(Inverse(position))); + CHECK(particle.getAcceleration().y == Approx(Inverse(acceleration))); + CHECK(particle.getVelocity().y == Approx(Inverse(velocity))); + REQUIRE(particle.getPosition().y == Approx(Inverse(position))); + } +} + +TEST_CASE("Benchmarking euler", "[benchmark]") { + const double EPSILON = 0.1; + Solver solver(EPSILON); + + ParticleBuilder particleBuilder; + std::vector particles = particleBuilder.build(100); + BENCHMARK("Benchmarking with 100 particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(200); + BENCHMARK("Benchmarking with 200 particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(400); + BENCHMARK("Benchmarking with 400 particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(800); + BENCHMARK("Benchmarking with 800 particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(1600); + BENCHMARK("Benchmarking with 1.6K particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(3200); + BENCHMARK("Benchmarking with 3.2K particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(6400); + BENCHMARK("Benchmarking with 6.4K particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(12800); + BENCHMARK("Benchmarking with 12.8K particles") { + particles = solver.solve(particles); + } + + particles = particleBuilder.build(25600); + BENCHMARK("Benchmarking with 25.6K particles") { + particles = solver.solve(particles); } } \ No newline at end of file diff --git a/solverTest/src/vectorTest.cpp b/solverTest/src/vectorTest.cpp new file mode 100644 index 0000000..ae5367c --- /dev/null +++ b/solverTest/src/vectorTest.cpp @@ -0,0 +1,73 @@ +#include +#include +#include + +TEST_CASE( "Test Vector2D operators and functions", "[vector]" ) { + const Vector2D vecA = { 1, 1 }; + const Vector2D vecB = { 3, 3 }; + + SECTION( "Comparing vectors" ) { + const Vector2D expected = vecA; + REQUIRE(vecA == expected); + } + + SECTION( "Length of a vector" ) { + REQUIRE(vecB.length() == Approx(sqrt(3*3 + 3*3))); + } + + SECTION( "Adding two vectors" ) { + const Vector2D expected = { 4, 4 }; + CHECK((vecA + vecB) == expected); + + Vector2D result = vecA; + result += vecB; + REQUIRE(result == expected); + } + + SECTION( "Subtracting two vectors" ) { + const Vector2D expected = { 2, 2 }; + CHECK((vecB - vecA) == expected); + + Vector2D result = vecB; + result -= vecA; + REQUIRE(result == expected); + } + + SECTION( "Multiplying vector with scalar" ) { + const Vector2D expected = { 2, 2 }; + CHECK((vecA * 2) == expected); + + Vector2D result = vecA; + result *= 2; + REQUIRE(result == expected); + } + + SECTION( "Dividing vector with scalar" ) { + const Vector2D expected = { 2, 2 }; + CHECK((vecA * 2) == expected); + + Vector2D result = vecA; + result *= 2; + REQUIRE(result == expected); + } +} + +TEST_CASE( "Negative test Vector2D operators and functions", "[vector]" ) { + const Vector2D vecA = { 1, 1 }; + + SECTION( "Comparing vectors" ) { + Vector2D error = { 0, 1 }; + CHECK_FALSE(vecA == error); + + error = { 1, 0}; + CHECK_FALSE(vecA == error); + + error = { 0.999, 0.999 }; + REQUIRE_FALSE(vecA == error); + } + + SECTION( "Dividing vector with scalar" ) { + //Vector2D failing = vecA / 0; + //Catching asserts not possible with CATCH2 https://github.com/catchorg/Catch2/issues/853 + } +} \ No newline at end of file