From 362ec38dd2e93749d05d51e2ab6b6a23cf239999 Mon Sep 17 00:00:00 2001 From: Jacob Date: Mon, 25 Oct 2021 15:45:29 +0200 Subject: [PATCH] Clang-format --- .clang-format | 146 ++++++++++++++------------ benchmark/bench.cpp | 33 +++--- src/Calculation.cpp | 48 +++++---- src/Calculation.h | 38 +++---- src/Compute.cpp | 163 +++++++++++++---------------- src/Compute.h | 54 ++++------ src/Integratoren2d_force.cpp | 138 ++++++++++++------------ src/Integratoren2d_force.h | 49 +++++---- src/Integratoren2d_forceless.cpp | 113 +++++++++++--------- src/Integratoren2d_forceless.h | 10 +- src/LiveAgg.cpp | 42 +++----- src/LiveAgg.hpp | 31 +++--- src/Rod2d.cpp | 57 ++++------ src/Rod2d.hpp | 13 +-- src/Simulation.cpp | 18 ++-- src/Simulation.h | 15 ++- src/main.cpp | 88 ++++++++++------ test/catch_main.cpp | 4 +- test/test_Calculation.cpp | 48 +++++---- test/test_Compute.cpp | 11 +- test/test_Integratoren.cpp | 173 +++++++++++++++++++------------ test/test_LiveAgg.cpp | 7 +- test/test_Rod2d.cpp | 57 +++++----- test/test_Simulation.cpp | 15 +-- 24 files changed, 708 insertions(+), 663 deletions(-) diff --git a/.clang-format b/.clang-format index 4ede542..807386f 100644 --- a/.clang-format +++ b/.clang-format @@ -1,12 +1,15 @@ -AccessModifierOffset: -2 -AlignAfterOpenBracket: DontAlign +--- +Language: Cpp +# BasedOnStyle: Google +AccessModifierOffset: -1 +AlignAfterOpenBracket: Align AlignConsecutiveAssignments: false AlignConsecutiveDeclarations: false AlignEscapedNewlines: Left -AlignOperands: true -AlignTrailingComments: false -AllowAllParametersOfDeclarationOnNextLine: false -AllowShortBlocksOnASingleLine: true +AlignOperands: true +AlignTrailingComments: true +AllowAllParametersOfDeclarationOnNextLine: true +AllowShortBlocksOnASingleLine: false AllowShortCaseLabelsOnASingleLine: false AllowShortFunctionsOnASingleLine: All AllowShortIfStatementsOnASingleLine: true @@ -14,85 +17,92 @@ AllowShortLoopsOnASingleLine: true AlwaysBreakAfterDefinitionReturnType: None AlwaysBreakAfterReturnType: None AlwaysBreakBeforeMultilineStrings: true -AlwaysBreakTemplateDeclarations: false -BinPackArguments: false -BinPackParameters: false -BraceWrapping: - AfterClass: true +AlwaysBreakTemplateDeclarations: true +BinPackArguments: true +BinPackParameters: true +BraceWrapping: + AfterClass: false AfterControlStatement: false - AfterEnum: false - AfterFunction: true - AfterNamespace: false + AfterEnum: false + AfterFunction: false + AfterNamespace: false AfterObjCDeclaration: false - AfterStruct: true - AfterUnion: false - BeforeCatch: false - BeforeElse: false - IndentBraces: false - SplitEmptyFunction: false - SplitEmptyNamespace: true + AfterStruct: false + AfterUnion: false + BeforeCatch: false + BeforeElse: false + IndentBraces: false + SplitEmptyFunction: true SplitEmptyRecord: true -BreakAfterJavaFieldAnnotations: true -BreakBeforeBinaryOperators: NonAssignment -BreakBeforeBraces: Custom -BreakBeforeInheritanceComma: true + SplitEmptyNamespace: true +BreakBeforeBinaryOperators: None +BreakBeforeBraces: Attach +BreakBeforeInheritanceComma: false BreakBeforeTernaryOperators: true -BreakConstructorInitializers: BeforeColon BreakConstructorInitializersBeforeComma: false +BreakConstructorInitializers: BeforeColon +BreakAfterJavaFieldAnnotations: false BreakStringLiterals: true -ColumnLimit: 0 -CommentPragmas: '^ IWYU pragma:' +ColumnLimit: 80 +CommentPragmas: '^ IWYU pragma:' CompactNamespaces: false -ConstructorInitializerAllOnOneLineOrOnePerLine: false -ConstructorInitializerIndentWidth: 2 -ContinuationIndentWidth: 2 -Cpp11BracedListStyle: false -DerivePointerAlignment: false -DisableFormat: false -ExperimentalAutoDetectBinPacking: true +ConstructorInitializerAllOnOneLineOrOnePerLine: true +ConstructorInitializerIndentWidth: 4 +ContinuationIndentWidth: 4 +Cpp11BracedListStyle: true +DerivePointerAlignment: true +DisableFormat: false +ExperimentalAutoDetectBinPacking: false FixNamespaceComments: true -ForEachMacros: -- foreach -- Q_FOREACH -- BOOST_FOREACH -IncludeCategories: -- Priority: 2 - Regex: ^"(llvm|llvm-c|clang|clang-c)/ -- Priority: 3 - Regex: ^(<|"(gtest|gmock|isl|json)/) -- Priority: 1 - Regex: .* -IncludeIsMainRegex: (Test)?$ -IndentCaseLabels: false -IndentWidth: 2 -IndentWrappedFunctionNames: true +ForEachMacros: + - foreach + - Q_FOREACH + - BOOST_FOREACH +IncludeCategories: + - Regex: '^<.*\.h>' + Priority: 1 + - Regex: '^<.*' + Priority: 2 + - Regex: '.*' + Priority: 3 +IncludeIsMainRegex: '([-_](test|unittest))?$' +IndentCaseLabels: true +IndentWidth: 4 +IndentWrappedFunctionNames: false JavaScriptQuotes: Leave JavaScriptWrapImports: true -KeepEmptyLinesAtTheStartOfBlocks: true -Language: Cpp +KeepEmptyLinesAtTheStartOfBlocks: false MacroBlockBegin: '' -MacroBlockEnd: '' -MaxEmptyLinesToKeep: 2 -NamespaceIndentation: Inner -ObjCBlockIndentWidth: 7 -ObjCSpaceAfterProperty: true +MacroBlockEnd: '' +MaxEmptyLinesToKeep: 1 +NamespaceIndentation: None +ObjCBlockIndentWidth: 4 +ObjCSpaceAfterProperty: false ObjCSpaceBeforeProtocolList: false -PointerAlignment: Right -ReflowComments: true -SortIncludes: false -SortUsingDeclarations: false +PenaltyBreakAssignment: 2 +PenaltyBreakBeforeFirstCallParameter: 1 +PenaltyBreakComment: 300 +PenaltyBreakFirstLessLess: 120 +PenaltyBreakString: 1000 +PenaltyExcessCharacter: 1000000 +PenaltyReturnTypeOnItsOwnLine: 200 +PointerAlignment: Left +ReflowComments: true +SortIncludes: true +SortUsingDeclarations: true SpaceAfterCStyleCast: false -SpaceAfterTemplateKeyword: false +SpaceAfterTemplateKeyword: true SpaceBeforeAssignmentOperators: true SpaceBeforeParens: ControlStatements SpaceInEmptyParentheses: false -SpacesBeforeTrailingComments: 0 -SpacesInAngles: false -SpacesInCStyleCastParentheses: false +SpacesBeforeTrailingComments: 2 +SpacesInAngles: false SpacesInContainerLiterals: true +SpacesInCStyleCastParentheses: false SpacesInParentheses: false SpacesInSquareBrackets: false -Standard: Cpp11 -TabWidth: 8 -UseTab: Never +Standard: Auto +TabWidth: 8 +UseTab: Never +... diff --git a/benchmark/bench.cpp b/benchmark/bench.cpp index cc11b46..528d1ea 100644 --- a/benchmark/bench.cpp +++ b/benchmark/bench.cpp @@ -6,27 +6,22 @@ #define CATCH_CONFIG_ENABLE_BENCHMARKING #include -#include "Rod2d.hpp" #include "Integratoren2d_forceless.h" +#include "Rod2d.hpp" -TEST_CASE("Euler - Baseline", "[benchmark]") -{ +TEST_CASE("Euler - Baseline", "[benchmark]") { Rod2d rod(1.0); Simulation sim(0.01, Catch::rngSeed()); - BENCHMARK("Euler without force") - { - Integratoren2d_forceless::Set1_Euler(rod, sim); - }; - BENCHMARK("Heun without force") - { - Integratoren2d_forceless::Set2_Heun(rod, sim); - }; - BENCHMARK("Exact without force") - { - Integratoren2d_forceless::Set3_Exact(rod, sim); - }; - BENCHMARK("BDAS without force") - { - Integratoren2d_forceless::Set4_BDAS(rod, sim); - }; + BENCHMARK("Euler without force") { + Integratoren2d_forceless::Set1_Euler(rod, sim); + }; + BENCHMARK("Heun without force") { + Integratoren2d_forceless::Set2_Heun(rod, sim); + }; + BENCHMARK("Exact without force") { + Integratoren2d_forceless::Set3_Exact(rod, sim); + }; + BENCHMARK("BDAS without force") { + Integratoren2d_forceless::Set4_BDAS(rod, sim); + }; } \ No newline at end of file diff --git a/src/Calculation.cpp b/src/Calculation.cpp index 8e6fe3f..9eaf12f 100644 --- a/src/Calculation.cpp +++ b/src/Calculation.cpp @@ -6,46 +6,48 @@ #include - void Calculation::run(size_t steps) { for (size_t step = 0; step < steps; ++step) { m_integrator(rod, sim); - for (auto &comp: computes) { comp.eval(rod); } + for (auto &comp : computes) { + comp.eval(rod); + } } } -const Rod2d &Calculation::getRod() const { - return rod; -} +const Rod2d &Calculation::getRod() const { return rod; } -Calculation::Calculation(std::function t_integrator, - std::initializer_list> t_computes, double deltaT, - size_t seed) - : sim(deltaT, seed), m_integrator( - std::move(t_integrator)) { - for (const auto &pair: t_computes) { +Calculation::Calculation( + std::function t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed) + : 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( - std::function, - std::function)> t_integrator, - std::initializer_list> t_computes, double deltaT, - size_t seed, std::function force, - std::function torque) - : 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) { + std::function< + void(Rod2d &, Simulation &, + std::function, + std::function)> + t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, + std::function force, + std::function torque) + : 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)); } } - const std::vector &Calculation::getComputes() const { return computes; } -const Simulation &Calculation::getSim() const { - return sim; -} +const Simulation &Calculation::getSim() const { return sim; } diff --git a/src/Calculation.h b/src/Calculation.h index aff6d62..0078f5b 100644 --- a/src/Calculation.h +++ b/src/Calculation.h @@ -2,35 +2,38 @@ // Created by jholder on 22.10.21. // -#ifndef MYPROJECT_CALCULATION_H -#define MYPROJECT_CALCULATION_H +#pragma once +#include +#include "Compute.h" #include "Rod2d.hpp" #include "Simulation.h" -#include "functional" -#include "Compute.h" - class Calculation { -private: - + private: Rod2d rod = Rod2d(1.0); Simulation sim; std::function m_integrator; std::vector computes; -public: + + public: const Simulation &getSim() const; Calculation( - std::function, - std::function)> t_integrator, - std::initializer_list> t_computes, double deltaT, size_t seed, - std::function force, - std::function torque); + std::function, + std::function)> + t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed, + std::function force, + std::function torque); - explicit Calculation(std::function t_integrator, - std::initializer_list> t_computes, double deltaT, - size_t seed); + explicit Calculation( + std::function t_integrator, + std::initializer_list> t_computes, + double deltaT, size_t seed); [[nodiscard]] const std::vector &getComputes() const; @@ -38,6 +41,3 @@ public: [[nodiscard]] const Rod2d &getRod() const; }; - - -#endif //MYPROJECT_CALCULATION_H diff --git a/src/Compute.cpp b/src/Compute.cpp index ce4e90e..8db6a05 100644 --- a/src/Compute.cpp +++ b/src/Compute.cpp @@ -4,114 +4,99 @@ #include "Compute.h" -#include #include +#include + #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) -{ - 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_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) -{ - 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_empYY(const Rod2d &rod2D) { + const auto &new_pos = rod2D.getPos(); + 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); + 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; +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; } - 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: { - 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; - } + : 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); -} \ No newline at end of file +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 95ffa2a..e121171 100644 --- a/src/Compute.h +++ b/src/Compute.h @@ -5,53 +5,43 @@ #ifndef MYPROJECT_COMPUTE_H #define MYPROJECT_COMPUTE_H - #include "LiveAgg.hpp" #include "Rod2d.hpp" #include "Simulation.h" -class Compute -{ +class Compute { + public: + enum class Type { msd, oaf, empxx, empyy }; -public: - 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; + private: + 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/src/Integratoren2d_force.cpp b/src/Integratoren2d_force.cpp index d92fd23..85e5687 100644 --- a/src/Integratoren2d_force.cpp +++ b/src/Integratoren2d_force.cpp @@ -7,19 +7,16 @@ constexpr double kBT = 1.0; using Vector = Eigen::Vector2d; - -Vector Integratoren2d_force::ortho(Vector e) { - return Vector{-e[1], e[0]}; -} +Vector Integratoren2d_force::ortho(Vector e) { return Vector{-e[1], e[0]}; } Eigen::Matrix2d Integratoren2d_force::MatTraf(Vector e) { Eigen::Matrix2d mat; - mat << e[0], -e[1], - e[1], e[0]; + mat << e[0], -e[1], e[1], e[0]; return mat; } -void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { +void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, + Simulation & /*sim*/) { auto trans_lab = Vector({0, 0.01}); rod2D.setPos(rod2D.getPos() + trans_lab); Eigen::Rotation2D rot(0.1); @@ -28,42 +25,47 @@ void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/ const Vector Integratoren2d_force::unitVec = {1, 0}; +void Integratoren2d_force::Set5_MBD(Rod2d & /*rod2D*/, Simulation & /*sim*/) {} -void Integratoren2d_force::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) { -} - -void Integratoren2d_force::Set1_Euler(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - //translation - Vector rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise - Vector trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System +void Integratoren2d_force::Set1_Euler( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { + auto std = sim.getSTD(); // sqrt(2*delta_T) + // translation + Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise + Vector trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Vector trans_lab = rod2D.getE_Base_matrix() * trans_part; - //Force + // Force Vector F_lab = force(rod2D.getPos(), rod2D.getE()); Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab; Vector F_trans = rod2D.getDiff_Sqrt() * F_part; F_trans *= sim.getMDeltaT() / kBT; - - //rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + // rotation + double rot = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Vector e = rod2D.getE(); - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / kBT * sim.getMDeltaT(); + double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / + kBT * sim.getMDeltaT(); Vector new_e = e + (rot + M_rot) * ortho(e); new_e.normalize(); - //apply + // apply rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); rod2D.setE(new_e); } -Vector Integratoren2d_force::Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim, - const std::function &force) { +Vector Integratoren2d_force::Heun_predictor_pos( + const Rod2d &rod2D, Simulation &sim, + const std::function &force) { auto standard_dev = sim.getSTD(); - Vector rand_pred = {sim.getNorm(standard_dev), sim.getNorm(standard_dev)}; //gaussian noise - Vector trans_pred_part = rod2D.getDiff_Sqrt() * rand_pred; //translation vector in Particle System + Vector rand_pred = {sim.getNorm(standard_dev), + sim.getNorm(standard_dev)}; // gaussian noise + Vector trans_pred_part = + rod2D.getDiff_Sqrt() * + rand_pred; // translation vector in Particle System Vector trans_pred_lab = rod2D.getE_Base_matrix() * trans_pred_part; Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE()); @@ -74,22 +76,24 @@ Vector Integratoren2d_force::Heun_predictor_pos(const Rod2d &rod2D, Simulation & return F_pred_trans + trans_pred_lab; } -Vector Integratoren2d_force::Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim, - const std::function &torque) { +Vector Integratoren2d_force::Heun_predictor_rot( + const Rod2d &rod2D, Simulation &sim, + const std::function &torque) { auto std = sim.getSTD(); - double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + double rot_predict = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Vector e = rod2D.getE(); - double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / kBT * sim.getMDeltaT(); + double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * + rod2D.getDRot() / kBT * sim.getMDeltaT(); Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e); return e_change_predict; } - -void Integratoren2d_force::Set2_Heun(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque) { - +void Integratoren2d_force::Set2_Heun( + Rod2d &rod2D, Simulation &sim, + const std::function &force, + const std::function &torque) { Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force); Vector pos_predictor = rod2D.getPos() + delta_pos_predictor; @@ -97,7 +101,6 @@ void Integratoren2d_force::Set2_Heun(Rod2d &rod2D, Simulation &sim, Vector e_predict = rod2D.getE() + delta_e_predictor; e_predict.normalize(); - Rod2d pred_rod = rod2D; pred_rod.setPos(pos_predictor); pred_rod.setE(e_predict); @@ -105,21 +108,25 @@ void Integratoren2d_force::Set2_Heun(Rod2d &rod2D, Simulation &sim, Vector delta_pos_future = Heun_predictor_pos(pred_rod, sim, force); Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque); - //integration - Vector pos_integrated = 0.5 * rod2D.getPos() + 0.5 * pos_predictor + 0.5 * delta_pos_future; - Vector e_integrated = rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); - //apply + // integration + Vector pos_integrated = + 0.5 * rod2D.getPos() + 0.5 * pos_predictor + 0.5 * delta_pos_future; + Vector e_integrated = + rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future); + // apply rod2D.setPos(pos_integrated); rod2D.setE(e_integrated.normalized()); } -void Integratoren2d_force::Set3_Exact(Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - //translation - Vector rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise - Vector trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System +void Integratoren2d_force::Set3_Exact( + Rod2d &rod2D, Simulation &sim, + std::function force, + std::function torque) { + auto std = sim.getSTD(); // sqrt(2*delta_T) + // translation + Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise + Vector trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Vector F_lab = force(rod2D.getPos(), rod2D.getE()); Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab; @@ -128,25 +135,30 @@ void Integratoren2d_force::Set3_Exact(Rod2d &rod2D, Simulation &sim, F_trans *= sim.getMDeltaT() / kBT; Vector trans_lab = rod2D.getE_Base_matrix() * trans_part; - //rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + // rotation + double rot = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Vector e = rod2D.getE(); - double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / kBT * sim.getMDeltaT(); - auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); + double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / + kBT * sim.getMDeltaT(); + auto correction = + -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e; new_e.normalize(); - //apply + // apply rod2D.setPos(rod2D.getPos() + trans_lab + F_trans); rod2D.setE(new_e); } -void Integratoren2d_force::Set4_BDAS(Rod2d &rod2D, Simulation &sim, - std::function force, - std::function /*torque*/) { +void Integratoren2d_force::Set4_BDAS( + Rod2d &rod2D, Simulation &sim, + std::function force, + std::function /*torque*/) { auto std = sim.getSTD(); - //translation + // translation auto rand = Vector(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + auto trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE()))); auto trans_lab = rotMat * trans_part; @@ -156,16 +168,14 @@ void Integratoren2d_force::Set4_BDAS(Rod2d &rod2D, Simulation &sim, Vector F_trans = rod2D.getDiff_Sqrt() * F_part; F_trans *= sim.getMDeltaT() / kBT; - - auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); + auto rot = + Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); Eigen::Rotation2Dd rotation(rot); auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); - // Normalisation should not be necessary if a proper angular representation is used. - // But with vector e it is done in case of numerical errors + // Normalisation should not be necessary if a proper angular representation + // is used. But with vector e it is done in case of numerical errors rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE(e_new); } - - diff --git a/src/Integratoren2d_force.h b/src/Integratoren2d_force.h index db22fba..d1c56fb 100644 --- a/src/Integratoren2d_force.h +++ b/src/Integratoren2d_force.h @@ -8,23 +8,28 @@ #include "Simulation.h" class Integratoren2d_force { -public: - static void Set1_Euler(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque); + public: + static void Set1_Euler( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque); + static void Set2_Heun( + Rod2d &rod2D, Simulation &sim, + const std::function + &force, + const std::function &torque); - static void Set2_Heun(Rod2d &rod2D, Simulation &sim, - const std::function &force, - const std::function &torque); + static void Set3_Exact( + Rod2d &rod2D, Simulation &sim, + std::function force, + std::function torque); - static void Set3_Exact(Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque); - - static void Set4_BDAS(Rod2d &rod2D, Simulation &sim, - std::function force, - std::function torque); + static void Set4_BDAS( + Rod2d &rod2D, Simulation &sim, + std::function force, + std::function torque); static const Eigen::Vector2d unitVec; @@ -32,17 +37,17 @@ public: static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); -private: + private: static Eigen::Vector2d ortho(Eigen::Vector2d e); static Eigen::Matrix2d MatTraf(Eigen::Matrix e); - static Eigen::Vector2d Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim, - const std::function &force); + static Eigen::Vector2d Heun_predictor_pos( + const Rod2d &rod2D, Simulation &sim, + const std::function + &force); - static Eigen::Vector2d - Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim, - const std::function &torque); + static Eigen::Vector2d Heun_predictor_rot( + const Rod2d &rod2D, Simulation &sim, + const std::function &torque); }; - diff --git a/src/Integratoren2d_forceless.cpp b/src/Integratoren2d_forceless.cpp index cbd4dbc..cbc2c99 100644 --- a/src/Integratoren2d_forceless.cpp +++ b/src/Integratoren2d_forceless.cpp @@ -8,7 +8,8 @@ Eigen::Vector2d Integratoren2d_forceless::ortho(Eigen::Vector2d e) { return Eigen::Vector2d{-e[1], e[0]}; } -void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) { +void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, + Simulation & /*sim*/) { auto trans_lab = Eigen::Vector2d({0, 0.01}); rod2D.setPos(rod2D.getPos() + trans_lab); Eigen::Rotation2D rot(0.1); @@ -16,26 +17,29 @@ void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*s } void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) { - auto std = sim.getSTD(); // sqrt(2*delta_T) - //translation - Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise - Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System - Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; + auto std = sim.getSTD(); // sqrt(2*delta_T) + // translation + Eigen::Vector2d rand = {sim.getNorm(std), + sim.getNorm(std)}; // gaussian noise + Eigen::Vector2d trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System + Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; /* TODO - Statt Zeile 22 und 23 kannst du trans_lab = paralleler_matrix_eintrag * e plus senkrechter_matrix_eintrag * ortho(e); + Statt Zeile 22 und 23 kannst du trans_lab = paralleler_matrix_eintrag * e + plus senkrechter_matrix_eintrag * ortho(e); */ - //rotation - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + // rotation + double rot = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Eigen::Vector2d e = rod2D.getE(); Eigen::Vector2d new_e = e + rot * ortho(e); new_e.normalize(); - //apply + // apply rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE(new_e); - } void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { @@ -43,31 +47,34 @@ void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + Eigen::Vector2d trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab; -/* + /* - TODO + TODO - Siehe Set_1: Hier auch über e und orth(e) + Siehe Set_1: Hier auch über e und orth(e) -*/ + */ - //rotation - double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden. + // rotation + double rot_predict = + sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden. Eigen::Vector2d e = rod2D.getE(); Eigen::Vector2d delta_e_predict = rot_predict * ortho(e); Eigen::Vector2d e_predict = (e + delta_e_predict).normalized(); - double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); - Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); - //TODO pos integration with future system ??? - // Eigen::Vector2d pos_integrated = rod2D.getPos() + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); + Eigen::Vector2d e_integrated = + e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); + // TODO pos integration with future system ??? + // Eigen::Vector2d pos_integrated = rod2D.getPos() + 0.5 * (rot * ortho(e) + // + rot * ortho(e_predict)); - //apply + // apply rod2D.setPos(pos_predictor); rod2D.setE(e_integrated.normalized()); } @@ -76,9 +83,11 @@ constexpr double kBT = 1.0; void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); - //translation - Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise - Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + // translation + Eigen::Vector2d rand = {sim.getNorm(std), + sim.getNorm(std)}; // gaussian noise + Eigen::Vector2d trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part; /* @@ -87,10 +96,11 @@ void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { */ - //rotation + // rotation double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt(); Eigen::Vector2d e = rod2D.getE(); - auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); + auto correction = + -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT(); Eigen::Vector2d delta_e = correction * e + rot * ortho(e); @@ -98,20 +108,21 @@ void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) { TODO Warum hier kBT definiert? */ - //apply + // apply rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE((e + delta_e).normalized()); - } const Eigen::Vector2d Integratoren2d_forceless::unitVec = {1, 0}; -//this is a slow implementation to keep the same data structure for performance analysis this must be altered +// this is a slow implementation to keep the same data structure for performance +// analysis this must be altered void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { auto std = sim.getSTD(); - //translation + // translation auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System + auto trans_part = + rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE()))); auto trans_lab = rotMat * trans_part; @@ -120,17 +131,18 @@ void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { */ - auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); + auto rot = + Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); Eigen::Rotation2Dd rotation(rot); auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized(); /* - TODO Warum Normierung? Ist der Fehler so gros? Wie siehts mit einer Winkelvariable aus, - dann musst du den Winkel nicht jedesmal berechnen? + TODO Warum Normierung? Ist der Fehler so gros? Wie siehts mit einer + Winkelvariable aus, dann musst du den Winkel nicht jedesmal berechnen? */ - // Normalisation should not be necessary if a proper angular representation is used. - // But with vector e it is done in case of numerical errors + // Normalisation should not be necessary if a proper angular representation + // is used. But with vector e it is done in case of numerical errors rod2D.setPos(rod2D.getPos() + trans_lab); rod2D.setE(e_new); @@ -138,27 +150,28 @@ void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) { Eigen::Matrix2d Integratoren2d_forceless::e_2_matrix(Eigen::Vector2d m_e) { Eigen::Matrix2d mat; - mat << m_e[0], -m_e[1], - m_e[1], m_e[0]; + mat << m_e[0], -m_e[1], m_e[1], m_e[0]; return mat; } -void Integratoren2d_forceless::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) { +void Integratoren2d_forceless::Set5_MBD(Rod2d & /*rod2D*/, + Simulation & /*sim*/) { /* auto std = sim.getSTD(); auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std)); - auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System - auto trans_lab = rod2D.getE_Base_matrix() * trans_part; - auto pos_predictor = rod2D.getPos() + trans_lab; + auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in + Particle System auto trans_lab = rod2D.getE_Base_matrix() * trans_part; auto + pos_predictor = rod2D.getPos() + trans_lab; //rotation - auto rot_predict = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt());// rotationsmatrix verwenden. - auto e = rod2D.getE(); - auto delta_e_predict = rot_predict * ortho(e); - auto e_predict = (e + delta_e_predict).normalized(); + auto rot_predict = Eigen::Rotation2D(sim.getNorm(std) * + rod2D.getDRot_Sqrt());// rotationsmatrix verwenden. auto e = rod2D.getE(); + auto delta_e_predict = rot_predict * ortho(e); auto e_predict = (e + + delta_e_predict).normalized(); - auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() * rod2D.getDiff(); //old + new diff + auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() * + rod2D.getDiff(); //old + new diff auto rot = Eigen::Rotation2D(sim.getNorm(std) * rod2D.getDRot_Sqrt()); auto e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict)); @@ -167,5 +180,3 @@ void Integratoren2d_forceless::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) { rod2D.setE(e_integrated.normalized()); */ } - - diff --git a/src/Integratoren2d_forceless.h b/src/Integratoren2d_forceless.h index 0033dd3..9ba8344 100644 --- a/src/Integratoren2d_forceless.h +++ b/src/Integratoren2d_forceless.h @@ -2,14 +2,13 @@ // Created by jholder on 21.10.21. // -#ifndef MYPROJECT_INTEGRATOREN2D_FORCELESS_H -#define MYPROJECT_INTEGRATOREN2D_FORCELESS_H +#pragma once #include "Rod2d.hpp" #include "Simulation.h" class Integratoren2d_forceless { -public: + public: static void Set1_Euler(Rod2d &rod2D, Simulation &sim); static void Set2_Heun(Rod2d &rod2D, Simulation &sim); @@ -26,9 +25,6 @@ public: static void Set0_deterministic(Rod2d &rod2D, Simulation &sim); -private: + private: static Eigen::Vector2d ortho(Eigen::Vector2d e); }; - - -#endif //MYPROJECT_INTEGRATOREN2D_FORCELESS_H diff --git a/src/LiveAgg.cpp b/src/LiveAgg.cpp index 90bf8a1..eda1b71 100644 --- a/src/LiveAgg.cpp +++ b/src/LiveAgg.cpp @@ -1,33 +1,21 @@ #include "LiveAgg.hpp" + #include -LiveAgg::LiveAgg() : num_of_data(0), mean(0.0), S(0.0) -{ +LiveAgg::LiveAgg() : num_of_data(0), mean(0.0), S(0.0) {} + +double LiveAgg::getSD() const noexcept { return S / (num_of_data - 1); } + +double LiveAgg::getSEM() const noexcept { + return S / (num_of_data - 1) / std::sqrt(num_of_data); } -double LiveAgg::getSD() const noexcept -{ - return S / (num_of_data - 1); -} +double LiveAgg::getMean() const noexcept { return mean; } -double LiveAgg::getSEM() const noexcept -{ - return S / (num_of_data - 1) / std::sqrt(num_of_data); -} - -double LiveAgg::getMean() const noexcept -{ - return mean; -} - -void LiveAgg::feed(double value) noexcept -{ - num_of_data += 1; - auto delta = value - mean; - mean += delta / num_of_data; - auto delta2 = value - mean; - S += delta * delta2; -} -int LiveAgg::getNumPoints() const noexcept -{ - return num_of_data; +void LiveAgg::feed(double value) noexcept { + num_of_data += 1; + auto delta = value - mean; + mean += delta / num_of_data; + auto delta2 = value - mean; + S += delta * delta2; } +int LiveAgg::getNumPoints() const noexcept { return num_of_data; } diff --git a/src/LiveAgg.hpp b/src/LiveAgg.hpp index e15271f..aa15eaf 100644 --- a/src/LiveAgg.hpp +++ b/src/LiveAgg.hpp @@ -1,20 +1,15 @@ -#ifndef LIVEAGG_H -#define LIVEAGG_H +#pragma once +class LiveAgg { + public: + LiveAgg(); + [[nodiscard]] double getSD() const noexcept; + [[nodiscard]] double getSEM() const noexcept; + [[nodiscard]] double getMean() const noexcept; + [[nodiscard]] int getNumPoints() const noexcept; + void feed(double value) noexcept; - -class LiveAgg -{ - public: - LiveAgg(); - [[nodiscard]] double getSD() const noexcept; - [[nodiscard]] double getSEM() const noexcept; - [[nodiscard]] double getMean() const noexcept; - [[nodiscard]] int getNumPoints() const noexcept; - void feed(double value) noexcept; - private: - int num_of_data; - double mean; - double S; + private: + int num_of_data; + double mean; + double S; }; - -#endif // LIVEAGG_H diff --git a/src/Rod2d.cpp b/src/Rod2d.cpp index 517de77..05e53e9 100644 --- a/src/Rod2d.cpp +++ b/src/Rod2d.cpp @@ -1,28 +1,30 @@ #include "Rod2d.hpp" + #include "Eigen/Dense" constexpr double M_Pl = 3.141592653589793238462643383279502884; /* pi */ -void sqrt(Eigen::Matrix2d &mat){ +void sqrt(Eigen::Matrix2d &mat) { const auto size = static_cast(mat.size()); for (size_t i = 0; i < size; i++) { mat.data()[i] = sqrt(mat.data()[i]); } } Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) { - assert(L>=1.0); + assert(L >= 1.0); if (L == 1.0) { m_D_rot = 3.0; - m_Diff << 1, 0, - 0, 1; + m_Diff << 1, 0, 0, 1; } else { const double D0 = 3.0 * M_Pl / L; const double p = L; - auto D_para = D0 / (M_Pl * 2.0) * (log(p) - 0.207 + 0.980 / p - 0.133 / (p * p)); - auto D_ortho = D0 / (M_Pl * 2.0) * (log(p) + 0.839 + 0.185 / p + 0.233 / (p * p)); - m_D_rot = 3 * D0 / (M_Pl * L * L) * (log(p) - 0.662 + 0.917 / p - 0.050 / (p * p)); - m_Diff << D_para, 0, - 0, D_ortho; + auto D_para = + D0 / (M_Pl * 2.0) * (log(p) - 0.207 + 0.980 / p - 0.133 / (p * p)); + auto D_ortho = + D0 / (M_Pl * 2.0) * (log(p) + 0.839 + 0.185 / p + 0.233 / (p * p)); + m_D_rot = 3 * D0 / (M_Pl * L * L) * + (log(p) - 0.662 + 0.917 / p - 0.050 / (p * p)); + m_Diff << D_para, 0, 0, D_ortho; } m_Diff_sqrt = m_Diff; m_D_rot_sqrt = sqrt(m_D_rot); @@ -34,43 +36,24 @@ void Rod2d::reset() { m_e = {1, 0}; } -void Rod2d::setPos(const Eigen::Vector2d &Pos) { - m_pos = Pos; -} +void Rod2d::setPos(const Eigen::Vector2d &Pos) { m_pos = Pos; } -double Rod2d::getDRot() const { - return m_D_rot; -} +double Rod2d::getDRot() const { return m_D_rot; } -const Eigen::Vector2d &Rod2d::getPos() const { - return m_pos; -} +const Eigen::Vector2d &Rod2d::getPos() const { return m_pos; } -const Eigen::Matrix2d &Rod2d::getDiff() const { - return m_Diff; -} +const Eigen::Matrix2d &Rod2d::getDiff() const { return m_Diff; } -const Eigen::Matrix2d &Rod2d::getDiff_Sqrt() const { - return m_Diff_sqrt; -} +const Eigen::Matrix2d &Rod2d::getDiff_Sqrt() const { return m_Diff_sqrt; } -double Rod2d::getDRot_Sqrt() const { - return m_D_rot_sqrt; -} +double Rod2d::getDRot_Sqrt() const { return m_D_rot_sqrt; } -const Eigen::Vector2d &Rod2d::getE() const { - return m_e; -} +const Eigen::Vector2d &Rod2d::getE() const { return m_e; } -void Rod2d::setE(const Eigen::Vector2d &mE) { - m_e = mE; -} +void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; } Eigen::Matrix2d Rod2d::getE_Base_matrix() const { Eigen::Matrix2d mat; - mat << m_e[0], -m_e[1], - m_e[1], m_e[0]; + mat << m_e[0], -m_e[1], m_e[1], m_e[0]; return mat; } - - diff --git a/src/Rod2d.hpp b/src/Rod2d.hpp index bc2beeb..f594b4c 100644 --- a/src/Rod2d.hpp +++ b/src/Rod2d.hpp @@ -1,9 +1,8 @@ #pragma once -#include #include -class Rod2d -{ -public: +#include +class Rod2d { + public: explicit Rod2d(double L); void reset(); @@ -25,13 +24,11 @@ public: Eigen::Matrix2d getE_Base_matrix() const; -private: + private: Eigen::Matrix2d m_Diff; Eigen::Matrix2d m_Diff_sqrt; double m_D_rot; double m_D_rot_sqrt; - Eigen::Vector2d m_pos; // position + Eigen::Vector2d m_pos; // position Eigen::Vector2d m_e; }; - - diff --git a/src/Simulation.cpp b/src/Simulation.cpp index 176e752..7d39a2d 100644 --- a/src/Simulation.cpp +++ b/src/Simulation.cpp @@ -5,17 +5,15 @@ #include "Simulation.h" double Simulation::getNorm(double t_norm) { - return t_norm*m_norm(m_generator); + return t_norm * m_norm(m_generator); } -Simulation::Simulation(double t_delta_T, size_t seed): m_delta_T(t_delta_T),m_std(std::sqrt(t_delta_T*2.0)), m_generator(seed), m_norm(0, 1.0) { +Simulation::Simulation(double t_delta_T, size_t seed) + : 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::getMDeltaT() const { return m_delta_T; } -double Simulation::getMDeltaT() const { - return m_delta_T; -} - -double Simulation::getSTD() const { - return m_std; -} +double Simulation::getSTD() const { return m_std; } diff --git a/src/Simulation.h b/src/Simulation.h index 8edc09a..aa0ed0d 100644 --- a/src/Simulation.h +++ b/src/Simulation.h @@ -2,25 +2,22 @@ // Created by jholder on 21.10.21. // -#ifndef MYPROJECT_SIMULATION_H -#define MYPROJECT_SIMULATION_H - +#pragma once #include class Simulation { -public: + public: explicit Simulation(double t_delta_T, size_t seed); double getNorm(double t_norm); -private: + + private: double m_delta_T; double m_std; std::mt19937_64 m_generator; std::normal_distribution m_norm; -public: + + public: [[nodiscard]] double getMDeltaT() const; [[nodiscard]] double getSTD() const; }; - - -#endif //MYPROJECT_SIMULATION_H diff --git a/src/main.cpp b/src/main.cpp index 8984f7d..df9b018 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -3,49 +3,73 @@ // #include + #include "Calculation.h" -#include "Integratoren2d_forceless.h" #include "Integratoren2d_force.h" +#include "Integratoren2d_forceless.h" int main() { - constexpr int numStep = 1000; - auto zero_Force = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; }; - auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> double { return 0.0; }; + constexpr int numStep = 1000000; + constexpr double k = 0.01; + [[maybe_unused]] auto harmonic_Force = + [](const Eigen::Vector2d &pos, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return -k * pos; + }; + [[maybe_unused]] auto zero_Force = + [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return {0.0, 0.0}; + }; + auto zero_Torque = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> double { + return 0.0; + }; /* auto euler_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set1_Euler(rod, sim, zero_Force, zero_Torque); + return Integratoren2d_force::Set1_Euler(rod, sim, zero_Force, + zero_Torque); };*/ auto heun_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set2_Heun(rod, sim, zero_Force, zero_Torque); + return Integratoren2d_force::Set1_Euler(rod, sim, harmonic_Force, + zero_Torque); }; /*auto exact_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set3_Exact(rod, sim, zero_Force, zero_Torque); - };*//* - auto bdas_Zero = [&](Rod2d &rod, Simulation &sim) { - return Integratoren2d_force::Set4_BDAS(rod, sim, zero_Force, zero_Torque); + return Integratoren2d_force::Set3_Exact(rod, sim, zero_Force, + zero_Torque); };*/ - Calculation euler(heun_Zero, {{Compute::Type::msd, 1}, - {Compute::Type::msd, 2}, - {Compute::Type::msd, 4}, - {Compute::Type::msd, 8}, - {Compute::Type::msd, 16}, - {Compute::Type::oaf, 1}, - {Compute::Type::oaf, 2}, - {Compute::Type::oaf, 4}, - {Compute::Type::oaf, 8}, - {Compute::Type::oaf, 16}}, 0.01, 12345); + /* + auto bdas_Zero = [&](Rod2d &rod, Simulation &sim) { + return Integratoren2d_force::Set4_BDAS(rod, sim, zero_Force, zero_Torque); + };*/ + Calculation euler(heun_Zero, + {{Compute::Type::msd, 1}, + {Compute::Type::msd, 2}, + {Compute::Type::msd, 4}, + {Compute::Type::msd, 8}, + {Compute::Type::msd, 16}, + {Compute::Type::msd, 32}, + {Compute::Type::msd, 64}, + {Compute::Type::msd, 128}, + {Compute::Type::msd, 256}, + {Compute::Type::msd, 512}, + {Compute::Type::msd, 1024}, + {Compute::Type::oaf, 1}, + {Compute::Type::oaf, 2}, + {Compute::Type::oaf, 4}, + {Compute::Type::oaf, 8}, + {Compute::Type::oaf, 16}}, + 0.01, 12345); euler.run(numStep); - for (const auto &com: euler.getComputes()) { - if (com.getType() != Compute::Type::msd) - continue; - std::cout << "MSD: " << com.getDifference() << " " << com.getAgg().getMean() - << " <-> " << com.getTarget() << std::endl; + for (const auto &com : euler.getComputes()) { + if (com.getType() != Compute::Type::msd) continue; + std::cout << "MSD: " << com.getDifference() << " " + << com.getAgg().getMean() << " <-> " << com.getTarget() + << std::endl; } - for (const auto &com: euler.getComputes()) { - if (com.getType() != Compute::Type::oaf) - continue; - std::cout << "OAF: " << com.getDifference() << " " << com.getAgg().getMean() - << " <-> "<< com.getTarget() << std::endl; + for (const auto &com : euler.getComputes()) { + if (com.getType() != Compute::Type::oaf) continue; + std::cout << "OAF: " << com.getDifference() << " " + << com.getAgg().getMean() << " <-> " << com.getTarget() + << std::endl; } } \ No newline at end of file diff --git a/test/catch_main.cpp b/test/catch_main.cpp index d8d2eca..db3af94 100644 --- a/test/catch_main.cpp +++ b/test/catch_main.cpp @@ -1,5 +1,3 @@ -#define CATCH_CONFIG_MAIN // This tells the catch header to generate a main +#define CATCH_CONFIG_MAIN // This tells the catch header to generate a main #include - - diff --git a/test/test_Calculation.cpp b/test/test_Calculation.cpp index 5679651..420feca 100644 --- a/test/test_Calculation.cpp +++ b/test/test_Calculation.cpp @@ -2,32 +2,41 @@ // Created by jholder on 22.10.21. // +#include "Calculation.h" + #include -#include -#include "LiveAgg.hpp" + #include "Integratoren2d_forceless.h" +#include "LiveAgg.hpp" TEST_CASE("Basic run of Calculation") { - Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); - Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); - Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); - Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 10}, - {Compute::Type::oaf, 10}}, 1, 1); + Calculation euler(Integratoren2d_forceless::Set1_Euler, + {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, + 1); + Calculation heun(Integratoren2d_forceless::Set2_Heun, + {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, + 1); + Calculation exact(Integratoren2d_forceless::Set3_Exact, + {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, + 1); + Calculation bdas(Integratoren2d_forceless::Set4_BDAS, + {{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1, + 1); SECTION("Euler") { euler.run(100); CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0)); CHECK(euler.getComputes()[0].getType() == Compute::Type::msd); CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf); - }SECTION("Heun") { + } + SECTION("Heun") { heun.run(100); CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - }SECTION("Exact") { + } + SECTION("Exact") { exact.run(100); CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0)); - }SECTION("Euler") { + } + SECTION("Euler") { bdas.run(100); CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0)); } @@ -36,17 +45,20 @@ TEST_CASE("Basic run of Calculation") { euler.run(100); CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10); CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10); - }SECTION("Deterministic") { + } + SECTION("Deterministic") { Calculation determ(Integratoren2d_forceless::Set0_deterministic, - {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 1, 1); + {{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, + 1, 1); determ.run(10); auto time = 1; auto targetMSD = pow(0.01, 2) * time; auto targetOAF = cos(0.1); - CHECK(determ.getComputes()[0].getAgg().getMean() == Catch::Detail::Approx(targetMSD)); - CHECK(determ.getComputes()[1].getAgg().getMean() == Catch::Detail::Approx(targetOAF)); + CHECK(determ.getComputes()[0].getAgg().getMean() == + Catch::Detail::Approx(targetMSD)); + CHECK(determ.getComputes()[1].getAgg().getMean() == + Catch::Detail::Approx(targetOAF)); } } diff --git a/test/test_Compute.cpp b/test/test_Compute.cpp index fbcceb3..6e06456 100644 --- a/test/test_Compute.cpp +++ b/test/test_Compute.cpp @@ -2,19 +2,20 @@ // Created by jholder on 22.10.21. // -#include -#include "catch2/catch.hpp" +#include "Compute.h" + +#include #include "Integratoren2d_forceless.h" TEST_CASE("Compute") { Rod2d rod(1.0); - Simulation sim(0.1,1); + Simulation sim(0.1, 1); auto com = Compute(rod, Compute::Type::msd, 10, sim); SECTION("Mean of same values") { for (int i = 0; i < 100; ++i) { com.eval(rod); } - CHECK(com.getAgg().getNumPoints()==10); - CHECK(com.getAgg().getMean()==0.0); + CHECK(com.getAgg().getNumPoints() == 10); + CHECK(com.getAgg().getMean() == 0.0); } } diff --git a/test/test_Integratoren.cpp b/test/test_Integratoren.cpp index 62bebd6..a618bd8 100644 --- a/test/test_Integratoren.cpp +++ b/test/test_Integratoren.cpp @@ -1,103 +1,146 @@ // // Created by jholder on 24.10.21. // +#include "Calculation.h" +#include "Compute.h" +#include "Integratoren2d_force.h" +#include "Integratoren2d_forceless.h" + #include -#include -#include -#include -#include TEST_CASE("Forceless Integratoren") { const size_t SEED = Catch::rngSeed(); - auto force = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; }; - auto torque = [](const Eigen::Vector2d & /*pos*/, - const Eigen::Vector2d & /*torque*/) -> double { return 0.0; }; + auto force = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { + return {0.0, 0.0}; + }; + auto torque = [](const Eigen::Vector2d & /*pos*/, + const Eigen::Vector2d & /*torque*/) -> double { + return 0.0; + }; SECTION("Euler") { - - Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED); + Calculation euler(Integratoren2d_forceless::Set1_Euler, + {{Compute::Type::msd, 1}, + {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); } + for (auto &c : euler.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } } - Calculation euler_force(Integratoren2d_force::Set1_Euler, { - {Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}}, 0.01, SEED, force, torque); + Calculation euler_force( + Integratoren2d_force::Set1_Euler, + {{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, 0.01, SEED, + force, torque); { euler_force.run(10000); - for (auto &c: euler_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + for (auto &c : euler_force.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } } - CHECK(euler.getComputes()[0].getAgg().getMean() == euler_force.getComputes()[0].getAgg().getMean()); - CHECK(euler.getComputes()[1].getAgg().getMean() == euler_force.getComputes()[1].getAgg().getMean()); - }SECTION("Heun") { - - Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED); + CHECK(euler.getComputes()[0].getAgg().getMean() == + euler_force.getComputes()[0].getAgg().getMean()); + CHECK(euler.getComputes()[1].getAgg().getMean() == + euler_force.getComputes()[1].getAgg().getMean()); + } + SECTION("Heun") { + Calculation heun(Integratoren2d_forceless::Set2_Heun, + {{Compute::Type::msd, 1}, + {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); } + 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}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED, force, - torque); + Calculation heun_force(Integratoren2d_force::Set2_Heun, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED, force, torque); { heun_force.run(10000); - for (auto &c: heun_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + for (auto &c : heun_force.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } } - CHECK(heun.getComputes()[0].getAgg().getMean() == heun_force.getComputes()[0].getAgg().getMean()); - CHECK(heun.getComputes()[1].getAgg().getMean() == heun_force.getComputes()[1].getAgg().getMean()); - }SECTION("Exact") { - - Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED); + CHECK(heun.getComputes()[0].getAgg().getMean() == + heun_force.getComputes()[0].getAgg().getMean()); + CHECK(heun.getComputes()[1].getAgg().getMean() == + heun_force.getComputes()[1].getAgg().getMean()); + } + SECTION("Exact") { + Calculation exact(Integratoren2d_forceless::Set3_Exact, + {{Compute::Type::msd, 1}, + {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); } + 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}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED, force, - torque); + Calculation exact_force(Integratoren2d_force::Set3_Exact, + {{Compute::Type::msd, 1}, + {Compute::Type::oaf, 1}, + {Compute::Type::empxx, 1}, + {Compute::Type::empyy, 1}}, + 0.01, SEED, force, torque); { exact_force.run(10000); - for (auto &c: exact_force.getComputes()) { CHECK(c.getDifference() <= 0.005); } + for (auto &c : exact_force.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } } - CHECK(exact.getComputes()[0].getAgg().getMean() == Approx(exact_force.getComputes()[0].getAgg().getMean()).epsilon(1e-10)); - CHECK(exact.getComputes()[1].getAgg().getMean() == exact_force.getComputes()[1].getAgg().getMean()); - }SECTION("BDAS") { - - Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1}, - {Compute::Type::oaf, 1}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED); + CHECK(exact.getComputes()[0].getAgg().getMean() == + Approx(exact_force.getComputes()[0].getAgg().getMean()) + .epsilon(1e-10)); + CHECK(exact.getComputes()[1].getAgg().getMean() == + exact_force.getComputes()[1].getAgg().getMean()); + } + SECTION("BDAS") { + Calculation bdas(Integratoren2d_forceless::Set4_BDAS, + {{Compute::Type::msd, 1}, + {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); } + 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}, - {Compute::Type::empxx, 1}, - {Compute::Type::empyy, 1}}, 0.01, SEED, force, - torque); + Calculation bdas_force(Integratoren2d_force::Set4_BDAS, + {{Compute::Type::msd, 1}, + {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); } + for (auto &c : bdas_force.getComputes()) { + CHECK(c.getDifference() <= 0.005); + } } - 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)); + 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_LiveAgg.cpp b/test/test_LiveAgg.cpp index 60c5021..2246153 100644 --- a/test/test_LiveAgg.cpp +++ b/test/test_LiveAgg.cpp @@ -3,14 +3,15 @@ // #include + #include "LiveAgg.hpp" TEST_CASE("LiveAgg") { auto ag = LiveAgg(); - SECTION("Mean of same values"){ + SECTION("Mean of same values") { ag.feed(1.0); ag.feed(1.0); - REQUIRE(ag.getMean()==1.0); - REQUIRE(ag.getNumPoints()==2); + REQUIRE(ag.getMean() == 1.0); + REQUIRE(ag.getNumPoints() == 2); } } diff --git a/test/test_Rod2d.cpp b/test/test_Rod2d.cpp index d233c41..3cbdffd 100644 --- a/test/test_Rod2d.cpp +++ b/test/test_Rod2d.cpp @@ -2,9 +2,10 @@ // Created by jholder on 21.10.21. // -#include "Rod2d.hpp" -#include #include +#include + +#include "Rod2d.hpp" TEST_CASE("Rods") { Rod2d sphere(1); @@ -14,45 +15,47 @@ TEST_CASE("Rods") { auto newPos = rod.getPos(); 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") { - SECTION("Forwards == e 90"){ + 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") { + SECTION("Forwards == e 90") { // // ===== // - rod.setE({0, 1}); // Lab system is 90 degree rotated - auto test = Eigen::Vector2d ({1,0}); - auto erg = (rod.getE_Base_matrix()*test).normalized(); + rod.setE({0, 1}); // Lab system is 90 degree rotated + auto test = Eigen::Vector2d({1, 0}); + auto erg = (rod.getE_Base_matrix() * test).normalized(); auto e = rod.getE(); REQUIRE(erg.isApprox(e)); - } - SECTION("Forwards == e 45"){ + SECTION("Forwards == e 45") { // o // o // o - rod.setE(Eigen::Vector2d ({1,1}).normalized()); // Lab system is 90 degree rotated - auto test = Eigen::Vector2d ({1,0}); - auto erg = (rod.getE_Base_matrix()*test).normalized(); + rod.setE(Eigen::Vector2d({1, 1}) + .normalized()); // Lab system is 90 degree rotated + auto test = Eigen::Vector2d({1, 0}); + auto erg = (rod.getE_Base_matrix() * test).normalized(); auto e = rod.getE(); REQUIRE(erg.isApprox(e)); } - - } - SECTION("Angle Stuff"){ + 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()); + 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)); + 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 diff --git a/test/test_Simulation.cpp b/test/test_Simulation.cpp index e750a7d..97f19a7 100644 --- a/test/test_Simulation.cpp +++ b/test/test_Simulation.cpp @@ -3,25 +3,26 @@ // #include + #include "Simulation.h" TEST_CASE("Simulation") { - auto sim = Simulation(1,0); - REQUIRE(sim.getSTD()==sqrt(1*2.0)); - SECTION("Symmetrie"){ + auto sim = Simulation(1, 0); + REQUIRE(sim.getSTD() == sqrt(1 * 2.0)); + SECTION("Symmetrie") { constexpr int num = 10000000; double sum = 0.0; for (int i = 0; i < num; ++i) { sum += sim.getNorm(1); } - CHECK(sum/num == Catch::Detail::Approx(0.0).margin(0.001)); + CHECK(sum / num == Catch::Detail::Approx(0.0).margin(0.001)); } - SECTION("STD"){ + SECTION("STD") { constexpr int num = 10000000; double sum = 0.0; for (int i = 0; i < num; ++i) { - sum += pow(sim.getNorm(1),2); + sum += pow(sim.getNorm(1), 2); } - CHECK(sum/num == Catch::Detail::Approx(1.0).margin(0.001)); + CHECK(sum / num == Catch::Detail::Approx(1.0).margin(0.001)); } } \ No newline at end of file