diff --git a/src/Compute.cpp b/src/Compute.cpp index 85b609c..ce4e90e 100644 --- a/src/Compute.cpp +++ b/src/Compute.cpp @@ -9,85 +9,109 @@ #include "Rod2d.hpp" #include "Simulation.h" -void Compute::evalMSD(const Rod2d &rod2D) { - const auto &new_pos = rod2D.getPos(); - auto old_pos = start_rod.getPos(); - auto msd = (new_pos - old_pos).squaredNorm(); - agg.feed(msd); +void Compute::evalMSD(const Rod2d &rod2D) +{ + const auto &new_pos = rod2D.getPos(); + auto old_pos = start_rod.getPos(); + auto msd = (new_pos - old_pos).squaredNorm(); + agg.feed(msd); } -void Compute::evalOAF(const Rod2d &rod2D) { - const auto &new_e = rod2D.getE(); - auto old_e = start_rod.getE(); - auto oaf = old_e.dot(new_e); - agg.feed(oaf); +void Compute::evalOAF(const Rod2d &rod2D) +{ + const auto &new_e = rod2D.getE(); + auto old_e = start_rod.getE(); + auto oaf = old_e.dot(new_e); + agg.feed(oaf); } -void Compute::eval_empXX(const Rod2d & /*rod2D*/) { - //TODO +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); + agg.feed(empxx); } -void Compute::eval_empYY(const Rod2d & /*rod2D*/) { - //TODO +void Compute::eval_empYY(const Rod2d &rod2D) +{ + const auto &new_pos = rod2D.getPos(); + auto old_pos = start_rod.getPos(); + auto old_e_ortho = start_rod.getE_ortho(); + double empxx = pow((new_pos - old_pos).dot(old_e_ortho), 2.0); + agg.feed(empxx); } -void Compute::eval(const Rod2d &rod2D) { - time_step++; - if (time_step % every == 0) { - switch (type) { - case Type::msd: - evalMSD(rod2D); - break; - case Type::oaf: - evalOAF(rod2D); - break; - case Type::empxx: - eval_empXX(rod2D); - break; - case Type::empyy: - eval_empYY(rod2D); - break; - } - start_rod = rod2D; +void Compute::eval(const Rod2d &rod2D) +{ + time_step++; + if (time_step % every == 0) { + switch (type) { + case Type::msd: + evalMSD(rod2D); + break; + case Type::oaf: + evalOAF(rod2D); + break; + case Type::empxx: + eval_empXX(rod2D); + break; + case Type::empyy: + eval_empYY(rod2D); + break; } - + start_rod = rod2D; + } } 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) { - auto time = sim.getMDeltaT() * static_cast(every); - switch (type) { - case Type::msd: - target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; - break; - case Type::oaf: - target = std::exp(-rod.getDRot() * time); - break; - case Type::empxx: - target = 0; - break; - case Type::empyy: - target = 0; - break; - } + : start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type) +{ + auto time = sim.getMDeltaT() * static_cast(every); + switch (type) { + case Type::msd: + target = 4.0 * 0.5 * (rod.getDiff().trace()) * time; + break; + case Type::oaf: + target = std::exp(-rod.getDRot() * time); + break; + case Type::empxx: { + const double Dmean = 0.5 * (rod.getDiff().trace()); + const double u = 4.0; + 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; + } break; + case Type::empyy: { + const double Dmean = 0.5 * (rod.getDiff().trace()); + const double u = 4.0; + 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; + } break; + } } -const LiveAgg &Compute::getAgg() const { - return agg; +const LiveAgg &Compute::getAgg() const +{ + return agg; } -Compute::Type Compute::getType() const { - return type; +Compute::Type Compute::getType() const +{ + return type; } -double Compute::getTime() const { - return static_cast(every); +double Compute::getTime() const +{ + return static_cast(every); } -double Compute::getTarget() const { - return target; +double Compute::getTarget() const +{ + return target; } -double Compute::getDifference() const { - return abs(agg.getMean() -target); +double Compute::getDifference() const +{ + return abs(agg.getMean() - target); } \ No newline at end of file diff --git a/src/Compute.h b/src/Compute.h index ee1d1ac..95ffa2a 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -10,45 +10,48 @@ #include "Rod2d.hpp" #include "Simulation.h" -class Compute { +class Compute +{ public: - enum class Type { - msd, oaf, empxx, empyy - }; + enum class Type { + msd, + oaf, + empxx, + empyy + }; - explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); + explicit Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim); - void eval(const Rod2d &rod2D); + void eval(const Rod2d &rod2D); - void evalMSD(const Rod2d &rod2D); + void evalMSD(const Rod2d &rod2D); - void evalOAF(const Rod2d &rod2D); + void evalOAF(const Rod2d &rod2D); - void eval_empXX(const Rod2d &rod2D); + void eval_empXX(const Rod2d &rod2D); - double getTarget() const; + double getTarget() const; - void eval_empYY(const Rod2d &rod2D); + void eval_empYY(const Rod2d &rod2D); - [[nodiscard]] const LiveAgg &getAgg() const; + [[nodiscard]] const LiveAgg &getAgg() const; - [[nodiscard]] Type getType() const; + [[nodiscard]] Type getType() const; - double getTime() const; + double getTime() const; - double getDifference() const; + double getDifference() const; private: - Rod2d start_rod; - LiveAgg agg; - size_t every; - size_t time_step; - Type type; - double target; - + Rod2d start_rod; + LiveAgg agg; + size_t every; + size_t time_step; + Type type; + double target; }; -#endif //MYPROJECT_COMPUTE_H +#endif// MYPROJECT_COMPUTE_H diff --git a/test/test_Integratoren.cpp b/test/test_Integratoren.cpp index 0276bfd..62bebd6 100644 --- a/test/test_Integratoren.cpp +++ b/test/test_Integratoren.cpp @@ -16,7 +16,9 @@ TEST_CASE("Forceless Integratoren") { SECTION("Euler") { Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED); + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED); { euler.run(10000); for (auto &c: euler.getComputes()) { CHECK(c.getDifference() <= 0.005); } @@ -34,14 +36,18 @@ TEST_CASE("Forceless Integratoren") { }SECTION("Heun") { Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED); + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED); { heun.run(10000); for (auto &c: heun.getComputes()) { CHECK(c.getDifference() <= 0.005); } } Calculation heun_force(Integratoren2d_force::Set2_Heun, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED, force, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED, force, torque); { heun_force.run(10000); @@ -52,14 +58,18 @@ TEST_CASE("Forceless Integratoren") { }SECTION("Exact") { Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED); + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED); { exact.run(10000); for (auto &c: exact.getComputes()) { CHECK(c.getDifference() <= 0.005); } } Calculation exact_force(Integratoren2d_force::Set3_Exact, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED, force, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED, force, torque); { exact_force.run(10000); @@ -70,20 +80,24 @@ TEST_CASE("Forceless Integratoren") { }SECTION("BDAS") { Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED); + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED); { bdas.run(10000); for (auto &c: bdas.getComputes()) { CHECK(c.getDifference() <= 0.005); } } Calculation bdas_force(Integratoren2d_force::Set4_BDAS, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED, force, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, 0.01, SEED, force, torque); { bdas_force.run(10000); for (auto &c: bdas_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } } - CHECK(bdas.getComputes()[0].getAgg().getMean() == bdas_force.getComputes()[0].getAgg().getMean()); - CHECK(bdas.getComputes()[1].getAgg().getMean() == bdas_force.getComputes()[1].getAgg().getMean()); + CHECK(bdas.getComputes()[0].getAgg().getMean() == Approx(bdas_force.getComputes()[0].getAgg().getMean()).epsilon(10e-10)); + CHECK(bdas.getComputes()[1].getAgg().getMean() == Approx(bdas_force.getComputes()[1].getAgg().getMean()).epsilon(10e-10)); } } diff --git a/test/test_Rod2d.cpp b/test/test_Rod2d.cpp index f3939d6..d233c41 100644 --- a/test/test_Rod2d.cpp +++ b/test/test_Rod2d.cpp @@ -43,5 +43,16 @@ TEST_CASE("Rods") { } + SECTION("Angle Stuff"){ + const Eigen::Vector2d unitVec = {1, 0}; + double f1 = GENERATE(take(10,random(-100,100))); + double f2 = GENERATE(take(10,random(0,100))); + sphere.setE(Eigen::Vector2d({f1,f2}).normalized()); + Eigen::Rotation2D rotMat(acos(unitVec.dot(sphere.getE()))); + CHECK(rotMat.toRotationMatrix()(0,0) == Approx(sphere.getE_Base_matrix()(0,0)).epsilon(10e-10)); + CHECK(rotMat.toRotationMatrix()(1,0) == Approx(sphere.getE_Base_matrix()(1,0)).epsilon(10e-10)); + CHECK(rotMat.toRotationMatrix()(0,1) == Approx(sphere.getE_Base_matrix()(0,1)).epsilon(10e-10)); + CHECK(rotMat.toRotationMatrix()(1,1) == Approx(sphere.getE_Base_matrix()(1,1)).epsilon(10e-10)); + } } \ No newline at end of file