stuff done

C++/CMakeLists.txt

@@ -16,6 +16,8 @@ if(CMAKE_CXX_COMPILER_ID MATCHES ".*Clang")
   endif()
 endif()
 
+
+
 # Link this 'library' to use the warnings specified in CompilerWarnings.cmake
 add_library(project_warnings INTERFACE)
 

C++/conanfile.py

@@ -6,7 +6,6 @@ class CppStarterProject(ConanFile):
     version = "0.1"
     requires = (
         "catch2/2.13.7",
-        "docopt.cpp/0.6.2",
         "fmt/8.0.1",
         "spdlog/1.9.2",
         "eigen/3.4.0",

C++/src/CMakeLists.txt

@@ -5,38 +5,44 @@ find_package(Eigen3)
 
 add_library(logging tools/logging.cpp)
 
-add_library(integratoren SHARED tools/LiveAgg.cpp Data/Rod2d.cpp Data/Simulation.cpp
-        Data/Simulation.h Integratoren/Integratoren2d.cpp
-        tools/Compute.cpp runner/konvergenz_runner.cpp runner/konvergenz_runner.h)
+add_library(integratoren SHARED 
+	tools/LiveAgg.cpp 
+	Data/Rod2d.cpp 
+	Data/Simulation.cpp
+	Data/Rod2dMulti.cpp        
+	Force/ForceMulti.cpp
+	Integratoren/Integratoren2d.cpp 
+	tools/Compute.cpp 
+	runner/konvergenz_runner.cpp)
 target_include_directories(integratoren PUBLIC .)
 target_link_libraries(
-        integratoren
-        PRIVATE project_options
-        project_warnings
-        fmt::fmt
-        spdlog::spdlog
-        Eigen3::Eigen3
-        logging
-    )
+	integratoren
+	PRIVATE project_options
+	project_warnings
+	fmt::fmt
+	spdlog::spdlog
+	Eigen3::Eigen3
+	logging
+	)
 
 target_include_directories(logging PUBLIC .)
 target_link_libraries(
-        logging
-        PRIVATE project_options
-        project_warnings
-        fmt::fmt
-        spdlog::spdlog)
+	logging
+	PRIVATE project_options
+	project_warnings
+	fmt::fmt
+	spdlog::spdlog)
 
 add_executable(main main.cpp)
 find_package(OpenMP)
 if(OpenMP_CXX_FOUND)
-    target_link_libraries(main PUBLIC OpenMP::OpenMP_CXX)
+	target_link_libraries(main PUBLIC OpenMP::OpenMP_CXX)
 endif()
 target_link_libraries(
-        main
-        PRIVATE project_options
-        project_warnings
-        integratoren
-        Eigen3::Eigen3
-        fmt::fmt
-)
+	main
+	PRIVATE project_options
+	project_warnings
+	integratoren
+	Eigen3::Eigen3
+	fmt::fmt
+	)

C++/src/Data/Rod2d.cpp

@@ -1,34 +1,23 @@
 #include "Rod2d.hpp"
 
 #include "Eigen/Dense"
-
-constexpr double M_Pl = 3.141592653589793238462643383279502884; /* pi */
-
-void sqrt(Eigen::Matrix2d &mat) {
-    const auto size = static_cast<size_t>(mat.size());
-    for (size_t i = 0; i < size; i++) {
-        mat.data()[i] = sqrt(mat.data()[i]);
-    }
-}
+#include "tools.hpp"
 Rod2d::Rod2d(double L) : m_pos({1, 1}), m_e({1, 0}) {
     assert(L >= 1.0);
     if (L == 1.0) {
         m_D_rot = 3.0;
         m_Diff << 1, 0, 0, 1;
     } else {
-        const double D0 = 3.0 * M_Pl / L;
+        const double D0 = 3.0 * Tools::M_Pl / L;
         const double p = L;
-        auto D_para =
-            D0 / (M_Pl * 2.0) * (log(p) - 0.1401 + 1.034 / p - 0.228 / (p * p));
-        auto D_ortho = D0 / (M_Pl * 4.0) *
-                       (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p));
-        m_D_rot = 3 * D0 / (M_Pl * L * L) *
-                  (log(p) - 0.3512 + 0.7804 / p - 0.09801 / (p * p));
+        auto D_para = Tools::D_para(D0, p);
+        auto D_ortho = Tools::D_ortho(D0, p);
+        m_D_rot = Tools::D_rot(D0, p);
         m_Diff << D_para, 0, 0, D_ortho;
     }
     m_Diff_sqrt = m_Diff;
     m_D_rot_sqrt = sqrt(m_D_rot);
-    sqrt(m_Diff_sqrt);
+    Tools::mat_sqrt(m_Diff_sqrt);
 }
 
 void Rod2d::setPos(const Eigen::Vector2d &Pos) { m_pos = Pos; }

C++/src/Data/Rod2d.hpp

@@ -1,6 +1,7 @@
 #pragma once
 #include <Eigen/Dense>
 #include <Eigen/Geometry>
+using pos_type = Eigen::Vector2d;
 class Rod2d {
    public:
     explicit Rod2d(double L);

C++/src/Data/Rod2dMulti.cpp

@@ -0,0 +1,46 @@
+#include "Rod2dMulti.hpp"
+#include "Eigen/Dense"
+#include "tools.hpp"
+
+Rod2dMulti::Rod2dMulti(double length, size_t num_particles)
+    : m_pos(num_particles, Eigen::Vector2d({1, 1})),
+      m_e(num_particles, {1, 0}) {
+    assert(L >= 1.0);
+    if (length == 1.0) {
+        m_D_rot = 3;
+        m_Diff << 1, 0, 0, 1;
+    } else {
+        const double D0 = 3.0 * Tools::M_Pl / length;
+        const double p = length;
+        auto D_para = Tools::D_para(D0, p);
+        auto D_ortho = Tools::D_ortho(D0, p);
+        m_D_rot = Tools::D_rot(D0, p);
+        m_Diff << D_para, 0, 0, D_ortho;
+    }
+    m_Diff_sqrt = m_Diff;
+    m_D_rot_sqrt = sqrt(m_D_rot);
+    Tools::mat_sqrt(m_Diff_sqrt);
+}
+
+void Rod2dMulti::setPos(const Eigen::Vector2d &Pos, size_t pos) {
+    assert(pos < m_pos.size());
+    m_pos[pos] = Pos;
+}
+
+double Rod2dMulti::getDRot() const { return m_D_rot; }
+
+const Eigen::Vector2d &Rod2dMulti::getPos(size_t pos) const {
+    return m_pos[pos];
+}
+
+const Eigen::Matrix2d &Rod2dMulti::getDiff() const { return m_Diff; }
+
+const Eigen::Matrix2d &Rod2dMulti::getDiff_Sqrt() const { return m_Diff_sqrt; }
+
+double Rod2dMulti::getDRot_Sqrt() const { return m_D_rot_sqrt; }
+
+const Eigen::Vector2d &Rod2dMulti::getE(size_t pos) const { return m_e[pos]; }
+size_t Rod2dMulti::num_parts() const{
+return m_e.size();
+}
+void Rod2dMulti::setE(const Eigen::Vector2d &mE, size_t pos) { m_e[pos] = mE; }

C++/src/Data/Rod2dMulti.hpp

@@ -0,0 +1,31 @@
+#pragma once
+#include <Eigen/Dense>
+#include <Eigen/Geometry>
+class Rod2dMulti {
+   public:
+    explicit Rod2dMulti(double length, size_t num_particles);
+
+    [[nodiscard]] const Eigen::Matrix2d &getDiff() const;
+
+    [[nodiscard]] double getDRot() const;
+
+    [[nodiscard]] const Eigen::Vector2d &getPos(size_t pos) const;
+
+    void setPos(const Eigen::Vector2d &Pos, size_t pos);
+
+    [[nodiscard]] const Eigen::Matrix2d &getDiff_Sqrt() const;
+
+    [[nodiscard]] const Eigen::Vector2d &getE(size_t pos) const;
+
+    void setE(const Eigen::Vector2d &mE, size_t pos);
+
+    [[nodiscard]] double getDRot_Sqrt() const;
+    [[nodiscard]] size_t num_parts()const;
+   private:
+    Eigen::Matrix2d m_Diff;
+    Eigen::Matrix2d m_Diff_sqrt;
+    double m_D_rot;
+    double m_D_rot_sqrt;
+    std::vector<Eigen::Vector2d> m_pos;  // position
+    std::vector<Eigen::Vector2d> m_e;
+};

C++/src/Data/Simulation.cpp

@@ -4,17 +4,53 @@
 
 #include "Simulation.h"
 
-double Simulation::getNorm(double t_norm) {
-    //return 0.0;
-    return t_norm * m_norm(m_generator);
-}
-
 Simulation::Simulation(double t_delta_T, size_t seed)
-    : m_delta_T(t_delta_T),
+    : m_xBorder(100),
+      m_yBorder(100),
+      m_delta_T(t_delta_T),
       m_std(std::sqrt(t_delta_T * 2.0)),
       m_generator(seed),
       m_norm(0, 1.0) {}
 
+double Simulation::getNorm(double t_norm) {
+    return t_norm * m_norm(m_generator);
+}
+
 double Simulation::getMDeltaT() const { return m_delta_T; }
 
 double Simulation::getSTD() const { return m_std; }
+
+void Simulation::set_xBorder(double border) {
+    m_enable_border = true;
+    m_xBorder = border;
+};
+void Simulation::set_yBorder(double border) {
+    m_enable_border = true;
+    m_yBorder = border;
+};
+pos_type Simulation::boundaryPos(pos_type pos) const {
+    if (not m_enable_border) {
+        return pos;
+    }
+    if (pos[0] > m_xBorder) {
+        pos[0] -= m_xBorder;
+    }
+    if (pos[1] > m_yBorder) {
+        pos[1] -= m_yBorder;
+    }
+    return pos;
+}
+
+pos_type Simulation::boundaryDistance(const pos_type& pos1, const pos_type& pos2) const {
+    pos_type distance = pos2 - pos1;
+    if (not m_enable_border) {
+        return distance;
+    }
+    if (distance[0] > m_xBorder / 2) {
+        distance[0] -= m_xBorder;
+    }
+    if (distance[1] > m_yBorder / 2) {
+        distance[1] -= m_yBorder;
+    }
+    return distance;
+}

C++/src/Data/Simulation.h

@@ -5,19 +5,26 @@
 #pragma once
 
 #include <random>
-
+#include "Rod2d.hpp"
 class Simulation {
    public:
     explicit Simulation(double t_delta_T, size_t seed);
     double getNorm(double t_norm);
+    void set_xBorder(double border);
+    void set_yBorder(double border);
 
    private:
+    bool m_enable_border = false;
+    double m_xBorder;
+    double m_yBorder;
     double m_delta_T;
     double m_std;
     std::mt19937_64 m_generator;
     std::normal_distribution<double> m_norm;
 
-public:
+   public:
     [[nodiscard]] double getMDeltaT() const;
     [[nodiscard]] double getSTD() const;
+    [[nodiscard]] pos_type boundaryPos(pos_type pos) const;
+    [[nodiscard]] pos_type boundaryDistance(const pos_type& pos1, const pos_type& pos2) const;
 };

C++/src/Data/tools.hpp

@@ -0,0 +1,27 @@
+#include <cmath>
+#include "Eigen/Dense"
+class Tools {
+   public:
+    constexpr static void mat_sqrt(Eigen::Matrix2d &mat) {
+        const auto size = static_cast<unsigned long>(mat.size());
+        for (unsigned long index = 0; index < size; index++) {
+            mat.data()[index] = sqrt(mat.data()[index]);
+        }
+    }
+
+    constexpr static double M_Pl =
+        3.141592653589793238462643383279502884; /* pi */
+
+    constexpr static double D_para(double D0, double p) {
+        return D0 / (M_Pl * 2.0) *
+               (log(p) - 0.1401 + 1.034 / p - 0.228 / (p * p));
+    }
+    constexpr static double D_ortho(double D0, double p) {
+        return D0 / (M_Pl * 4.0) *
+               (log(p) + 0.8369 + 0.5551 / p - 0.06066 / (p * p));
+    }
+    constexpr static double D_rot(double D0, double p) {
+        return 3 * D0 / (M_Pl * p * p) *
+               (log(p) - 0.3512 + 0.7804 / p - 0.09801 / (p * p));
+    }
+};      // namespace Tools

C++/src/Force/ForceMulti.cpp

@@ -0,0 +1,41 @@
+#include "ForceMulti.h"
+ForceMulti::ForceMulti(size_t num_particles)
+    : m_force(num_particles, {0.0, 0.0}), m_torque(num_particles, 0.0){};
+
+Eigen::Vector2d ForceMulti::force(Eigen::Vector2d pos1, Eigen::Vector2d pos2,
+                                  const Simulation& sim) {
+    auto distance = sim.boundaryDistance(pos1, pos2);
+    auto r = distance.norm();
+    auto e = distance.normalized();
+    return r * e;
+}
+
+double ForceMulti::torque(Eigen::Vector2d pos1, Eigen::Vector2d e1,
+                          Eigen::Vector2d pos2, Eigen::Vector2d e2,
+                          const Simulation& sim) {
+    auto distance = sim.boundaryDistance(pos1, pos2);
+    auto r = distance.norm();
+    auto e = distance.normalized();
+    return r;
+}
+
+void ForceMulti::force_torque(const Rod2dMulti& rods, const Simulation& sim) {
+    for (size_t index_1 = 0; index_1 < rods.num_parts(); ++index_1) {
+        for (size_t index_2 = index_1 + 1; index_2 < rods.num_parts();
+             ++index_2) {
+            auto f = force(rods.getPos(index_1), rods.getPos(index_2), sim);
+            m_force[index_1] += f;
+            m_force[index_2] -= f;
+
+            auto t = torque(rods.getPos(index_1), rods.getE(index_1),
+                            rods.getPos(index_2), rods.getE(index_2), sim);
+            m_torque[index_1] += t;
+            m_torque[index_2] -= t;
+        }
+    }
+};
+
+const std::vector<double>& ForceMulti::get_torque() const { return m_torque; }
+const std::vector<Eigen::Vector2d>& ForceMulti::get_force() const {
+    return m_force;
+}

C++/src/Force/ForceMulti.h

@@ -0,0 +1,20 @@
+#pragma once
+#include <Eigen/Dense>
+#include "Data/Rod2dMulti.hpp"
+#include "Data/Simulation.h"
+class ForceMulti {
+   private:
+    std::vector<Eigen::Vector2d> m_force;
+    std::vector<double> m_torque;
+
+   public:
+    explicit ForceMulti(size_t num_particles);
+    Eigen::Vector2d force(Eigen::Vector2d pos1, Eigen::Vector2d pos2,
+                          const Simulation& sim);
+    double torque(Eigen::Vector2d pos1, Eigen::Vector2d e1,
+                  Eigen::Vector2d pos2, Eigen::Vector2d e2,
+                  const Simulation& sim);
+    void force_torque(const Rod2dMulti& rods, const Simulation &sim);
+    [[nodiscard]] const std::vector<double>& get_torque() const;
+    [[nodiscard]] const std::vector<Eigen::Vector2d>& get_force() const;
+};

C++/src/Force/force_lambdas.h

@@ -5,21 +5,21 @@
 #pragma once
 [[maybe_unused]] const static auto harmonic_Force =
         [](const Eigen::Vector2d &pos,
-           const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
+           const Eigen::Vector2d & /*e_vec*/) -> Eigen::Vector2d {
             return -1.0 * pos;
         };
 [[maybe_unused]] const static auto const_Force =
         [](const Eigen::Vector2d & /*pos*/,
-           const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
+           const Eigen::Vector2d & /*e_vec*/) -> Eigen::Vector2d {
             return {1,1};
         };
 [[maybe_unused]] const static auto zero_Force =
         [](const Eigen::Vector2d & /*pos*/,
-           const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
+           const Eigen::Vector2d & /*e_vec*/) -> Eigen::Vector2d {
             return {0.0, 0.0};
         };
 [[maybe_unused]] const static auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
-                                       const Eigen::Vector2d & /*torque*/) -> double {
+                                       const Eigen::Vector2d & /*e_vec*/) -> double {
     return 0.0;
 };
 

C++/src/main.cpp

@@ -67,17 +67,17 @@ constexpr int num_threads = 5;
 int main() {
     // TIME
     std::cout<<"no_force"<<std::endl;
-    time_no_force(1.0);
-    time_no_force(3.0);
+    //time_no_force(1.0);
+    //time_no_force(3.0);
     std::cout<<"const_force"<<std::endl;
-    time_const(1.0);
-    time_const(3.0);
+    //time_const(1.0);
+    //time_const(3.0);
     std::cout<<"harmonic_force"<<std::endl;
-    time_harmonic(1.0);
-    time_harmonic(3.0);
+    //time_harmonic(1.0);
+    //time_harmonic(3.0);
 
     std::cout<<"convergence"<<std::endl;
-    oaf();
+    //oaf();
     msd();
     return EXIT_SUCCESS;
 }

C++/src/runner/konvergenz_runner.cpp

@@ -5,7 +5,7 @@
 #include <fstream>
 #include <iostream>
 #include <array>
-#include "force_lambdas.h"
+#include "Force/force_lambdas.h"
 #include <indicators/progress_bar.hpp>
 
 

C++/src/tools/Compute.cpp

@@ -111,7 +111,7 @@ Compute::Compute(Rod2d rod, Type t_type, Force force, size_t t_every,
                 case Force::harmonic_F:
                     target =
                         2 * (1 - 2 * std::exp(-time) + std::exp(-2 * time)) +
-                        2 * (1 - std::exp(-2 * time));
+                        2 * (1 - std::exp(-2 * time))   - 2 * (rod.getDiff().trace()) * time;
                     break;
             }
             break;

C++/test/CMakeLists.txt

@@ -9,8 +9,11 @@ target_link_libraries(catch_main PUBLIC Catch2::Catch2 logging)
 target_link_libraries(catch_main PRIVATE project_options)
 
 set(TEST_FILES
-        test_LiveAgg.cpp
+	test_ForceMulti.cpp
+	test_Simulation.cpp
+	test_LiveAgg.cpp
         test_Rod2d.cpp
+        test_Rod2d_multi.cpp
         test_Compute.cpp)
 
 add_executable(tests ${TEST_FILES})

C++/test/test_ForceMulti.cpp

@@ -0,0 +1,38 @@
+#include <catch2/catch.hpp>
+#include "Data/Rod2dMulti.hpp"
+#include "Force/ForceMulti.h"
+
+TEST_CASE("Forces") {
+    ForceMulti force(100);
+    Rod2dMulti rod(1, 100);
+    Simulation sim(1, 123);
+    sim.set_xBorder(100);
+    Eigen::Vector2d zero_pos{0.0, 0.0};
+    Eigen::Vector2d one_pos{1.0, 0.0};
+    Eigen::Vector2d per_zero_pos{100, 100};
+    SECTION("Basic Force") {
+        auto F00 = force.force(zero_pos, zero_pos, sim);
+        auto F11 = force.force(one_pos, one_pos, sim);
+
+        auto F01 = force.force(zero_pos, one_pos, sim);
+        auto F10 = force.force(one_pos, zero_pos, sim);
+
+        CHECK(F11 == F00);
+        CHECK(F01 == -F10);
+    }
+    SECTION("Basic Torque") {}
+    SECTION("Peridoic Torque") {}
+    SECTION("Complete Calculation") {
+        force.force_torque(rod, sim);
+        const auto F = force.get_force();
+        auto F0 = F[0];
+        for (const auto &f : F) {
+            CHECK(F0 == f);
+        }
+        const auto T = force.get_torque();
+        auto T0 = T[0];
+        for (const auto &f : T) {
+            CHECK(T0 == f);
+        }
+    }
+}

C++/test/test_Rod2d_multi.cpp

@@ -0,0 +1,34 @@
+//
+// Created by jholder on 21.10.21.
+//
+
+#include <Eigen/Dense>
+#include <catch2/catch.hpp>
+
+#include "Data/Rod2dMulti.hpp"
+#include "tools/vec_trafos.h"
+
+TEST_CASE("Multi Rods") {
+    Rod2dMulti sphere(1, 100);
+    Rod2dMulti rod(2, 100);
+    SECTION("Checking translation") {
+        rod.setPos({0, 0}, 0);
+        rod.setPos(rod.getPos(0) + Eigen::Vector2d(1, -1), 0);
+        auto newPos = rod.getPos(0);
+        REQUIRE(newPos[0] == 1);
+        REQUIRE(newPos[1] == -1);
+        REQUIRE(sphere.getDiff_Sqrt()(0, 0) == 1.0);
+        REQUIRE(sphere.getDiff_Sqrt()(0, 1) == 0.0);
+        REQUIRE(sphere.getDiff_Sqrt()(1, 0) == 0.0);
+        REQUIRE(sphere.getDiff_Sqrt()(1, 1) == 1.0);
+    }
+    SECTION("Checking rotation") {
+        double f1 = GENERATE(take(10, random(-100, 100)));
+        double f2 = GENERATE(take(10, random(-100, 100)));
+        rod.setE(Eigen::Vector2d({f1, f2}).normalized(), 0);
+        auto test = Eigen::Vector2d({1, 0});
+        auto e = rod.getE(0);
+        auto erg = (rotation_Matrix(e) * test).normalized();
+        REQUIRE(erg.isApprox(e));
+    }
+}

C++/test/test_Simulation.cpp

@@ -0,0 +1,30 @@
+//
+// Created by jholder on 21.10.21.
+//
+
+#include <Eigen/Dense>
+#include <catch2/catch.hpp>
+
+#include "Data/Simulation.h"
+#include "tools/vec_trafos.h"
+
+TEST_CASE("Simulation") {
+    Simulation sim(0.1, 1234);
+    SECTION("Checking Simulation with non periodic settings") {
+    	pos_type pos_Inside {10,10};
+    	pos_type pos_Outside {200,200};
+	CHECK(pos_Inside == sim.boundaryPos(pos_Inside));
+	CHECK(pos_Outside == sim.boundaryPos(pos_Outside));
+	CHECK(pos_Outside-pos_Inside == sim.boundaryDistance(pos_Inside,pos_Outside));
+    }
+    SECTION("Checking Simulation with periodic settings") {
+    	pos_type pos_Inside {10,10};
+    	pos_type pos_Outside {200,200};
+    	pos_type per_pos_Outside {100,100};
+	sim.set_xBorder(100);
+	sim.set_yBorder(100);
+	CHECK(pos_Inside == sim.boundaryPos(pos_Inside));
+	CHECK(per_pos_Outside == sim.boundaryPos(pos_Outside));
+	//CHECK(pos_Outside-pos_Inside == sim.boundaryDistance(pos_Inside,pos_Outside));
+    }
+}