2d working without MBD

src/Calculation.cpp

@@ -9,7 +9,7 @@
 void Calculation::run(size_t steps) {
     for (size_t step = 0; step < steps; ++step) {
         m_integrator(rod, sim);
-        for (auto &comp : computes) {
+        for (auto &comp: computes) {
             comp.eval(rod);
         }
     }
@@ -20,28 +20,21 @@ const Rod2d &Calculation::getRod() const { return rod; }
 Calculation::Calculation(
         std::function<void(Rod2d &, Simulation &)> t_integrator,
         std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
-    double deltaT, size_t seed,double length)
+        double deltaT, size_t seed, double length)
-    : rod(length), sim(deltaT, seed),m_integrator(std::move(t_integrator)) {
+        : rod(length), sim(deltaT, seed), m_integrator(std::move(t_integrator)) {
-    for (const auto &pair : t_computes) {
+    for (const auto &pair: t_computes) {
         computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
     }
 }
 
-Calculation::Calculation(
+Calculation::Calculation(inte_force_type t_integrator,
-    std::function<
-        void(Rod2d &, Simulation &,
-             std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>,
-             std::function<double(Eigen::Vector2d, Eigen::Vector2d)>)>
-        t_integrator,
                          std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
-    double deltaT, size_t seed,
+                         double deltaT, size_t seed, force_type force, torque_type torque, double length)
-    std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
+        : rod(length), sim(deltaT, seed) {
-    std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque,  double length)
-    :rod(length), sim(deltaT, seed) {
     m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) {
         t_integrator(t_rod, t_sim, force, torque);
     };
-    for (const auto &pair : t_computes) {
+    for (const auto &pair: t_computes) {
         computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
     }
 }
@@ -51,3 +44,24 @@ const std::vector<Compute> &Calculation::getComputes() const {
 }
 
 const Simulation &Calculation::getSim() const { return sim; }
+
+Calculation::Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator, const std::vector<std::pair<Compute::Type,size_t>>&  t_computes,
+                         double deltaT, size_t seed, double length) : rod(length), sim(deltaT, seed),
+                                                                      m_integrator(std::move(t_integrator)){
+    for (const auto &pair: t_computes) {
+        computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
+    }
+}
+
+Calculation::Calculation(Calculation::inte_force_type t_integrator, std::vector<std::pair<Compute::Type,size_t>> t_computes, double deltaT,
+                         size_t seed, Calculation::force_type force, Calculation::torque_type torque, double length)
+        : rod(length), sim(deltaT, seed) {
+    m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) {
+        t_integrator(t_rod, t_sim, force, torque);
+    };
+    for (const auto &pair: t_computes) {
+        computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
+    }
+    rod.setPos({1000,1000});
+
+}

src/Calculation.h

@@ -10,30 +10,32 @@
 #include "Simulation.h"
 
 class Calculation {
-   private:
+private:
     Rod2d rod;
     Simulation sim;
     std::function<void(Rod2d &, Simulation &)> m_integrator;
     std::vector<Compute> computes;
 
-   public:
+public:
     [[nodiscard]] const Simulation &getSim() const;
+
     using force_type = std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>;
     using torque_type = std::function<double(Eigen::Vector2d, Eigen::Vector2d)>;
     using inte_force_type = std::function<void(Rod2d &, Simulation &, force_type, torque_type)>;
-    Calculation(
-      inte_force_type
-            t_integrator,
-        std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
-        double deltaT, size_t seed,
-        force_type force,
-        torque_type torque, double length = 1.0);
 
-    explicit Calculation(
+    Calculation(inte_force_type t_integrator, std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
-        std::function<void(Rod2d &, Simulation &)> t_integrator,
+                double deltaT, size_t seed, force_type force, torque_type torque, double length = 1.0);
+
+    Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator,
                 std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
                 double deltaT, size_t seed, double length = 1.0);
 
+    Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator,
+                const std::vector<std::pair<Compute::Type,size_t>>&  t_computes, double deltaT, size_t seed, double length = 1.0);
+
+    Calculation(inte_force_type t_integrator, std::vector<std::pair<Compute::Type,size_t>>  t_computes,
+                double deltaT, size_t seed, force_type force, torque_type torque, double length = 1.0);
+
     [[nodiscard]] const std::vector<Compute> &getComputes() const;
 
     void run(size_t steps);

src/Compute.cpp

@@ -28,7 +28,7 @@ void Compute::eval_empXX(const Rod2d &rod2D) {
     const auto &new_pos = rod2D.getPos();
     auto old_pos = start_rod.getPos();
     auto old_e = start_rod.getE();
-    double empxx = pow((new_pos - old_pos).dot(old_e), 2.0)/2.0 / time;
+    double empxx = pow((new_pos - old_pos).dot(old_e), 2.0) / 2.0 / time;
     agg.feed(empxx);
 }
 
@@ -37,7 +37,7 @@ void Compute::eval_empYY(const Rod2d &rod2D) {
     auto old_pos = start_rod.getPos();
     auto old_e = start_rod.getE();
     Eigen::Vector2d old_e_ortho = {-old_e[1], old_e[0]};
-    double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0)/2.0 / time;
+    double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0) / 2.0 / time;
     agg.feed(empxx);
 }
 
@@ -57,6 +57,9 @@ void Compute::eval(const Rod2d &rod2D) {
             case Type::empyy:
                 eval_empYY(rod2D);
                 break;
+            case Type::msd_harmonic_k_10:
+                evalMSD(rod2D);
+                break;
         }
         start_rod = rod2D;
     }
@@ -66,10 +69,11 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
         : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type) {
     time = sim.getMDeltaT() * static_cast<double>(every);
     switch (type) {
-        case Type::msd:{
+        case Type::msd: {
             type_str = "msd";
             target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
-            break;}
+            break;
+        }
         case Type::oaf: {
             type_str = "oaf";
             target = std::exp(-rod.getDRot() * time);
@@ -83,7 +87,8 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
             target = Dmean - deltaD / 2.0 *
                              (1 - exp(-u * rod.getDRot() * time)) / u /
                              rod.getDRot() / time;
-        } break;
+        }
+            break;
         case Type::empyy: {
             type_str = "empyy";
             const double Dmean = 0.5 * (rod.getDiff().trace());
@@ -91,8 +96,14 @@ Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
             const double deltaD = rod.getDiff()(1, 1) - rod.getDiff()(0, 0);
             target = Dmean + deltaD / 2.0 *
                              (1 - exp(-u * rod.getDRot() * time)) / u /
-                                 rod.getDRot() / time ;
+                             rod.getDRot() / time;
-        } break;
+        }
+            break;
+        case Type::msd_harmonic_k_10: {
+            type_str = "msd_harmonic";
+            target = 2.0/10.0*(1.0 - std::exp(-20 * time));
+        }
+            break;
     }
 }
 

src/Compute.h

@@ -12,7 +12,7 @@
 
 class Compute {
    public:
-    enum class Type { msd, oaf, empxx, empyy };
+    enum class Type { msd, oaf, empxx, empyy, msd_harmonic_k_10 };
 
     explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim);
 

src/Integratoren2d_force.cpp

@@ -32,7 +32,7 @@ void Integratoren2d_force::Set1_Euler(
     // Force
     Vector F_lab = force(rod2D.getPos(), e);
     Vector F_part = rotation_Matrix(e).inverse() * F_lab;
-    Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
+    Vector F_trans = rotation_Matrix(e)*rod2D.getDiff_Sqrt() * F_part;
     F_trans *= sim.getMDeltaT();
 
     // rotation
@@ -59,7 +59,7 @@ Vector Integratoren2d_force::Heun_predictor_pos(
 
     Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE());
     Vector F_pred_part = rotation_Matrix(rod2D.getE()).inverse() * F_pred_lab;
-    Vector F_pred_trans = rod2D.getDiff_Sqrt() * F_pred_part;
+    Vector F_pred_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_pred_part;
     F_pred_trans *= sim.getMDeltaT() ;
 
     return F_pred_trans + trans_pred_lab;
@@ -96,7 +96,6 @@ void Integratoren2d_force::Set2_Heun(
 
     Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque);
 
-
     Vector e_integrated =
         rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future);
     // apply
@@ -117,7 +116,7 @@ void Integratoren2d_force::Set3_Exact(
     Vector F_lab = force(rod2D.getPos(), rod2D.getE());
     Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
 
-    Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
+    Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part;
     F_trans *= sim.getMDeltaT();
     Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
 
@@ -149,7 +148,7 @@ void Integratoren2d_force::Set4_BDAS(
     // Force
     Vector F_lab = force(rod2D.getPos(), rod2D.getE());
     Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
-    Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
+    Vector F_trans = rotation_Matrix(rod2D.getE())*rod2D.getDiff_Sqrt() * F_part;
     F_trans *= sim.getMDeltaT();
 
     // rotation

src/main.cpp

@@ -3,6 +3,7 @@
 //
 
 #include <iostream>
+#include <fstream>
 
 #include "Calculation.h"
 #include "Integratoren2d_force.h"
@@ -10,7 +11,7 @@
 
 int main() {
     constexpr int numStep = 10000000;
-    constexpr double k = 0.01;
+    constexpr double k = 10;
     [[maybe_unused]] auto harmonic_Force =
             [](const Eigen::Vector2d &pos,
                const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
@@ -25,27 +26,40 @@ int main() {
                                            const Eigen::Vector2d & /*torque*/) -> double {
         return 0.0;
     };
+    std::vector<std::pair<Compute::Type, size_t>> computes;
+    for (int i = 1; i < 100; ++i) {
+        computes.emplace_back(Compute::Type::msd_harmonic_k_10, i * 1);
+//computes.emplace_back(Compute::Type::oaf, i * 1);
+        //computes.emplace_back(Compute::Type::empxx, i * 1);
+        //computes.emplace_back(Compute::Type::empyy, i * 1);
+    }
+
+    Calculation euler(Integratoren2d_force::Set4_BDAS,
+                      computes, 0.01, 12345, harmonic_Force, zero_Torque, 1.0);
 
-    constexpr Compute::Type compute = Compute::Type::oaf;
-    Calculation euler(Integratoren2d_force::Set2_Heun,
-                      {{compute, 1},
-                       {compute, 2},
-                       {compute, 4},
-                       {compute, 8},
-                       {compute, 16},
-                       {compute, 32},
-                       {compute, 64},
-                       {compute, 128},
-                       {compute, 256},
-                       {compute, 512},
-                       {compute, 1024}},
-                      0.01, 12345, zero_Force,zero_Torque,1.0);
     euler.run(numStep);
-    for (const auto &com : euler.getComputes()) {
+
-        if (com.getType() != compute) { continue;
+    std::ofstream msdFile("msd.dat");
-}
+    std::ofstream msd_harmFile("msd_harm.dat");
-        std::cout
+    std::ofstream oafFile("oaf.dat");
-                  << com.getAgg().getMean() << " , " << com.getTarget()
+    std::ofstream empxxFile("empxx.dat");
-                  << std::endl;
+    std::ofstream empyyFile("empyy.dat");
+
+    for (const auto &com: euler.getComputes()) {
+        if (com.getType() == Compute::Type::msd) {
+            msdFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
+        }
+        if (com.getType() == Compute::Type::oaf) {
+            oafFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
+        }
+        if (com.getType() == Compute::Type::empxx) {
+            empxxFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
+        }
+        if (com.getType() == Compute::Type::empyy) {
+            empyyFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
+        }
+        if (com.getType() == Compute::Type::msd_harmonic_k_10) {
+            msd_harmFile << com.getTime() << ", " << com.getAgg().getMean() << ", " << com.getTarget() << std::endl;
+        }
     }
 }