Compare commits
4 Commits
81f4796eeb
...
8af7eee267
Author | SHA1 | Date | |
---|---|---|---|
8af7eee267 | |||
362ec38dd2 | |||
07b1f194a6 | |||
7eac132a5f |
146
.clang-format
146
.clang-format
@ -1,12 +1,15 @@
|
|||||||
AccessModifierOffset: -2
|
---
|
||||||
AlignAfterOpenBracket: DontAlign
|
Language: Cpp
|
||||||
|
# BasedOnStyle: Google
|
||||||
|
AccessModifierOffset: -1
|
||||||
|
AlignAfterOpenBracket: Align
|
||||||
AlignConsecutiveAssignments: false
|
AlignConsecutiveAssignments: false
|
||||||
AlignConsecutiveDeclarations: false
|
AlignConsecutiveDeclarations: false
|
||||||
AlignEscapedNewlines: Left
|
AlignEscapedNewlines: Left
|
||||||
AlignOperands: true
|
AlignOperands: true
|
||||||
AlignTrailingComments: false
|
AlignTrailingComments: true
|
||||||
AllowAllParametersOfDeclarationOnNextLine: false
|
AllowAllParametersOfDeclarationOnNextLine: true
|
||||||
AllowShortBlocksOnASingleLine: true
|
AllowShortBlocksOnASingleLine: false
|
||||||
AllowShortCaseLabelsOnASingleLine: false
|
AllowShortCaseLabelsOnASingleLine: false
|
||||||
AllowShortFunctionsOnASingleLine: All
|
AllowShortFunctionsOnASingleLine: All
|
||||||
AllowShortIfStatementsOnASingleLine: true
|
AllowShortIfStatementsOnASingleLine: true
|
||||||
@ -14,85 +17,92 @@ AllowShortLoopsOnASingleLine: true
|
|||||||
AlwaysBreakAfterDefinitionReturnType: None
|
AlwaysBreakAfterDefinitionReturnType: None
|
||||||
AlwaysBreakAfterReturnType: None
|
AlwaysBreakAfterReturnType: None
|
||||||
AlwaysBreakBeforeMultilineStrings: true
|
AlwaysBreakBeforeMultilineStrings: true
|
||||||
AlwaysBreakTemplateDeclarations: false
|
AlwaysBreakTemplateDeclarations: true
|
||||||
BinPackArguments: false
|
BinPackArguments: true
|
||||||
BinPackParameters: false
|
BinPackParameters: true
|
||||||
BraceWrapping:
|
BraceWrapping:
|
||||||
AfterClass: true
|
AfterClass: false
|
||||||
AfterControlStatement: false
|
AfterControlStatement: false
|
||||||
AfterEnum: false
|
AfterEnum: false
|
||||||
AfterFunction: true
|
AfterFunction: false
|
||||||
AfterNamespace: false
|
AfterNamespace: false
|
||||||
AfterObjCDeclaration: false
|
AfterObjCDeclaration: false
|
||||||
AfterStruct: true
|
AfterStruct: false
|
||||||
AfterUnion: false
|
AfterUnion: false
|
||||||
BeforeCatch: false
|
BeforeCatch: false
|
||||||
BeforeElse: false
|
BeforeElse: false
|
||||||
IndentBraces: false
|
IndentBraces: false
|
||||||
SplitEmptyFunction: false
|
SplitEmptyFunction: true
|
||||||
SplitEmptyNamespace: true
|
|
||||||
SplitEmptyRecord: true
|
SplitEmptyRecord: true
|
||||||
BreakAfterJavaFieldAnnotations: true
|
SplitEmptyNamespace: true
|
||||||
BreakBeforeBinaryOperators: NonAssignment
|
BreakBeforeBinaryOperators: None
|
||||||
BreakBeforeBraces: Custom
|
BreakBeforeBraces: Attach
|
||||||
BreakBeforeInheritanceComma: true
|
BreakBeforeInheritanceComma: false
|
||||||
BreakBeforeTernaryOperators: true
|
BreakBeforeTernaryOperators: true
|
||||||
BreakConstructorInitializers: BeforeColon
|
|
||||||
BreakConstructorInitializersBeforeComma: false
|
BreakConstructorInitializersBeforeComma: false
|
||||||
|
BreakConstructorInitializers: BeforeColon
|
||||||
|
BreakAfterJavaFieldAnnotations: false
|
||||||
BreakStringLiterals: true
|
BreakStringLiterals: true
|
||||||
ColumnLimit: 0
|
ColumnLimit: 80
|
||||||
CommentPragmas: '^ IWYU pragma:'
|
CommentPragmas: '^ IWYU pragma:'
|
||||||
CompactNamespaces: false
|
CompactNamespaces: false
|
||||||
ConstructorInitializerAllOnOneLineOrOnePerLine: false
|
ConstructorInitializerAllOnOneLineOrOnePerLine: true
|
||||||
ConstructorInitializerIndentWidth: 2
|
ConstructorInitializerIndentWidth: 4
|
||||||
ContinuationIndentWidth: 2
|
ContinuationIndentWidth: 4
|
||||||
Cpp11BracedListStyle: false
|
Cpp11BracedListStyle: true
|
||||||
DerivePointerAlignment: false
|
DerivePointerAlignment: true
|
||||||
DisableFormat: false
|
DisableFormat: false
|
||||||
ExperimentalAutoDetectBinPacking: true
|
ExperimentalAutoDetectBinPacking: false
|
||||||
FixNamespaceComments: true
|
FixNamespaceComments: true
|
||||||
ForEachMacros:
|
ForEachMacros:
|
||||||
- foreach
|
- foreach
|
||||||
- Q_FOREACH
|
- Q_FOREACH
|
||||||
- BOOST_FOREACH
|
- BOOST_FOREACH
|
||||||
IncludeCategories:
|
IncludeCategories:
|
||||||
- Priority: 2
|
- Regex: '^<.*\.h>'
|
||||||
Regex: ^"(llvm|llvm-c|clang|clang-c)/
|
Priority: 1
|
||||||
- Priority: 3
|
- Regex: '^<.*'
|
||||||
Regex: ^(<|"(gtest|gmock|isl|json)/)
|
Priority: 2
|
||||||
- Priority: 1
|
- Regex: '.*'
|
||||||
Regex: .*
|
Priority: 3
|
||||||
IncludeIsMainRegex: (Test)?$
|
IncludeIsMainRegex: '([-_](test|unittest))?$'
|
||||||
IndentCaseLabels: false
|
IndentCaseLabels: true
|
||||||
IndentWidth: 2
|
IndentWidth: 4
|
||||||
IndentWrappedFunctionNames: true
|
IndentWrappedFunctionNames: false
|
||||||
JavaScriptQuotes: Leave
|
JavaScriptQuotes: Leave
|
||||||
JavaScriptWrapImports: true
|
JavaScriptWrapImports: true
|
||||||
KeepEmptyLinesAtTheStartOfBlocks: true
|
KeepEmptyLinesAtTheStartOfBlocks: false
|
||||||
Language: Cpp
|
|
||||||
MacroBlockBegin: ''
|
MacroBlockBegin: ''
|
||||||
MacroBlockEnd: ''
|
MacroBlockEnd: ''
|
||||||
MaxEmptyLinesToKeep: 2
|
MaxEmptyLinesToKeep: 1
|
||||||
NamespaceIndentation: Inner
|
NamespaceIndentation: None
|
||||||
ObjCBlockIndentWidth: 7
|
ObjCBlockIndentWidth: 4
|
||||||
ObjCSpaceAfterProperty: true
|
ObjCSpaceAfterProperty: false
|
||||||
ObjCSpaceBeforeProtocolList: false
|
ObjCSpaceBeforeProtocolList: false
|
||||||
PointerAlignment: Right
|
PenaltyBreakAssignment: 2
|
||||||
ReflowComments: true
|
PenaltyBreakBeforeFirstCallParameter: 1
|
||||||
SortIncludes: false
|
PenaltyBreakComment: 300
|
||||||
SortUsingDeclarations: false
|
PenaltyBreakFirstLessLess: 120
|
||||||
|
PenaltyBreakString: 1000
|
||||||
|
PenaltyExcessCharacter: 1000000
|
||||||
|
PenaltyReturnTypeOnItsOwnLine: 200
|
||||||
|
PointerAlignment: Left
|
||||||
|
ReflowComments: true
|
||||||
|
SortIncludes: true
|
||||||
|
SortUsingDeclarations: true
|
||||||
SpaceAfterCStyleCast: false
|
SpaceAfterCStyleCast: false
|
||||||
SpaceAfterTemplateKeyword: false
|
SpaceAfterTemplateKeyword: true
|
||||||
SpaceBeforeAssignmentOperators: true
|
SpaceBeforeAssignmentOperators: true
|
||||||
SpaceBeforeParens: ControlStatements
|
SpaceBeforeParens: ControlStatements
|
||||||
SpaceInEmptyParentheses: false
|
SpaceInEmptyParentheses: false
|
||||||
SpacesBeforeTrailingComments: 0
|
SpacesBeforeTrailingComments: 2
|
||||||
SpacesInAngles: false
|
SpacesInAngles: false
|
||||||
SpacesInCStyleCastParentheses: false
|
|
||||||
SpacesInContainerLiterals: true
|
SpacesInContainerLiterals: true
|
||||||
|
SpacesInCStyleCastParentheses: false
|
||||||
SpacesInParentheses: false
|
SpacesInParentheses: false
|
||||||
SpacesInSquareBrackets: false
|
SpacesInSquareBrackets: false
|
||||||
Standard: Cpp11
|
Standard: Auto
|
||||||
TabWidth: 8
|
TabWidth: 8
|
||||||
UseTab: Never
|
UseTab: Never
|
||||||
|
...
|
||||||
|
|
||||||
|
36
benchmark.81f4796.fourier
Normal file
36
benchmark.81f4796.fourier
Normal file
@ -0,0 +1,36 @@
|
|||||||
|
<?xml version="1.0" encoding="UTF-8"?>
|
||||||
|
<Catch name="benchmark">
|
||||||
|
<Group name="benchmark">
|
||||||
|
<TestCase name="Euler - Baseline" tags="[benchmark]" filename="/home/watt/jholder/Integratoren/benchmark/bench.cpp" line="12">
|
||||||
|
<BenchmarkResults name="Euler without force" samples="100" resamples="100000" iterations="179" clockResolution="19.6142" estimatedDuration="1.9511e+06">
|
||||||
|
<!--All values in nano seconds-->
|
||||||
|
<mean value="109.097" lowerBound="108.49" upperBound="109.829" ci="0.95"/>
|
||||||
|
<standardDeviation value="3.40566" lowerBound="2.90384" upperBound="3.93169" ci="0.95"/>
|
||||||
|
<outliers variance="0.267429" lowMild="0" lowSevere="0" highMild="15" highSevere="3"/>
|
||||||
|
</BenchmarkResults>
|
||||||
|
<BenchmarkResults name="Heun without force" samples="100" resamples="100000" iterations="128" clockResolution="19.6142" estimatedDuration="1.9584e+06">
|
||||||
|
<!--All values in nano seconds-->
|
||||||
|
<mean value="152.946" lowerBound="152.282" upperBound="153.785" ci="0.95"/>
|
||||||
|
<standardDeviation value="3.77535" lowerBound="3.05802" upperBound="4.579" ci="0.95"/>
|
||||||
|
<outliers variance="0.180743" lowMild="0" lowSevere="0" highMild="8" highSevere="3"/>
|
||||||
|
</BenchmarkResults>
|
||||||
|
<BenchmarkResults name="Exact without force" samples="100" resamples="100000" iterations="166" clockResolution="19.6142" estimatedDuration="1.9588e+06">
|
||||||
|
<!--All values in nano seconds-->
|
||||||
|
<mean value="118.824" lowerBound="117.587" upperBound="121.925" ci="0.95"/>
|
||||||
|
<standardDeviation value="9.21278" lowerBound="2.96214" upperBound="17.3383" ci="0.95"/>
|
||||||
|
<outliers variance="0.696897" lowMild="0" lowSevere="0" highMild="0" highSevere="7"/>
|
||||||
|
</BenchmarkResults>
|
||||||
|
<BenchmarkResults name="BDAS without force" samples="100" resamples="100000" iterations="114" clockResolution="19.6142" estimatedDuration="1.9722e+06">
|
||||||
|
<!--All values in nano seconds-->
|
||||||
|
<mean value="173.079" lowerBound="171.528" upperBound="176.479" ci="0.95"/>
|
||||||
|
<standardDeviation value="11.1774" lowerBound="6.53959" upperBound="22.4088" ci="0.95"/>
|
||||||
|
<outliers variance="0.605598" lowMild="0" lowSevere="0" highMild="0" highSevere="1"/>
|
||||||
|
</BenchmarkResults>
|
||||||
|
<OverallResult success="true"/>
|
||||||
|
</TestCase>
|
||||||
|
<OverallResults successes="0" failures="0" expectedFailures="0"/>
|
||||||
|
<OverallResultsCases successes="1" failures="0" expectedFailures="0"/>
|
||||||
|
</Group>
|
||||||
|
<OverallResults successes="0" failures="0" expectedFailures="0"/>
|
||||||
|
<OverallResultsCases successes="1" failures="0" expectedFailures="0"/>
|
||||||
|
</Catch>
|
@ -6,27 +6,25 @@
|
|||||||
#define CATCH_CONFIG_ENABLE_BENCHMARKING
|
#define CATCH_CONFIG_ENABLE_BENCHMARKING
|
||||||
|
|
||||||
#include <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
#include "Rod2d.hpp"
|
|
||||||
#include "Integratoren2d_forceless.h"
|
#include "Integratoren2d_forceless.h"
|
||||||
|
#include "Rod2d.hpp"
|
||||||
|
|
||||||
TEST_CASE("Euler - Baseline", "[benchmark]")
|
TEST_CASE("Euler - Baseline", "[benchmark]") {
|
||||||
{
|
|
||||||
Rod2d rod(1.0);
|
Rod2d rod(1.0);
|
||||||
Simulation sim(0.01, Catch::rngSeed());
|
Simulation sim(0.01, Catch::rngSeed());
|
||||||
BENCHMARK("Euler without force")
|
BENCHMARK("Euler without force") {
|
||||||
{
|
Integratoren2d_forceless::Set1_Euler(rod, sim);
|
||||||
Integratoren2d_forceless::Set1_Euler(rod, sim);
|
};
|
||||||
};
|
BENCHMARK("Heun without force") {
|
||||||
BENCHMARK("Heun without force")
|
Integratoren2d_forceless::Set2_Heun(rod, sim);
|
||||||
{
|
};
|
||||||
Integratoren2d_forceless::Set2_Heun(rod, sim);
|
BENCHMARK("Exact without force") {
|
||||||
};
|
Integratoren2d_forceless::Set3_Exact(rod, sim);
|
||||||
BENCHMARK("Exact without force")
|
};
|
||||||
{
|
BENCHMARK("BDAS without force") {
|
||||||
Integratoren2d_forceless::Set3_Exact(rod, sim);
|
Integratoren2d_forceless::Set4_BDAS(rod, sim);
|
||||||
};
|
};
|
||||||
BENCHMARK("BDAS without force")
|
BENCHMARK("MBD without force") {
|
||||||
{
|
Integratoren2d_forceless::Set5_MBD(rod, sim);
|
||||||
Integratoren2d_forceless::Set4_BDAS(rod, sim);
|
};
|
||||||
};
|
|
||||||
}
|
}
|
@ -6,46 +6,48 @@
|
|||||||
|
|
||||||
#include <utility>
|
#include <utility>
|
||||||
|
|
||||||
|
|
||||||
void Calculation::run(size_t steps) {
|
void Calculation::run(size_t steps) {
|
||||||
for (size_t step = 0; step < steps; ++step) {
|
for (size_t step = 0; step < steps; ++step) {
|
||||||
m_integrator(rod, sim);
|
m_integrator(rod, sim);
|
||||||
for (auto &comp: computes) { comp.eval(rod); }
|
for (auto &comp : computes) {
|
||||||
|
comp.eval(rod);
|
||||||
|
}
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
const Rod2d &Calculation::getRod() const {
|
const Rod2d &Calculation::getRod() const { return rod; }
|
||||||
return rod;
|
|
||||||
}
|
|
||||||
|
|
||||||
Calculation::Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator,
|
Calculation::Calculation(
|
||||||
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes, double deltaT,
|
std::function<void(Rod2d &, Simulation &)> t_integrator,
|
||||||
size_t seed)
|
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
|
||||||
: sim(deltaT, seed), m_integrator(
|
double deltaT, size_t seed,double length)
|
||||||
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));
|
computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
Calculation::Calculation(
|
Calculation::Calculation(
|
||||||
std::function<void(Rod2d &, Simulation &, std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>,
|
std::function<
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)>)> t_integrator,
|
void(Rod2d &, Simulation &,
|
||||||
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes, double deltaT,
|
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>,
|
||||||
size_t seed, std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
std::function<double(Eigen::Vector2d, Eigen::Vector2d)>)>
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque)
|
t_integrator,
|
||||||
: sim(deltaT, seed) {
|
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
|
||||||
m_integrator = [&](Rod2d &t_rod, Simulation &t_sim) { t_integrator(t_rod, t_sim, force, torque); };
|
double deltaT, size_t seed,
|
||||||
for (const auto &pair: t_computes) {
|
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
||||||
|
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) {
|
||||||
computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
|
computes.emplace_back(Compute(rod, pair.first, pair.second, sim));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
const std::vector<Compute> &Calculation::getComputes() const {
|
const std::vector<Compute> &Calculation::getComputes() const {
|
||||||
return computes;
|
return computes;
|
||||||
}
|
}
|
||||||
|
|
||||||
const Simulation &Calculation::getSim() const {
|
const Simulation &Calculation::getSim() const { return sim; }
|
||||||
return sim;
|
|
||||||
}
|
|
||||||
|
@ -2,35 +2,37 @@
|
|||||||
// Created by jholder on 22.10.21.
|
// Created by jholder on 22.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
#ifndef MYPROJECT_CALCULATION_H
|
#pragma once
|
||||||
#define MYPROJECT_CALCULATION_H
|
|
||||||
|
|
||||||
|
#include <functional>
|
||||||
|
#include "Compute.h"
|
||||||
#include "Rod2d.hpp"
|
#include "Rod2d.hpp"
|
||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
#include "functional"
|
|
||||||
#include "Compute.h"
|
|
||||||
|
|
||||||
|
|
||||||
class Calculation {
|
class Calculation {
|
||||||
private:
|
private:
|
||||||
|
Rod2d rod;
|
||||||
Rod2d rod = Rod2d(1.0);
|
|
||||||
Simulation sim;
|
Simulation sim;
|
||||||
std::function<void(Rod2d &, Simulation &)> m_integrator;
|
std::function<void(Rod2d &, Simulation &)> m_integrator;
|
||||||
std::vector<Compute> computes;
|
std::vector<Compute> computes;
|
||||||
public:
|
|
||||||
const Simulation &getSim() const;
|
|
||||||
|
|
||||||
|
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(
|
Calculation(
|
||||||
std::function<void(Rod2d &, Simulation &, std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>,
|
inte_force_type
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)>)> t_integrator,
|
t_integrator,
|
||||||
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes, double deltaT, size_t seed,
|
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
|
||||||
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
double deltaT, size_t seed,
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque);
|
force_type force,
|
||||||
|
torque_type torque, double length = 1.0);
|
||||||
|
|
||||||
explicit Calculation(std::function<void(Rod2d &, Simulation &)> t_integrator,
|
explicit Calculation(
|
||||||
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes, double deltaT,
|
std::function<void(Rod2d &, Simulation &)> t_integrator,
|
||||||
size_t seed);
|
std::initializer_list<std::pair<Compute::Type, size_t>> t_computes,
|
||||||
|
double deltaT, size_t seed, double length = 1.0);
|
||||||
|
|
||||||
[[nodiscard]] const std::vector<Compute> &getComputes() const;
|
[[nodiscard]] const std::vector<Compute> &getComputes() const;
|
||||||
|
|
||||||
@ -38,6 +40,3 @@ public:
|
|||||||
|
|
||||||
[[nodiscard]] const Rod2d &getRod() const;
|
[[nodiscard]] const Rod2d &getRod() const;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
#endif //MYPROJECT_CALCULATION_H
|
|
||||||
|
@ -4,8 +4,9 @@
|
|||||||
|
|
||||||
#include "Compute.h"
|
#include "Compute.h"
|
||||||
|
|
||||||
#include <utility>
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
#include <utility>
|
||||||
|
|
||||||
#include "Rod2d.hpp"
|
#include "Rod2d.hpp"
|
||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
@ -23,12 +24,21 @@ void Compute::evalOAF(const Rod2d &rod2D) {
|
|||||||
agg.feed(oaf);
|
agg.feed(oaf);
|
||||||
}
|
}
|
||||||
|
|
||||||
void Compute::eval_empXX(const Rod2d & /*rod2D*/) {
|
void Compute::eval_empXX(const Rod2d &rod2D) {
|
||||||
//TODO
|
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;
|
||||||
|
agg.feed(empxx);
|
||||||
}
|
}
|
||||||
|
|
||||||
void Compute::eval_empYY(const Rod2d & /*rod2D*/) {
|
void Compute::eval_empYY(const Rod2d &rod2D) {
|
||||||
//TODO
|
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)/2.0 / time;
|
||||||
|
agg.feed(empxx);
|
||||||
}
|
}
|
||||||
|
|
||||||
void Compute::eval(const Rod2d &rod2D) {
|
void Compute::eval(const Rod2d &rod2D) {
|
||||||
@ -50,44 +60,53 @@ void Compute::eval(const Rod2d &rod2D) {
|
|||||||
}
|
}
|
||||||
start_rod = rod2D;
|
start_rod = rod2D;
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
Compute::Compute(Rod2d rod, Type t_type, size_t t_every, Simulation &sim)
|
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) {
|
: start_rod(std::move(rod)), every(t_every), time_step(0), type(t_type) {
|
||||||
auto time = sim.getMDeltaT() * static_cast<double>(every);
|
time = sim.getMDeltaT() * static_cast<double>(every);
|
||||||
switch (type) {
|
switch (type) {
|
||||||
case Type::msd:
|
case Type::msd:{
|
||||||
|
type_str = "msd";
|
||||||
target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
|
target = 4.0 * 0.5 * (rod.getDiff().trace()) * time;
|
||||||
break;
|
break;}
|
||||||
case Type::oaf:
|
case Type::oaf: {
|
||||||
|
type_str = "oaf";
|
||||||
target = std::exp(-rod.getDRot() * time);
|
target = std::exp(-rod.getDRot() * time);
|
||||||
break;
|
break;
|
||||||
case Type::empxx:
|
}
|
||||||
target = 0;
|
case Type::empxx: {
|
||||||
break;
|
type_str = "empxx";
|
||||||
case Type::empyy:
|
const double Dmean = 0.5 * (rod.getDiff().trace());
|
||||||
target = 0;
|
const double u = 4.0;
|
||||||
break;
|
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: {
|
||||||
|
type_str = "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 {
|
const LiveAgg &Compute::getAgg() const { return agg; }
|
||||||
return agg;
|
|
||||||
}
|
|
||||||
|
|
||||||
Compute::Type Compute::getType() const {
|
Compute::Type Compute::getType() const { return type; }
|
||||||
return type;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Compute::getTime() const {
|
double Compute::getTime() const { return static_cast<double>(every); }
|
||||||
return static_cast<double>(every);
|
|
||||||
}
|
|
||||||
|
|
||||||
double Compute::getTarget() const {
|
double Compute::getTarget() const { return target; }
|
||||||
return target;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Compute::getDifference() const {
|
double Compute::getDifference() const { return abs(agg.getMean() - target); }
|
||||||
return abs(agg.getMean() -target);
|
|
||||||
}
|
std::ostream &operator<<(std::ostream &os, const Compute &compute) {
|
||||||
|
os << "agg: " << compute.agg << " type: " << compute.type_str;
|
||||||
|
return os;
|
||||||
|
}
|
||||||
|
@ -5,17 +5,14 @@
|
|||||||
#ifndef MYPROJECT_COMPUTE_H
|
#ifndef MYPROJECT_COMPUTE_H
|
||||||
#define MYPROJECT_COMPUTE_H
|
#define MYPROJECT_COMPUTE_H
|
||||||
|
|
||||||
|
#include <ostream>
|
||||||
#include "LiveAgg.hpp"
|
#include "LiveAgg.hpp"
|
||||||
#include "Rod2d.hpp"
|
#include "Rod2d.hpp"
|
||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
class Compute {
|
class Compute {
|
||||||
|
public:
|
||||||
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);
|
||||||
|
|
||||||
@ -39,6 +36,7 @@ public:
|
|||||||
|
|
||||||
double getDifference() const;
|
double getDifference() const;
|
||||||
|
|
||||||
|
friend std::ostream &operator<<(std::ostream &os, const Compute &compute);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
Rod2d start_rod;
|
Rod2d start_rod;
|
||||||
@ -47,8 +45,8 @@ private:
|
|||||||
size_t time_step;
|
size_t time_step;
|
||||||
Type type;
|
Type type;
|
||||||
double target;
|
double target;
|
||||||
|
double time;
|
||||||
|
std::string type_str;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
#endif // MYPROJECT_COMPUTE_H
|
||||||
#endif //MYPROJECT_COMPUTE_H
|
|
||||||
|
@ -3,23 +3,13 @@
|
|||||||
//
|
//
|
||||||
|
|
||||||
#include "Integratoren2d_force.h"
|
#include "Integratoren2d_force.h"
|
||||||
|
#include "vec_trafos.h"
|
||||||
|
|
||||||
constexpr double kBT = 1.0;
|
|
||||||
using Vector = Eigen::Vector2d;
|
using Vector = Eigen::Vector2d;
|
||||||
|
using Matrix = Eigen::Matrix2d;
|
||||||
|
|
||||||
|
void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D,
|
||||||
Vector Integratoren2d_force::ortho(Vector e) {
|
Simulation & /*sim*/) {
|
||||||
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];
|
|
||||||
return mat;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) {
|
|
||||||
auto trans_lab = Vector({0, 0.01});
|
auto trans_lab = Vector({0, 0.01});
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
rod2D.setPos(rod2D.getPos() + trans_lab);
|
||||||
Eigen::Rotation2D<double> rot(0.1);
|
Eigen::Rotation2D<double> rot(0.1);
|
||||||
@ -28,144 +18,177 @@ void Integratoren2d_force::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/
|
|||||||
|
|
||||||
const Vector Integratoren2d_force::unitVec = {1, 0};
|
const Vector Integratoren2d_force::unitVec = {1, 0};
|
||||||
|
|
||||||
|
void Integratoren2d_force::Set1_Euler(
|
||||||
|
Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Vector(Vector, Vector)> &force,
|
||||||
|
const std::function<double(Vector, Vector)> &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;
|
||||||
|
Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
|
||||||
|
|
||||||
void Integratoren2d_force::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) {
|
// Force
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren2d_force::Set1_Euler(Rod2d &rod2D, Simulation &sim,
|
|
||||||
const std::function<Vector(Vector, Vector)> &force,
|
|
||||||
const std::function<double(Vector, Vector)> &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
|
|
||||||
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
|
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
|
||||||
Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab;
|
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
|
||||||
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
||||||
F_trans *= sim.getMDeltaT() / kBT;
|
F_trans *= sim.getMDeltaT();
|
||||||
|
|
||||||
|
// rotation
|
||||||
//rotation
|
double rot =
|
||||||
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden.
|
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
|
||||||
Vector e = rod2D.getE();
|
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() * sim.getMDeltaT();
|
||||||
Vector new_e = e + (rot + M_rot) * ortho(e);
|
Vector new_e = e + (rot + M_rot) * ortho(e);
|
||||||
new_e.normalize();
|
new_e.normalize();
|
||||||
//apply
|
// apply
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
|
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
|
||||||
rod2D.setE(new_e);
|
rod2D.setE(new_e);
|
||||||
}
|
}
|
||||||
|
|
||||||
Vector Integratoren2d_force::Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim,
|
Vector Integratoren2d_force::Heun_predictor_pos(
|
||||||
const std::function<Vector(Vector, Vector)> &force) {
|
const Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Vector(Vector, Vector)> &force) {
|
||||||
auto standard_dev = sim.getSTD();
|
auto standard_dev = sim.getSTD();
|
||||||
Vector rand_pred = {sim.getNorm(standard_dev), sim.getNorm(standard_dev)}; //gaussian noise
|
Vector rand_pred = {sim.getNorm(standard_dev),
|
||||||
Vector trans_pred_part = rod2D.getDiff_Sqrt() * rand_pred; //translation vector in Particle System
|
sim.getNorm(standard_dev)}; // gaussian noise
|
||||||
Vector trans_pred_lab = rod2D.getE_Base_matrix() * trans_pred_part;
|
Vector trans_pred_part =
|
||||||
|
rod2D.getDiff_Sqrt() *
|
||||||
|
rand_pred; // translation vector in Particle System
|
||||||
|
Vector trans_pred_lab = rotation_Matrix(rod2D.getE()) * trans_pred_part;
|
||||||
|
|
||||||
Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE());
|
Vector F_pred_lab = force(rod2D.getPos(), rod2D.getE());
|
||||||
Vector F_pred_part = rod2D.getE_Base_matrix().inverse() * F_pred_lab;
|
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 = rod2D.getDiff_Sqrt() * F_pred_part;
|
||||||
F_pred_trans *= sim.getMDeltaT() / kBT;
|
F_pred_trans *= sim.getMDeltaT() ;
|
||||||
|
|
||||||
return F_pred_trans + trans_pred_lab;
|
return F_pred_trans + trans_pred_lab;
|
||||||
}
|
}
|
||||||
|
|
||||||
Vector Integratoren2d_force::Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim,
|
Vector Integratoren2d_force::Heun_predictor_rot(
|
||||||
const std::function<double(Vector, Vector)> &torque) {
|
const Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<double(Vector, Vector)> &torque) {
|
||||||
auto std = sim.getSTD();
|
auto std = sim.getSTD();
|
||||||
double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden.
|
double rot_predict =
|
||||||
Vector e = rod2D.getE();
|
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
|
||||||
double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() / kBT * sim.getMDeltaT();
|
const Vector& e = rod2D.getE();
|
||||||
|
double M_predict_rot = torque(rod2D.getPos(), rod2D.getE()) *
|
||||||
|
rod2D.getDRot() * sim.getMDeltaT();
|
||||||
Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e);
|
Vector e_change_predict = (rot_predict + M_predict_rot) * ortho(e);
|
||||||
|
|
||||||
return e_change_predict;
|
return e_change_predict;
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void Integratoren2d_force::Set2_Heun(
|
||||||
void Integratoren2d_force::Set2_Heun(Rod2d &rod2D, Simulation &sim,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
const std::function<Vector(Vector, Vector)> &force,
|
const std::function<Vector(Vector, Vector)> &force,
|
||||||
const std::function<double(Vector, Vector)> &torque) {
|
const std::function<double(Vector, Vector)> &torque) {
|
||||||
|
|
||||||
Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force);
|
Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force);
|
||||||
Vector pos_predictor = rod2D.getPos() + delta_pos_predictor;
|
Vector pos_predictor = rod2D.getPos() + delta_pos_predictor;
|
||||||
|
|
||||||
Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque);
|
Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque);
|
||||||
Vector e_predict = rod2D.getE() + delta_e_predictor;
|
Vector e_predictor = rod2D.getE() + delta_e_predictor;
|
||||||
e_predict.normalize();
|
e_predictor.normalize();
|
||||||
|
|
||||||
|
|
||||||
Rod2d pred_rod = rod2D;
|
Rod2d pred_rod = rod2D;
|
||||||
pred_rod.setPos(pos_predictor);
|
pred_rod.setPos(pos_predictor);
|
||||||
pred_rod.setE(e_predict);
|
pred_rod.setE(e_predictor);
|
||||||
|
|
||||||
Vector delta_pos_future = Heun_predictor_pos(pred_rod, sim, force);
|
|
||||||
Vector delta_e_future = Heun_predictor_rot(pred_rod, sim, torque);
|
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 =
|
||||||
Vector e_integrated = rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future);
|
rod2D.getE() + 0.5 * (delta_e_predictor + delta_e_future);
|
||||||
//apply
|
// apply
|
||||||
rod2D.setPos(pos_integrated);
|
rod2D.setPos(pos_predictor);
|
||||||
rod2D.setE(e_integrated.normalized());
|
rod2D.setE(e_integrated.normalized());
|
||||||
}
|
}
|
||||||
|
|
||||||
void Integratoren2d_force::Set3_Exact(Rod2d &rod2D, Simulation &sim,
|
void Integratoren2d_force::Set3_Exact(
|
||||||
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque) {
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
auto std = sim.getSTD(); // sqrt(2*delta_T)
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque) {
|
||||||
//translation
|
auto std = sim.getSTD(); // sqrt(2*delta_T)
|
||||||
Vector rand = {sim.getNorm(std), sim.getNorm(std)}; //gaussian noise
|
// translation
|
||||||
Vector trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System
|
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_lab = force(rod2D.getPos(), rod2D.getE());
|
||||||
Vector F_part = rod2D.getE_Base_matrix().inverse() * F_lab;
|
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
|
||||||
|
|
||||||
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
||||||
F_trans *= sim.getMDeltaT() / kBT;
|
F_trans *= sim.getMDeltaT();
|
||||||
Vector trans_lab = rod2D.getE_Base_matrix() * trans_part;
|
Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
|
||||||
|
|
||||||
//rotation
|
// rotation
|
||||||
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden.
|
double rot =
|
||||||
|
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
|
||||||
Vector e = rod2D.getE();
|
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() * sim.getMDeltaT();
|
||||||
auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT();
|
auto correction =
|
||||||
|
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
|
||||||
Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e;
|
Vector new_e = e + (rot + M_rot) * ortho(e) + correction * e;
|
||||||
new_e.normalize();
|
new_e.normalize();
|
||||||
//apply
|
// apply
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
|
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
|
||||||
rod2D.setE(new_e);
|
rod2D.setE(new_e);
|
||||||
}
|
}
|
||||||
|
|
||||||
void Integratoren2d_force::Set4_BDAS(Rod2d &rod2D, Simulation &sim,
|
void Integratoren2d_force::Set4_BDAS(
|
||||||
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> /*torque*/) {
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
auto std = sim.getSTD();
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque) {
|
||||||
//translation
|
auto std = sim.getSTD(); // sqrt(2*delta_T)
|
||||||
auto rand = Vector(sim.getNorm(std), sim.getNorm(std));
|
Vector e = rod2D.getE();
|
||||||
auto trans_part = rod2D.getDiff_Sqrt() * rand; //translation vector in Particle System
|
// translation
|
||||||
Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE())));
|
Vector rand = {sim.getNorm(std), sim.getNorm(std)}; // gaussian noise
|
||||||
auto trans_lab = rotMat * trans_part;
|
Vector trans_part = rod2D.getDiff_Sqrt() * rand;
|
||||||
|
Vector trans_lab = rotation_Matrix(rod2D.getE()) * trans_part;
|
||||||
|
|
||||||
|
// Force
|
||||||
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
|
Vector F_lab = force(rod2D.getPos(), rod2D.getE());
|
||||||
Vector F_part = rotMat.inverse() * F_lab;
|
Vector F_part = rotation_Matrix(rod2D.getE()).inverse() * F_lab;
|
||||||
|
|
||||||
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
Vector F_trans = rod2D.getDiff_Sqrt() * F_part;
|
||||||
F_trans *= sim.getMDeltaT() / kBT;
|
F_trans *= sim.getMDeltaT();
|
||||||
|
|
||||||
|
// rotation
|
||||||
|
double rot =
|
||||||
|
sim.getNorm(std) * rod2D.getDRot_Sqrt(); // rotationsmatrix verwenden.
|
||||||
|
|
||||||
auto rot = Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
|
double M_rot = torque(rod2D.getPos(), rod2D.getE()) * rod2D.getDRot() * sim.getMDeltaT();
|
||||||
Eigen::Rotation2Dd rotation(rot);
|
Vector new_e = Eigen::Rotation2D<double>(rot + M_rot) * e;
|
||||||
auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized();
|
new_e.normalize();
|
||||||
|
// apply
|
||||||
// Normalisation should not be necessary if a proper angular representation is used.
|
rod2D.setPos(rod2D.getPos() + trans_lab + F_trans);
|
||||||
// But with vector e it is done in case of numerical errors
|
rod2D.setE(new_e);
|
||||||
|
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
|
||||||
rod2D.setE(e_new);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
void Integratoren2d_force::Set5_MBD(Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque) {
|
||||||
|
|
||||||
|
Vector delta_pos_predictor = Heun_predictor_pos(rod2D, sim, force);
|
||||||
|
[[maybe_unused]] Vector pos_predictor = rod2D.getPos() + delta_pos_predictor;
|
||||||
|
|
||||||
|
Vector delta_e_predictor = Heun_predictor_rot(rod2D, sim, torque);
|
||||||
|
Vector e_predictor = rod2D.getE() + delta_e_predictor;
|
||||||
|
e_predictor.normalize();
|
||||||
|
|
||||||
|
|
||||||
|
auto std = sim.getSTD();
|
||||||
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
|
|
||||||
|
|
||||||
|
Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predictor)*rod2D.getDiff());
|
||||||
|
Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL();
|
||||||
|
|
||||||
|
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
|
||||||
|
Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand;
|
||||||
|
|
||||||
|
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
|
Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predictor));
|
||||||
|
//apply
|
||||||
|
rod2D.setPos(pos_integrated);
|
||||||
|
rod2D.setE(e_integrated.normalized());
|
||||||
|
}
|
||||||
|
@ -8,41 +8,44 @@
|
|||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
class Integratoren2d_force {
|
class Integratoren2d_force {
|
||||||
public:
|
public:
|
||||||
static void Set1_Euler(Rod2d &rod2D, Simulation &sim,
|
static void Set1_Euler(
|
||||||
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
|
||||||
|
&force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
||||||
|
|
||||||
|
static void Set2_Heun(
|
||||||
|
Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
|
||||||
|
&force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
||||||
|
|
||||||
static void Set2_Heun(Rod2d &rod2D, Simulation &sim,
|
static void Set3_Exact(
|
||||||
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> &force,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque);
|
||||||
|
|
||||||
static void Set3_Exact(Rod2d &rod2D, Simulation &sim,
|
static void Set4_BDAS(
|
||||||
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
Rod2d &rod2D, Simulation &sim,
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque);
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque);
|
||||||
static void Set4_BDAS(Rod2d &rod2D, Simulation &sim,
|
|
||||||
std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)> force,
|
|
||||||
std::function<double(Eigen::Vector2d, Eigen::Vector2d)> torque);
|
|
||||||
|
|
||||||
static const Eigen::Vector2d unitVec;
|
static const Eigen::Vector2d unitVec;
|
||||||
|
|
||||||
static void Set5_MBD(Rod2d &rod2D, Simulation &sim);
|
static void Set5_MBD(Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>& force,
|
||||||
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)>& torque);
|
||||||
|
|
||||||
static void Set0_deterministic(Rod2d &rod2D, Simulation &sim);
|
static void Set0_deterministic(Rod2d &rod2D, Simulation &sim);
|
||||||
|
|
||||||
private:
|
private:
|
||||||
static Eigen::Vector2d ortho(Eigen::Vector2d e);
|
static Eigen::Vector2d Heun_predictor_pos(
|
||||||
|
const Rod2d &rod2D, Simulation &sim,
|
||||||
|
const std::function<Eigen::Vector2d(Eigen::Vector2d, Eigen::Vector2d)>
|
||||||
|
&force);
|
||||||
|
|
||||||
static Eigen::Matrix2d MatTraf(Eigen::Matrix<double, 2, 1> e);
|
static Eigen::Vector2d Heun_predictor_rot(
|
||||||
|
const Rod2d &rod2D, Simulation &sim,
|
||||||
static Eigen::Vector2d Heun_predictor_pos(const Rod2d &rod2D, Simulation &sim,
|
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
||||||
const std::function<Eigen::Vector2d(Eigen::Vector2d,
|
|
||||||
Eigen::Vector2d)> &force);
|
|
||||||
|
|
||||||
static Eigen::Vector2d
|
|
||||||
Heun_predictor_rot(const Rod2d &rod2D, Simulation &sim,
|
|
||||||
const std::function<double(Eigen::Vector2d, Eigen::Vector2d)> &torque);
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
@ -3,12 +3,10 @@
|
|||||||
//
|
//
|
||||||
|
|
||||||
#include "Integratoren2d_forceless.h"
|
#include "Integratoren2d_forceless.h"
|
||||||
|
#include "vec_trafos.h"
|
||||||
|
|
||||||
Eigen::Vector2d Integratoren2d_forceless::ortho(Eigen::Vector2d e) {
|
void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D,
|
||||||
return Eigen::Vector2d{-e[1], e[0]};
|
Simulation & /*sim*/) {
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*sim*/) {
|
|
||||||
auto trans_lab = Eigen::Vector2d({0, 0.01});
|
auto trans_lab = Eigen::Vector2d({0, 0.01});
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
rod2D.setPos(rod2D.getPos() + trans_lab);
|
||||||
Eigen::Rotation2D<double> rot(0.1);
|
Eigen::Rotation2D<double> rot(0.1);
|
||||||
@ -16,156 +14,110 @@ void Integratoren2d_forceless::Set0_deterministic(Rod2d &rod2D, Simulation & /*s
|
|||||||
}
|
}
|
||||||
|
|
||||||
void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) {
|
void Integratoren2d_forceless::Set1_Euler(Rod2d &rod2D, Simulation &sim) {
|
||||||
auto std = sim.getSTD(); // sqrt(2*delta_T)
|
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);
|
|
||||||
|
|
||||||
*/
|
|
||||||
|
|
||||||
//rotation
|
|
||||||
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();// rotationsmatrix verwenden.
|
|
||||||
Eigen::Vector2d e = rod2D.getE();
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
Eigen::Vector2d new_e = e + rot * ortho(e);
|
// translation
|
||||||
new_e.normalize();
|
Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)};
|
||||||
//apply
|
Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand;
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
|
||||||
rod2D.setE(new_e);
|
|
||||||
|
|
||||||
|
// rotation
|
||||||
|
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
|
Eigen::Vector2d new_e = e + rot * ortho(e);
|
||||||
|
|
||||||
|
// apply
|
||||||
|
rod2D.setPos(rod2D.getPos() + trans_lab);
|
||||||
|
rod2D.setE(new_e.normalized());
|
||||||
}
|
}
|
||||||
|
|
||||||
void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) {
|
void Integratoren2d_forceless::Set2_Heun(Rod2d &rod2D, Simulation &sim) {
|
||||||
// Predict with Euler
|
// Predict with Euler
|
||||||
auto std = sim.getSTD();
|
auto std = sim.getSTD();
|
||||||
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
|
|
||||||
Eigen::Vector2d rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
|
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 =
|
||||||
Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part;
|
rod2D.getDiff_Sqrt() * rand;
|
||||||
|
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
|
||||||
Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab;
|
Eigen::Vector2d pos_predictor = rod2D.getPos() + trans_lab;
|
||||||
|
|
||||||
/*
|
// rotation
|
||||||
|
double rot_predict = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
TODO
|
|
||||||
|
|
||||||
Siehe Set_1: Hier auch über e und orth(e)
|
|
||||||
|
|
||||||
*/
|
|
||||||
|
|
||||||
//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 delta_e_predict = rot_predict * ortho(e);
|
||||||
Eigen::Vector2d e_predict = (e + delta_e_predict).normalized();
|
Eigen::Vector2d e_predict = (e + delta_e_predict).normalized();
|
||||||
|
|
||||||
|
|
||||||
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
|
Eigen::Vector2d e_integrated =
|
||||||
//TODO pos integration with future system ???
|
e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
|
||||||
// Eigen::Vector2d pos_integrated = rod2D.getPos() + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
|
|
||||||
|
|
||||||
//apply
|
// apply
|
||||||
rod2D.setPos(pos_predictor);
|
rod2D.setPos(pos_predictor); // pos with Euler
|
||||||
rod2D.setE(e_integrated.normalized());
|
rod2D.setE(e_integrated.normalized());
|
||||||
}
|
}
|
||||||
|
|
||||||
constexpr double kBT = 1.0;
|
|
||||||
|
|
||||||
void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) {
|
void Integratoren2d_forceless::Set3_Exact(Rod2d &rod2D, Simulation &sim) {
|
||||||
auto std = sim.getSTD();
|
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
|
|
||||||
Eigen::Vector2d trans_lab = rod2D.getE_Base_matrix() * trans_part;
|
|
||||||
|
|
||||||
/*
|
|
||||||
|
|
||||||
TODO: Siehe Set1
|
|
||||||
|
|
||||||
*/
|
|
||||||
|
|
||||||
//rotation
|
|
||||||
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
|
||||||
Eigen::Vector2d e = rod2D.getE();
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
auto correction = -0.5 * pow(rod2D.getDRot(), 2) / pow(kBT, 2) * sim.getMDeltaT();
|
// translation
|
||||||
|
Eigen::Vector2d rand = {sim.getNorm(std), sim.getNorm(std)};
|
||||||
|
Eigen::Vector2d trans_part = rod2D.getDiff_Sqrt() * rand;
|
||||||
|
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
|
||||||
|
|
||||||
|
// rotation
|
||||||
|
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
|
auto correction =
|
||||||
|
-0.5 * pow(rod2D.getDRot(), 2) * sim.getMDeltaT();
|
||||||
Eigen::Vector2d delta_e = correction * e + rot * ortho(e);
|
Eigen::Vector2d delta_e = correction * e + rot * ortho(e);
|
||||||
|
|
||||||
/*
|
// apply
|
||||||
TODO Warum hier kBT definiert?
|
|
||||||
*/
|
|
||||||
|
|
||||||
//apply
|
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
rod2D.setPos(rod2D.getPos() + trans_lab);
|
||||||
rod2D.setE((e + delta_e).normalized());
|
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
|
// Normalisation should not be necessary if a proper angular representation
|
||||||
|
// is used. But with vector e it is done in case of numerical errors
|
||||||
void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) {
|
void Integratoren2d_forceless::Set4_BDAS(Rod2d &rod2D, Simulation &sim) {
|
||||||
auto std = sim.getSTD();
|
auto std = sim.getSTD();
|
||||||
//translation
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
|
// translation
|
||||||
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
|
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 =
|
||||||
Eigen::Rotation2D rotMat(acos(unitVec.dot(rod2D.getE())));
|
rod2D.getDiff_Sqrt() * rand; // translation vector in Particle System
|
||||||
auto trans_lab = rotMat * trans_part;
|
Eigen::Vector2d trans_lab = rotation_Matrix(e) * trans_part;
|
||||||
|
|
||||||
/*
|
auto rot =
|
||||||
TODO: Warum hier nicht die Matrix von oben?
|
Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
|
||||||
|
|
||||||
*/
|
|
||||||
|
|
||||||
auto rot = Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
|
|
||||||
Eigen::Rotation2Dd rotation(rot);
|
Eigen::Rotation2Dd rotation(rot);
|
||||||
auto e_new = (rotation.toRotationMatrix() * rod2D.getE()).normalized();
|
auto e_new = (rotation.toRotationMatrix() * e).normalized();
|
||||||
|
|
||||||
/*
|
|
||||||
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
|
|
||||||
|
|
||||||
rod2D.setPos(rod2D.getPos() + trans_lab);
|
rod2D.setPos(rod2D.getPos() + trans_lab);
|
||||||
rod2D.setE(e_new);
|
rod2D.setE(e_new);
|
||||||
}
|
}
|
||||||
|
|
||||||
Eigen::Matrix2d Integratoren2d_forceless::e_2_matrix(Eigen::Vector2d m_e) {
|
void Integratoren2d_forceless::Set5_MBD(Rod2d & rod2D,
|
||||||
Eigen::Matrix2d mat;
|
Simulation & sim) {
|
||||||
mat << m_e[0], -m_e[1],
|
|
||||||
m_e[1], m_e[0];
|
|
||||||
return mat;
|
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren2d_forceless::Set5_MBD(Rod2d &/*rod2D*/, Simulation &/*sim*/) {
|
|
||||||
/*
|
|
||||||
auto std = sim.getSTD();
|
auto std = sim.getSTD();
|
||||||
|
Eigen::Vector2d e = rod2D.getE();
|
||||||
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;
|
|
||||||
|
|
||||||
//rotation
|
//rotation
|
||||||
auto rot_predict = Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());// rotationsmatrix verwenden.
|
double rot_predict = sim.getNorm(std) *
|
||||||
auto e = rod2D.getE();
|
rod2D.getDRot_Sqrt();
|
||||||
auto delta_e_predict = rot_predict * ortho(e);
|
auto delta_e_predict = rot_predict * ortho(e);
|
||||||
auto e_predict = (e + delta_e_predict).normalized();
|
auto e_predict = (e + delta_e_predict).normalized();
|
||||||
|
|
||||||
auto std_combined = rod2D.getDiff() + rod2D.getE_Base_matrix() * rod2D.getDiff(); //old + new diff
|
|
||||||
|
|
||||||
auto rot = Eigen::Rotation2D<double>(sim.getNorm(std) * rod2D.getDRot_Sqrt());
|
Eigen::Matrix2d Diff_combined = 0.5*(rod2D.getDiff() + rotation_Matrix(e).inverse() * rotation_Matrix(e_predict)*rod2D.getDiff());
|
||||||
auto e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
|
Eigen::Matrix2d D_combined_sqrt = Diff_combined.llt().matrixL();
|
||||||
|
|
||||||
|
auto rand = Eigen::Vector2d(sim.getNorm(std), sim.getNorm(std));
|
||||||
|
Eigen::Vector2d pos_integrated =rod2D.getPos() + rotation_Matrix(e) * D_combined_sqrt * rand;
|
||||||
|
|
||||||
|
double rot = sim.getNorm(std) * rod2D.getDRot_Sqrt();
|
||||||
|
Eigen::Vector2d e_integrated = e + 0.5 * (rot * ortho(e) + rot * ortho(e_predict));
|
||||||
//apply
|
//apply
|
||||||
rod2D.setPos(pos_predictor);
|
rod2D.setPos(pos_integrated);
|
||||||
rod2D.setE(e_integrated.normalized());
|
rod2D.setE(e_integrated.normalized());
|
||||||
*/
|
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
@ -2,14 +2,13 @@
|
|||||||
// Created by jholder on 21.10.21.
|
// Created by jholder on 21.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
#ifndef MYPROJECT_INTEGRATOREN2D_FORCELESS_H
|
#pragma once
|
||||||
#define MYPROJECT_INTEGRATOREN2D_FORCELESS_H
|
|
||||||
|
|
||||||
#include "Rod2d.hpp"
|
#include "Rod2d.hpp"
|
||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
class Integratoren2d_forceless {
|
class Integratoren2d_forceless {
|
||||||
public:
|
public:
|
||||||
static void Set1_Euler(Rod2d &rod2D, Simulation &sim);
|
static void Set1_Euler(Rod2d &rod2D, Simulation &sim);
|
||||||
|
|
||||||
static void Set2_Heun(Rod2d &rod2D, Simulation &sim);
|
static void Set2_Heun(Rod2d &rod2D, Simulation &sim);
|
||||||
@ -20,15 +19,7 @@ public:
|
|||||||
|
|
||||||
static const Eigen::Vector2d unitVec;
|
static const Eigen::Vector2d unitVec;
|
||||||
|
|
||||||
static Eigen::Matrix2d e_2_matrix(Eigen::Vector2d m_e);
|
|
||||||
|
|
||||||
static void Set5_MBD(Rod2d &rod2D, Simulation &sim);
|
static void Set5_MBD(Rod2d &rod2D, Simulation &sim);
|
||||||
|
|
||||||
static void Set0_deterministic(Rod2d &rod2D, Simulation &sim);
|
static void Set0_deterministic(Rod2d &rod2D, Simulation &sim);
|
||||||
|
|
||||||
private:
|
|
||||||
static Eigen::Vector2d ortho(Eigen::Vector2d e);
|
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
#endif //MYPROJECT_INTEGRATOREN2D_FORCELESS_H
|
|
||||||
|
@ -1,33 +1,26 @@
|
|||||||
#include "LiveAgg.hpp"
|
#include "LiveAgg.hpp"
|
||||||
|
|
||||||
#include <cmath>
|
#include <cmath>
|
||||||
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
|
double LiveAgg::getMean() const noexcept { return mean; }
|
||||||
{
|
|
||||||
return S / (num_of_data - 1);
|
|
||||||
}
|
|
||||||
|
|
||||||
double LiveAgg::getSEM() const noexcept
|
void LiveAgg::feed(double value) noexcept {
|
||||||
{
|
num_of_data += 1;
|
||||||
return S / (num_of_data - 1) / std::sqrt(num_of_data);
|
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; }
|
||||||
|
|
||||||
double LiveAgg::getMean() const noexcept
|
std::ostream &operator<<(std::ostream &os, const LiveAgg &agg) {
|
||||||
{
|
os << "num_of_data: " << agg.num_of_data << " mean: " << agg.mean << " S: " << agg.S;
|
||||||
return mean;
|
return os;
|
||||||
}
|
|
||||||
|
|
||||||
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;
|
|
||||||
}
|
}
|
||||||
|
@ -1,20 +1,20 @@
|
|||||||
#ifndef LIVEAGG_H
|
#pragma once
|
||||||
#define LIVEAGG_H
|
|
||||||
|
|
||||||
|
#include <ostream>
|
||||||
|
|
||||||
class LiveAgg
|
class LiveAgg {
|
||||||
{
|
public:
|
||||||
public:
|
LiveAgg();
|
||||||
LiveAgg();
|
[[nodiscard]] double getSD() const noexcept;
|
||||||
[[nodiscard]] double getSD() const noexcept;
|
[[nodiscard]] double getSEM() const noexcept;
|
||||||
[[nodiscard]] double getSEM() const noexcept;
|
[[nodiscard]] double getMean() const noexcept;
|
||||||
[[nodiscard]] double getMean() const noexcept;
|
[[nodiscard]] int getNumPoints() const noexcept;
|
||||||
[[nodiscard]] int getNumPoints() const noexcept;
|
void feed(double value) noexcept;
|
||||||
void feed(double value) noexcept;
|
|
||||||
private:
|
friend std::ostream &operator<<(std::ostream &os, const LiveAgg &agg);
|
||||||
int num_of_data;
|
|
||||||
double mean;
|
private:
|
||||||
double S;
|
int num_of_data;
|
||||||
|
double mean;
|
||||||
|
double S;
|
||||||
};
|
};
|
||||||
|
|
||||||
#endif // LIVEAGG_H
|
|
||||||
|
@ -1,28 +1,30 @@
|
|||||||
#include "Rod2d.hpp"
|
#include "Rod2d.hpp"
|
||||||
|
|
||||||
#include "Eigen/Dense"
|
#include "Eigen/Dense"
|
||||||
|
|
||||||
constexpr double M_Pl = 3.141592653589793238462643383279502884; /* pi */
|
constexpr double M_Pl = 3.141592653589793238462643383279502884; /* pi */
|
||||||
|
|
||||||
void sqrt(Eigen::Matrix2d &mat){
|
void sqrt(Eigen::Matrix2d &mat) {
|
||||||
const auto size = static_cast<size_t>(mat.size());
|
const auto size = static_cast<size_t>(mat.size());
|
||||||
for (size_t i = 0; i < size; i++) {
|
for (size_t i = 0; i < size; i++) {
|
||||||
mat.data()[i] = sqrt(mat.data()[i]);
|
mat.data()[i] = sqrt(mat.data()[i]);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) {
|
Rod2d::Rod2d(double L) : m_pos({0, 0}), m_e({1, 0}) {
|
||||||
assert(L>=1.0);
|
assert(L >= 1.0);
|
||||||
if (L == 1.0) {
|
if (L == 1.0) {
|
||||||
m_D_rot = 3.0;
|
m_D_rot = 3.0;
|
||||||
m_Diff << 1, 0,
|
m_Diff << 1, 0, 0, 1;
|
||||||
0, 1;
|
|
||||||
} else {
|
} else {
|
||||||
const double D0 = 3.0 * M_Pl / L;
|
const double D0 = 3.0 * M_Pl / L;
|
||||||
const double p = 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_para =
|
||||||
auto D_ortho = D0 / (M_Pl * 2.0) * (log(p) + 0.839 + 0.185 / p + 0.233 / (p * p));
|
D0 / (M_Pl * 2.0) * (log(p) - 0.207 + 0.980 / p - 0.133 / (p * p));
|
||||||
m_D_rot = 3 * D0 / (M_Pl * L * L) * (log(p) - 0.662 + 0.917 / p - 0.050 / (p * p));
|
auto D_ortho =
|
||||||
m_Diff << D_para, 0,
|
D0 / (M_Pl * 2.0) * (log(p) + 0.839 + 0.185 / p + 0.233 / (p * p));
|
||||||
0, D_ortho;
|
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_Diff_sqrt = m_Diff;
|
||||||
m_D_rot_sqrt = sqrt(m_D_rot);
|
m_D_rot_sqrt = sqrt(m_D_rot);
|
||||||
@ -34,43 +36,19 @@ void Rod2d::reset() {
|
|||||||
m_e = {1, 0};
|
m_e = {1, 0};
|
||||||
}
|
}
|
||||||
|
|
||||||
void Rod2d::setPos(const Eigen::Vector2d &Pos) {
|
void Rod2d::setPos(const Eigen::Vector2d &Pos) { m_pos = Pos; }
|
||||||
m_pos = Pos;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Rod2d::getDRot() const {
|
double Rod2d::getDRot() const { return m_D_rot; }
|
||||||
return m_D_rot;
|
|
||||||
}
|
|
||||||
|
|
||||||
const Eigen::Vector2d &Rod2d::getPos() const {
|
const Eigen::Vector2d &Rod2d::getPos() const { return m_pos; }
|
||||||
return m_pos;
|
|
||||||
}
|
|
||||||
|
|
||||||
const Eigen::Matrix2d &Rod2d::getDiff() const {
|
const Eigen::Matrix2d &Rod2d::getDiff() const { return m_Diff; }
|
||||||
return m_Diff;
|
|
||||||
}
|
|
||||||
|
|
||||||
const Eigen::Matrix2d &Rod2d::getDiff_Sqrt() const {
|
const Eigen::Matrix2d &Rod2d::getDiff_Sqrt() const { return m_Diff_sqrt; }
|
||||||
return m_Diff_sqrt;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Rod2d::getDRot_Sqrt() const {
|
double Rod2d::getDRot_Sqrt() const { return m_D_rot_sqrt; }
|
||||||
return m_D_rot_sqrt;
|
|
||||||
}
|
|
||||||
|
|
||||||
const Eigen::Vector2d &Rod2d::getE() const {
|
const Eigen::Vector2d &Rod2d::getE() const { return m_e; }
|
||||||
return m_e;
|
|
||||||
}
|
|
||||||
|
|
||||||
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];
|
|
||||||
return mat;
|
|
||||||
}
|
|
||||||
|
|
||||||
|
void Rod2d::setE(const Eigen::Vector2d &mE) { m_e = mE; }
|
||||||
|
|
||||||
|
@ -1,9 +1,8 @@
|
|||||||
#pragma once
|
#pragma once
|
||||||
#include <Eigen/Geometry>
|
|
||||||
#include <Eigen/Dense>
|
#include <Eigen/Dense>
|
||||||
class Rod2d
|
#include <Eigen/Geometry>
|
||||||
{
|
class Rod2d {
|
||||||
public:
|
public:
|
||||||
explicit Rod2d(double L);
|
explicit Rod2d(double L);
|
||||||
void reset();
|
void reset();
|
||||||
|
|
||||||
@ -23,15 +22,11 @@ public:
|
|||||||
|
|
||||||
[[nodiscard]] double getDRot_Sqrt() const;
|
[[nodiscard]] double getDRot_Sqrt() const;
|
||||||
|
|
||||||
Eigen::Matrix2d getE_Base_matrix() const;
|
private:
|
||||||
|
|
||||||
private:
|
|
||||||
Eigen::Matrix2d m_Diff;
|
Eigen::Matrix2d m_Diff;
|
||||||
Eigen::Matrix2d m_Diff_sqrt;
|
Eigen::Matrix2d m_Diff_sqrt;
|
||||||
double m_D_rot;
|
double m_D_rot;
|
||||||
double m_D_rot_sqrt;
|
double m_D_rot_sqrt;
|
||||||
Eigen::Vector2d m_pos; // position
|
Eigen::Vector2d m_pos; // position
|
||||||
Eigen::Vector2d m_e;
|
Eigen::Vector2d m_e;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
|
@ -5,17 +5,15 @@
|
|||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
double Simulation::getNorm(double t_norm) {
|
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 {
|
double Simulation::getSTD() const { return m_std; }
|
||||||
return m_delta_T;
|
|
||||||
}
|
|
||||||
|
|
||||||
double Simulation::getSTD() const {
|
|
||||||
return m_std;
|
|
||||||
}
|
|
||||||
|
@ -2,25 +2,22 @@
|
|||||||
// Created by jholder on 21.10.21.
|
// Created by jholder on 21.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
#ifndef MYPROJECT_SIMULATION_H
|
#pragma once
|
||||||
#define MYPROJECT_SIMULATION_H
|
|
||||||
|
|
||||||
|
|
||||||
#include <random>
|
#include <random>
|
||||||
|
|
||||||
class Simulation {
|
class Simulation {
|
||||||
public:
|
public:
|
||||||
explicit Simulation(double t_delta_T, size_t seed);
|
explicit Simulation(double t_delta_T, size_t seed);
|
||||||
double getNorm(double t_norm);
|
double getNorm(double t_norm);
|
||||||
private:
|
|
||||||
|
private:
|
||||||
double m_delta_T;
|
double m_delta_T;
|
||||||
double m_std;
|
double m_std;
|
||||||
std::mt19937_64 m_generator;
|
std::mt19937_64 m_generator;
|
||||||
std::normal_distribution<double> m_norm;
|
std::normal_distribution<double> m_norm;
|
||||||
public:
|
|
||||||
|
public:
|
||||||
[[nodiscard]] double getMDeltaT() const;
|
[[nodiscard]] double getMDeltaT() const;
|
||||||
[[nodiscard]] double getSTD() const;
|
[[nodiscard]] double getSTD() const;
|
||||||
};
|
};
|
||||||
|
|
||||||
|
|
||||||
#endif //MYPROJECT_SIMULATION_H
|
|
||||||
|
@ -1,130 +0,0 @@
|
|||||||
#include "Integratoren.hpp"
|
|
||||||
constexpr double kBT = 1.0;
|
|
||||||
|
|
||||||
void Integratoren::integrateITOEuler_1(Rod2d &rod)
|
|
||||||
{
|
|
||||||
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
|
|
||||||
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
|
|
||||||
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
|
|
||||||
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
|
|
||||||
|
|
||||||
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double tmpe1 = rod.e1 - rodOrtho * rod.e2;
|
|
||||||
double tmpe2 = rod.e2 + rodOrtho * rod.e1;
|
|
||||||
|
|
||||||
rod.x1 = tmpx1;
|
|
||||||
rod.x2 = tmpx2;
|
|
||||||
rod.e1 = tmpe1;
|
|
||||||
rod.e2 = tmpe2;
|
|
||||||
rod.normalize();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren::integrateITOHeun_2(Rod2d &rod)
|
|
||||||
{
|
|
||||||
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
|
|
||||||
|
|
||||||
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
|
|
||||||
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
|
|
||||||
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
|
|
||||||
|
|
||||||
|
|
||||||
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double e1 = rod.e1 - rodOrtho * rod.e1;
|
|
||||||
double e2 = rod.e2 + rodOrtho * rod.e2;
|
|
||||||
double norm = std::sqrt(e1 * e1 + e2 * e2);
|
|
||||||
e1 /= norm;
|
|
||||||
e2 /= norm;
|
|
||||||
|
|
||||||
double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT);
|
|
||||||
double rodOrthoRand = rod.dis(rod.generator);
|
|
||||||
double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrthoRand;//TODO Check if 0.5 faktor
|
|
||||||
double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrthoRand;
|
|
||||||
|
|
||||||
rod.x1 = tmpx1;
|
|
||||||
rod.x2 = tmpx2;
|
|
||||||
rod.e1 = tmpe1;
|
|
||||||
rod.e2 = tmpe2;
|
|
||||||
rod.normalize();
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
void Integratoren::integrateExact_3(Rod2d &rod)
|
|
||||||
{
|
|
||||||
//DLOG_SCOPE_F(1,"Integrateing Rod2d:");
|
|
||||||
|
|
||||||
double movePara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double moveOrtho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
|
|
||||||
double tmpx1 = rod.x1 + movePara * rod.e1 - moveOrtho * rod.e2;
|
|
||||||
double tmpx2 = rod.x2 + movePara * rod.e2 + moveOrtho * rod.e1;
|
|
||||||
|
|
||||||
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double korr = -0.5 * rod.Drot * rod.Drot / kBT / kBT * rod.deltaT;
|
|
||||||
double tmpe1 = rod.e1 - rodOrtho * rod.e2 + korr * rod.e1;
|
|
||||||
double tmpe2 = rod.e2 + rodOrtho * rod.e1 + korr * rod.e2;
|
|
||||||
|
|
||||||
rod.x1 = tmpx1;
|
|
||||||
rod.x2 = tmpx2;
|
|
||||||
rod.e1 = tmpe1;
|
|
||||||
rod.e2 = tmpe2;
|
|
||||||
rod.normalize();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren::integrateBDAS_4(Rod2d &rod)
|
|
||||||
{
|
|
||||||
double deltaRTpara = std::sqrt(2.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double deltaRTortho = std::sqrt(2.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double tmpx1 = rod.x1 + deltaRTpara * rod.e1 - deltaRTortho * rod.e2;
|
|
||||||
double tmpx2 = rod.x2 + deltaRTpara * rod.e2 + deltaRTortho * rod.e1;
|
|
||||||
|
|
||||||
double deltaPhi = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double tmpe1 = rod.e1 * std::cos(deltaPhi) - rod.e2 * std::sin(deltaPhi);
|
|
||||||
double tmpe2 = rod.e1 * std::sin(deltaPhi) + rod.e2 * std::cos(deltaPhi);
|
|
||||||
|
|
||||||
rod.x1 = tmpx1;
|
|
||||||
rod.x2 = tmpx2;
|
|
||||||
rod.e1 = tmpe1;
|
|
||||||
rod.e2 = tmpe2;
|
|
||||||
rod.normalize();
|
|
||||||
}
|
|
||||||
|
|
||||||
void Integratoren::integrateMBD_5(Rod2d &rod)
|
|
||||||
{
|
|
||||||
/*
|
|
||||||
double movePara = std::sqrt(2.0*rod.Dpara*rod.deltaT)*rod.dis(rod.generator); // Taken from euler
|
|
||||||
double moveOrtho = std::sqrt(2.0*rod.Dortho*rod.deltaT)*rod.dis(rod.generator); // taken from euler
|
|
||||||
double x1 = rod.x1 + movePara*rod.e1- moveOrtho*rod.e2;
|
|
||||||
double x2 = rod.x2 + movePara*rod.e2+ moveOrtho*rod.e1;
|
|
||||||
*/
|
|
||||||
// unused without force
|
|
||||||
double rodOrtho = std::sqrt(2.0 * rod.Drot * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double e1 = rod.e1 - rodOrtho * rod.e2;
|
|
||||||
double e2 = rod.e2 + rodOrtho * rod.e1;
|
|
||||||
double norm = std::sqrt(e1 * e1 + e2 * e2);
|
|
||||||
e1 /= norm;
|
|
||||||
e2 /= norm;
|
|
||||||
|
|
||||||
double movePara1 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D
|
|
||||||
double moveOrtho1 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double movePara2 = std::sqrt(1.0 * rod.Dpara * rod.deltaT) * rod.dis(rod.generator);// If difussion would be time dependend, average of the predicted D and now D
|
|
||||||
double moveOrtho2 = std::sqrt(1.0 * rod.Dortho * rod.deltaT) * rod.dis(rod.generator);
|
|
||||||
double tmpx1 = rod.x1 + (movePara1 * rod.e1 - moveOrtho1 * rod.e2) + (movePara2 * e1 - moveOrtho2 * e2);
|
|
||||||
double tmpx2 = rod.x2 + (movePara1 * rod.e2 + moveOrtho1 * rod.e1) + (movePara2 * e2 + moveOrtho2 * e1);
|
|
||||||
|
|
||||||
|
|
||||||
double rodOrtho1 = std::sqrt(2.0 * rod.Drot * rod.deltaT);
|
|
||||||
double rodOrtho2 = rod.dis(rod.generator);
|
|
||||||
|
|
||||||
double tmpe1 = rod.e1 - 0.5 * (rodOrtho1 * rod.e2 + rodOrtho1 * e2) * rodOrtho2;
|
|
||||||
double tmpe2 = rod.e2 + 0.5 * (rodOrtho1 * rod.e1 + rodOrtho1 * e1) * rodOrtho2;
|
|
||||||
|
|
||||||
rod.x1 = tmpx1;
|
|
||||||
rod.x2 = tmpx2;
|
|
||||||
rod.e1 = tmpe1;
|
|
||||||
rod.e2 = tmpe2;
|
|
||||||
rod.normalize();
|
|
||||||
}
|
|
@ -1,19 +0,0 @@
|
|||||||
#ifndef INTEGRATOREN_HPP
|
|
||||||
#define INTEGRATOREN_HPP
|
|
||||||
#include "Rod2d.hpp"
|
|
||||||
|
|
||||||
|
|
||||||
class Integratoren
|
|
||||||
{
|
|
||||||
public:
|
|
||||||
static void integrateITOEuler_1(Rod2d &rod);
|
|
||||||
|
|
||||||
static void integrateITOHeun_2(Rod2d &rod);
|
|
||||||
|
|
||||||
static void integrateExact_3(Rod2d &rod);
|
|
||||||
|
|
||||||
static void integrateBDAS_4(Rod2d &rod);
|
|
||||||
|
|
||||||
static void integrateMBD_5(Rod2d &rod);
|
|
||||||
};
|
|
||||||
#endif// INTEGRATOREN_HPP
|
|
@ -1,6 +0,0 @@
|
|||||||
#include "Simulation2d.hpp"
|
|
||||||
#include <vector>
|
|
||||||
#include<iostream>
|
|
||||||
#include "LiveAgg.hpp"
|
|
||||||
template <void (*T)(Rod2d &)>
|
|
||||||
std::vector<double> simulate(double deltaT, int seed, int meanSteps);
|
|
@ -1,157 +0,0 @@
|
|||||||
#ifndef SIMULATION2D_HPP
|
|
||||||
#define SIMULATION2D_HPP
|
|
||||||
#include "Rod2d.hpp"
|
|
||||||
#include "LiveAgg.hpp"
|
|
||||||
#include <fstream>
|
|
||||||
#include <iostream>
|
|
||||||
/* This is used to calculate the convergence order of the particles
|
|
||||||
*/
|
|
||||||
|
|
||||||
inline void printProgress(double progress)
|
|
||||||
{
|
|
||||||
const int barWidth = 70;
|
|
||||||
std::cout << "[";
|
|
||||||
const int pos = static_cast<int>(barWidth * progress);
|
|
||||||
for (int i = 0; i < barWidth; ++i) {
|
|
||||||
if (i < pos) {
|
|
||||||
std::cout << "=";
|
|
||||||
} else if (i == pos) {
|
|
||||||
std::cout << ">";
|
|
||||||
} else {
|
|
||||||
std::cout << " ";
|
|
||||||
}
|
|
||||||
}
|
|
||||||
std::cout << "] " << int(progress * 100.0) << " % \r";
|
|
||||||
std::cout.flush();
|
|
||||||
}
|
|
||||||
|
|
||||||
template<void (*T)(Rod2d &)>
|
|
||||||
void simulate(double L, int seed, const std::string &filename)
|
|
||||||
{
|
|
||||||
// we now simulate the msd after 10 steps
|
|
||||||
|
|
||||||
|
|
||||||
Rod2d rod(L, 1e-5, seed);
|
|
||||||
const int sampleWidth = 1000;
|
|
||||||
std::vector<LiveAgg> msd(sampleWidth, LiveAgg());
|
|
||||||
std::vector<LiveAgg> oaf(sampleWidth, LiveAgg());
|
|
||||||
std::vector<LiveAgg> empxx(sampleWidth, LiveAgg());
|
|
||||||
std::vector<LiveAgg> empyy(sampleWidth, LiveAgg());
|
|
||||||
std::vector<double> xPos(sampleWidth, 0.0);
|
|
||||||
std::vector<double> yPos(sampleWidth, 0.0);
|
|
||||||
std::vector<double> e1(sampleWidth, 0.0);
|
|
||||||
std::vector<double> e2(sampleWidth, 0.0);
|
|
||||||
std::vector<double> targetMSD(sampleWidth, 0.0);
|
|
||||||
std::vector<double> targetOAF(sampleWidth, 0.0);
|
|
||||||
std::vector<double> targetEMPXX(sampleWidth, 0.0);
|
|
||||||
std::vector<double> targetEMPYY(sampleWidth, 0.0);
|
|
||||||
std::vector<size_t> samplePoints(sampleWidth, 1);
|
|
||||||
size_t temp = 1;
|
|
||||||
const double Dmean = 0.5 * (rod.Dpara + rod.Dortho);
|
|
||||||
const double deltaD = rod.Dortho - rod.Dpara;
|
|
||||||
for (size_t i = 0; i < samplePoints.size(); ++i) {
|
|
||||||
samplePoints[i] = temp;
|
|
||||||
temp += 100;
|
|
||||||
const double time = static_cast<double>(temp) * rod.deltaT;
|
|
||||||
targetMSD[i] = 4.0 * 0.5 * (rod.Dpara + rod.Dortho) * time;
|
|
||||||
targetOAF[i] = std::exp(-rod.Drot * time);
|
|
||||||
const double u = 4.0;
|
|
||||||
targetEMPXX[i] = Dmean - deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
|
|
||||||
targetEMPYY[i] = Dmean + deltaD / 2.0 * (1 - exp(-u * rod.Drot * time)) / u / rod.Drot / time;
|
|
||||||
}
|
|
||||||
unsigned long long int iteration = 0;
|
|
||||||
//const unsigned long long int maxItr = ((long int)10000)*1000000;
|
|
||||||
const unsigned long long int maxItr = 10000000;
|
|
||||||
const unsigned long long int updateInter = 10000;
|
|
||||||
while (iteration < maxItr) {
|
|
||||||
T(rod);
|
|
||||||
iteration++;
|
|
||||||
for (size_t i = 0; i < samplePoints.size(); ++i) {
|
|
||||||
const size_t sample = samplePoints[i];
|
|
||||||
if (iteration % sample == 0) {
|
|
||||||
msd[i].feed((rod.x1 - xPos[i]) * (rod.x1 - xPos[i]) + (rod.x2 - yPos[i]) * (rod.x2 - yPos[i]));
|
|
||||||
oaf[i].feed(rod.e1 * e1[i] + rod.e2 * e2[i]);
|
|
||||||
empxx[i].feed(std::pow((rod.x1 - xPos[i]) * e1[i] + (rod.x2 - yPos[i]) * e2[i], 2));
|
|
||||||
empyy[i].feed(std::pow((rod.x1 - xPos[i]) * e2[i] - (rod.x2 - yPos[i]) * e1[i], 2));
|
|
||||||
xPos[i] = rod.x1;
|
|
||||||
yPos[i] = rod.x2;
|
|
||||||
e1[i] = rod.e1;
|
|
||||||
e2[i] = rod.e2;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
if (iteration % updateInter == 0) {
|
|
||||||
printProgress(static_cast<double>(iteration)/ maxItr);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
std::cout << "\n";
|
|
||||||
|
|
||||||
std::ofstream fileMSD(filename + "_" + std::to_string(L) + "_MSD.out");
|
|
||||||
std::ofstream fileOAF(filename + "_" + std::to_string(L) + "_OAF.out");
|
|
||||||
std::ofstream fileXX(filename + "_" + std::to_string(L) + "_XX.out");
|
|
||||||
std::ofstream fileYY(filename + "_" + std::to_string(L) + "_YY.out");
|
|
||||||
|
|
||||||
for (size_t i = 0; i < samplePoints.size(); ++i) {
|
|
||||||
fileMSD << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetMSD[i] << "\t" << msd[i].getMean() << "\t" << msd[i].getSEM() << "\t" << msd[i].getSD() << "\t" << msd[i].getNumPoints() << "\n";
|
|
||||||
fileOAF << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetOAF[i] << "\t" << oaf[i].getMean() << "\t" << oaf[i].getSEM() << "\t" << oaf[i].getSD() << "\t" << oaf[i].getNumPoints() << "\n";
|
|
||||||
fileXX << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPXX[i] << "\t" << empxx[i].getMean() << "\t" << empxx[i].getSEM() << "\t" << empxx[i].getSD() << "\t" << empxx[i].getNumPoints() << "\n";
|
|
||||||
fileYY << static_cast<double>(samplePoints[i]) * rod.deltaT << "\t" << targetEMPYY[i] << "\t" << empyy[i].getMean() << "\t" << empyy[i].getSEM() << "\t" << empyy[i].getSD() << "\t" << empyy[i].getNumPoints() << "\n";
|
|
||||||
}
|
|
||||||
fileMSD.close();
|
|
||||||
fileOAF.close();
|
|
||||||
fileXX.close();
|
|
||||||
fileYY.close();
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
template<void (*T)(Rod2d &)>
|
|
||||||
void konvergenz(double L, int seed, double deltaT, const std::string &filename)
|
|
||||||
{
|
|
||||||
// we now simulate the msd after 10 steps
|
|
||||||
Rod2d rod(L, deltaT, seed);
|
|
||||||
const double time = 2;
|
|
||||||
const int steps = static_cast<int>(time / deltaT);
|
|
||||||
const int minSteps = 1000;
|
|
||||||
|
|
||||||
//const double targetMSD = 4.0*0.5*(rod.Dpara+rod.Dortho)*time;
|
|
||||||
const double targetOAF = std::exp(-rod.Drot * time);
|
|
||||||
/*const double u = 4.0;
|
|
||||||
const double Dmean = 0.5*(rod.Dpara+rod.Dortho);
|
|
||||||
const double deltaD = rod.Dortho - rod.Dpara;
|
|
||||||
const double targetEMPXX = Dmean-deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;
|
|
||||||
const double targetEMPYY = Dmean+deltaD/2.0*(1 - exp(-u*rod.Drot*time))/u/rod.Drot / time;*/
|
|
||||||
//LiveAgg msd;
|
|
||||||
LiveAgg oaf;
|
|
||||||
//LiveAgg empxx;
|
|
||||||
//LiveAgg empyy;
|
|
||||||
int meanCount = 0;
|
|
||||||
while (meanCount < minSteps or std::abs(targetOAF - oaf.getMean()) < oaf.getSEM() * 10) {
|
|
||||||
//reset the position
|
|
||||||
rod.reset();
|
|
||||||
for (int i = 0; i < steps; ++i) {
|
|
||||||
T(rod);
|
|
||||||
}
|
|
||||||
//msd.feed(rod.x1*rod.x1+rod.x2*rod.x2);
|
|
||||||
oaf.feed(rod.e1);
|
|
||||||
//.feed(rod.x1*rod.x1/2/time);
|
|
||||||
//empyy.feed(rod.x2*rod.x2/2/time);
|
|
||||||
if (meanCount % 100000 == 0) {
|
|
||||||
//std::cout<<msd.getSEM()<<" "<<std::abs(targetMSD-msd.getMean());
|
|
||||||
std::cout << "\t" << oaf.getSEM() << " " << std::abs(targetOAF - oaf.getMean()) << std::endl;
|
|
||||||
//std::cout<<"\t"<<empxx.getSEM()<<" "<<std::abs(targetEMPXX-empxx.getMean());
|
|
||||||
//std::cout<<"\t"<<empyy.getSEM()<<" "<<std::abs(targetEMPYY-empyy.getMean())<<std::endl;;
|
|
||||||
}
|
|
||||||
if (meanCount % 1000000 == 0) {
|
|
||||||
std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
|
|
||||||
fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
|
|
||||||
fileOAF.close();
|
|
||||||
}
|
|
||||||
meanCount++;
|
|
||||||
}
|
|
||||||
std::ofstream fileOAF(filename + "_length=" + std::to_string(L) + "_dt=" + std::to_string(rod.deltaT) + "_OAF.out");
|
|
||||||
fileOAF << rod.deltaT << "\t" << std::abs(targetOAF - oaf.getMean()) << "\t" << oaf.getSEM() << "\n";
|
|
||||||
fileOAF.close();
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
#endif// SIMULATION2D_HPP
|
|
@ -1,72 +0,0 @@
|
|||||||
#include <random>
|
|
||||||
#include <span>
|
|
||||||
#include <string>
|
|
||||||
#include "legacy/Integratoren.hpp"
|
|
||||||
#include "legacy/Simulation2d.hpp"
|
|
||||||
const int seed = 12345;
|
|
||||||
[[maybe_unused]] void overTime(int sim)
|
|
||||||
{
|
|
||||||
const std::vector<double> Length = { 1.0, 1.3974, 20.0 };
|
|
||||||
const int numCases = 5;
|
|
||||||
int i = 0;
|
|
||||||
for (auto L : Length) {
|
|
||||||
if (sim == numCases * i) {
|
|
||||||
simulate<Integratoren::integrateITOEuler_1>(L, seed, "1");
|
|
||||||
}
|
|
||||||
if (sim == numCases * i + 1) {
|
|
||||||
simulate<Integratoren::integrateITOHeun_2>(L, seed, "2");
|
|
||||||
}
|
|
||||||
if (sim == numCases * i + 2) {
|
|
||||||
simulate<Integratoren::integrateExact_3>(L, seed, "3");
|
|
||||||
}
|
|
||||||
if (sim == numCases * i + 3) {
|
|
||||||
simulate<Integratoren::integrateBDAS_4>(L, seed, "4");
|
|
||||||
}
|
|
||||||
if (sim == numCases * i + 4) {
|
|
||||||
simulate<Integratoren::integrateMBD_5>(L, seed, "5");
|
|
||||||
}
|
|
||||||
i++;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
[[maybe_unused]]void konvergenz(int sim)
|
|
||||||
{
|
|
||||||
const std::vector<double> Length = { 1.0, 1.3974, 20.0 };
|
|
||||||
const std::vector<double> DELTAT = { 2.0, 1.0, 0.5, 0.25, 0.2, 0.1, 0.05, 0.025, 0.02, 0.01, 0.005, 0.0025, 0.200, 0.001, 0.0005, 0.00025 };
|
|
||||||
int i = 0;
|
|
||||||
for (auto dt : DELTAT) {
|
|
||||||
|
|
||||||
for (auto L : Length) {
|
|
||||||
if (sim == i++) {
|
|
||||||
konvergenz<Integratoren::integrateITOEuler_1>(L, seed, dt, "1");
|
|
||||||
}
|
|
||||||
if (sim == i++) {
|
|
||||||
konvergenz<Integratoren::integrateITOHeun_2>(L, seed, dt, "2");
|
|
||||||
}
|
|
||||||
if (sim == i++) {
|
|
||||||
konvergenz<Integratoren::integrateExact_3>(L, seed, dt, "3");
|
|
||||||
}
|
|
||||||
if (sim == i++) {
|
|
||||||
konvergenz<Integratoren::integrateBDAS_4>(L, seed, dt, "4");
|
|
||||||
}
|
|
||||||
if (sim == i++) {
|
|
||||||
konvergenz<Integratoren::integrateMBD_5>(L, seed, dt, "5");
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
int main(int argc, char *argv[])
|
|
||||||
{
|
|
||||||
int sim;
|
|
||||||
auto args = std::span(argv, size_t(argc));
|
|
||||||
if (argc > 1) {
|
|
||||||
sim = std::stoi(args[1]) - 1;
|
|
||||||
} else {
|
|
||||||
sim = 0;
|
|
||||||
}
|
|
||||||
overTime(sim);
|
|
||||||
//konvergenz(sim);
|
|
||||||
return 0;
|
|
||||||
}
|
|
76
src/main.cpp
76
src/main.cpp
@ -3,49 +3,49 @@
|
|||||||
//
|
//
|
||||||
|
|
||||||
#include <iostream>
|
#include <iostream>
|
||||||
|
|
||||||
#include "Calculation.h"
|
#include "Calculation.h"
|
||||||
#include "Integratoren2d_forceless.h"
|
|
||||||
#include "Integratoren2d_force.h"
|
#include "Integratoren2d_force.h"
|
||||||
|
#include "Integratoren2d_forceless.h"
|
||||||
|
|
||||||
int main() {
|
int main() {
|
||||||
constexpr int numStep = 1000;
|
constexpr int numStep = 10000000;
|
||||||
auto zero_Force = [](const Eigen::Vector2d & /*pos*/,
|
constexpr double k = 0.01;
|
||||||
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; };
|
[[maybe_unused]] auto harmonic_Force =
|
||||||
auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
|
[](const Eigen::Vector2d &pos,
|
||||||
const Eigen::Vector2d & /*torque*/) -> double { return 0.0; };
|
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
|
||||||
/* auto euler_Zero = [&](Rod2d &rod, Simulation &sim) {
|
return -k * pos;
|
||||||
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);
|
|
||||||
};
|
};
|
||||||
/*auto exact_Zero = [&](Rod2d &rod, Simulation &sim) {
|
[[maybe_unused]] auto zero_Force =
|
||||||
return Integratoren2d_force::Set3_Exact(rod, sim, zero_Force, zero_Torque);
|
[](const Eigen::Vector2d & /*pos*/,
|
||||||
};*//*
|
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
|
||||||
auto bdas_Zero = [&](Rod2d &rod, Simulation &sim) {
|
return {0.0, 0.0};
|
||||||
return Integratoren2d_force::Set4_BDAS(rod, sim, zero_Force, zero_Torque);
|
};
|
||||||
};*/
|
[[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
|
||||||
Calculation euler(heun_Zero, {{Compute::Type::msd, 1},
|
const Eigen::Vector2d & /*torque*/) -> double {
|
||||||
{Compute::Type::msd, 2},
|
return 0.0;
|
||||||
{Compute::Type::msd, 4},
|
};
|
||||||
{Compute::Type::msd, 8},
|
|
||||||
{Compute::Type::msd, 16},
|
constexpr Compute::Type compute = Compute::Type::oaf;
|
||||||
{Compute::Type::oaf, 1},
|
Calculation euler(Integratoren2d_force::Set2_Heun,
|
||||||
{Compute::Type::oaf, 2},
|
{{compute, 1},
|
||||||
{Compute::Type::oaf, 4},
|
{compute, 2},
|
||||||
{Compute::Type::oaf, 8},
|
{compute, 4},
|
||||||
{Compute::Type::oaf, 16}}, 0.01, 12345);
|
{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);
|
euler.run(numStep);
|
||||||
for (const auto &com: euler.getComputes()) {
|
for (const auto &com : euler.getComputes()) {
|
||||||
if (com.getType() != Compute::Type::msd)
|
if (com.getType() != compute) { continue;
|
||||||
continue;
|
}
|
||||||
std::cout << "MSD: " << com.getDifference() << " " << com.getAgg().getMean()
|
std::cout
|
||||||
<< " <-> " << com.getTarget() << std::endl;
|
<< 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;
|
|
||||||
}
|
}
|
||||||
}
|
}
|
13
src/vec_trafos.h
Normal file
13
src/vec_trafos.h
Normal file
@ -0,0 +1,13 @@
|
|||||||
|
//
|
||||||
|
// Created by jholder on 10/25/21.
|
||||||
|
//
|
||||||
|
#pragma once
|
||||||
|
#include <Eigen/Dense>
|
||||||
|
inline static Eigen::Vector2d ortho(Eigen::Vector2d e) {
|
||||||
|
return Eigen::Vector2d{-e[1], e[0]};
|
||||||
|
}
|
||||||
|
inline static Eigen::Matrix2d rotation_Matrix(Eigen::Vector2d e) {
|
||||||
|
Eigen::Matrix2d mat;
|
||||||
|
mat << e[0], -e[1], e[1], e[0];
|
||||||
|
return mat;
|
||||||
|
}
|
@ -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 <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
|
|
||||||
|
|
||||||
|
@ -2,32 +2,41 @@
|
|||||||
// Created by jholder on 22.10.21.
|
// Created by jholder on 22.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
|
#include "Calculation.h"
|
||||||
|
|
||||||
#include <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
#include <Calculation.h>
|
|
||||||
#include "LiveAgg.hpp"
|
|
||||||
#include "Integratoren2d_forceless.h"
|
#include "Integratoren2d_forceless.h"
|
||||||
|
#include "LiveAgg.hpp"
|
||||||
|
|
||||||
TEST_CASE("Basic run of Calculation") {
|
TEST_CASE("Basic run of Calculation") {
|
||||||
Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 10},
|
Calculation euler(Integratoren2d_forceless::Set1_Euler,
|
||||||
{Compute::Type::oaf, 10}}, 1, 1);
|
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
|
||||||
Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 10},
|
1);
|
||||||
{Compute::Type::oaf, 10}}, 1, 1);
|
Calculation heun(Integratoren2d_forceless::Set2_Heun,
|
||||||
Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 10},
|
{{Compute::Type::msd, 10}, {Compute::Type::oaf, 10}}, 1,
|
||||||
{Compute::Type::oaf, 10}}, 1, 1);
|
1);
|
||||||
Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 10},
|
Calculation exact(Integratoren2d_forceless::Set3_Exact,
|
||||||
{Compute::Type::oaf, 10}}, 1, 1);
|
{{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") {
|
SECTION("Euler") {
|
||||||
euler.run(100);
|
euler.run(100);
|
||||||
CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
CHECK(euler.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
||||||
CHECK(euler.getComputes()[0].getType() == Compute::Type::msd);
|
CHECK(euler.getComputes()[0].getType() == Compute::Type::msd);
|
||||||
CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf);
|
CHECK(euler.getComputes()[1].getType() == Compute::Type::oaf);
|
||||||
}SECTION("Heun") {
|
}
|
||||||
|
SECTION("Heun") {
|
||||||
heun.run(100);
|
heun.run(100);
|
||||||
CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
CHECK(heun.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
||||||
}SECTION("Exact") {
|
}
|
||||||
|
SECTION("Exact") {
|
||||||
exact.run(100);
|
exact.run(100);
|
||||||
CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
CHECK(exact.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
||||||
}SECTION("Euler") {
|
}
|
||||||
|
SECTION("Euler") {
|
||||||
bdas.run(100);
|
bdas.run(100);
|
||||||
CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
CHECK(bdas.getRod().getE().norm() == Catch::Detail::Approx(1.0));
|
||||||
}
|
}
|
||||||
@ -36,17 +45,20 @@ TEST_CASE("Basic run of Calculation") {
|
|||||||
euler.run(100);
|
euler.run(100);
|
||||||
CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10);
|
CHECK(euler.getComputes()[0].getAgg().getNumPoints() == 10);
|
||||||
CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10);
|
CHECK(euler.getComputes()[1].getAgg().getNumPoints() == 10);
|
||||||
}SECTION("Deterministic") {
|
}
|
||||||
|
SECTION("Deterministic") {
|
||||||
Calculation determ(Integratoren2d_forceless::Set0_deterministic,
|
Calculation determ(Integratoren2d_forceless::Set0_deterministic,
|
||||||
{{Compute::Type::msd, 1},
|
{{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}},
|
||||||
{Compute::Type::oaf, 1}}, 1, 1);
|
1, 1);
|
||||||
determ.run(10);
|
determ.run(10);
|
||||||
auto time = 1;
|
auto time = 1;
|
||||||
|
|
||||||
auto targetMSD = pow(0.01, 2) * time;
|
auto targetMSD = pow(0.01, 2) * time;
|
||||||
auto targetOAF = cos(0.1);
|
auto targetOAF = cos(0.1);
|
||||||
|
|
||||||
CHECK(determ.getComputes()[0].getAgg().getMean() == Catch::Detail::Approx(targetMSD));
|
CHECK(determ.getComputes()[0].getAgg().getMean() ==
|
||||||
CHECK(determ.getComputes()[1].getAgg().getMean() == Catch::Detail::Approx(targetOAF));
|
Catch::Detail::Approx(targetMSD));
|
||||||
|
CHECK(determ.getComputes()[1].getAgg().getMean() ==
|
||||||
|
Catch::Detail::Approx(targetOAF));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -2,19 +2,20 @@
|
|||||||
// Created by jholder on 22.10.21.
|
// Created by jholder on 22.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
#include <Compute.h>
|
#include "Compute.h"
|
||||||
#include "catch2/catch.hpp"
|
|
||||||
|
#include <catch2/catch.hpp>
|
||||||
#include "Integratoren2d_forceless.h"
|
#include "Integratoren2d_forceless.h"
|
||||||
|
|
||||||
TEST_CASE("Compute") {
|
TEST_CASE("Compute") {
|
||||||
Rod2d rod(1.0);
|
Rod2d rod(1.0);
|
||||||
Simulation sim(0.1,1);
|
Simulation sim(0.1, 1);
|
||||||
auto com = Compute(rod, Compute::Type::msd, 10, sim);
|
auto com = Compute(rod, Compute::Type::msd, 10, sim);
|
||||||
SECTION("Mean of same values") {
|
SECTION("Mean of same values") {
|
||||||
for (int i = 0; i < 100; ++i) {
|
for (int i = 0; i < 100; ++i) {
|
||||||
com.eval(rod);
|
com.eval(rod);
|
||||||
}
|
}
|
||||||
CHECK(com.getAgg().getNumPoints()==10);
|
CHECK(com.getAgg().getNumPoints() == 10);
|
||||||
CHECK(com.getAgg().getMean()==0.0);
|
CHECK(com.getAgg().getMean() == 0.0);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -1,89 +1,357 @@
|
|||||||
//
|
//
|
||||||
// Created by jholder on 24.10.21.
|
// Created by jholder on 24.10.21.
|
||||||
//
|
//
|
||||||
|
#include "Calculation.h"
|
||||||
|
#include "Compute.h"
|
||||||
|
#include "Integratoren2d_force.h"
|
||||||
|
#include "Integratoren2d_forceless.h"
|
||||||
|
|
||||||
#include <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
#include <Integratoren2d_forceless.h>
|
|
||||||
#include <Calculation.h>
|
|
||||||
#include <Compute.h>
|
|
||||||
#include <Integratoren2d_force.h>
|
|
||||||
|
|
||||||
TEST_CASE("Forceless Integratoren") {
|
TEST_CASE("Forceless Integratoren") {
|
||||||
const size_t SEED = Catch::rngSeed();
|
const size_t SEED = Catch::rngSeed();
|
||||||
auto force = [](const Eigen::Vector2d & /*pos*/,
|
const double length = GENERATE(1.0,1.4,2.0);
|
||||||
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d { return {0.0, 0.0}; };
|
constexpr int steps = 10000;
|
||||||
auto torque = [](const Eigen::Vector2d & /*pos*/,
|
constexpr double delta = 0.1;
|
||||||
const Eigen::Vector2d & /*torque*/) -> double { return 0.0; };
|
|
||||||
|
|
||||||
|
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,length);
|
||||||
|
euler.run(steps);
|
||||||
|
|
||||||
|
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,length);
|
||||||
|
heun.run(steps);
|
||||||
|
|
||||||
|
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,length);
|
||||||
|
exact.run(steps);
|
||||||
|
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,length);
|
||||||
|
bdas.run(steps);
|
||||||
|
Calculation mbd(Integratoren2d_forceless::Set5_MBD,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,length);
|
||||||
|
mbd.run(steps);
|
||||||
|
SECTION("ForceLess Euler") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: euler.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("ForceLess Heun") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: heun.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("ForceLess Exact") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: exact.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("ForceLess BDAS") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: bdas.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("ForceLess MBD") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: mbd.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
[[maybe_unused]] auto zero_Force =
|
||||||
|
[](const Eigen::Vector2d & /*pos*/,
|
||||||
|
const Eigen::Vector2d & /*torque*/) -> Eigen::Vector2d {
|
||||||
|
return {0.0, 0.0};
|
||||||
|
};
|
||||||
|
[[maybe_unused]] auto zero_Torque = [](const Eigen::Vector2d & /*pos*/,
|
||||||
|
const Eigen::Vector2d & /*torque*/) -> double {
|
||||||
|
return 0.0;
|
||||||
|
};
|
||||||
|
Calculation euler_zero(Integratoren2d_force::Set1_Euler,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,zero_Force,zero_Torque,length);
|
||||||
|
euler_zero.run(steps);
|
||||||
|
|
||||||
|
Calculation heun_zero(Integratoren2d_force::Set2_Heun,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,zero_Force,zero_Torque,length);
|
||||||
|
heun_zero.run(steps);
|
||||||
|
|
||||||
|
Calculation exact_zero(Integratoren2d_force::Set3_Exact,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,zero_Force,zero_Torque,length);
|
||||||
|
exact_zero.run(steps);
|
||||||
|
|
||||||
|
Calculation bdas_zero(Integratoren2d_force::Set4_BDAS,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,zero_Force,zero_Torque,length);
|
||||||
|
bdas_zero.run(steps);
|
||||||
|
|
||||||
|
Calculation mbd_zero(Integratoren2d_force::Set5_MBD,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED,zero_Force,zero_Torque,length);
|
||||||
|
mbd_zero.run(steps);
|
||||||
|
SECTION("Force Euler") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: euler_zero.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("Force Heun") {
|
||||||
|
CAPTURE(length);
|
||||||
|
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: heun_zero.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("Force Exact") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: exact_zero.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("Force BDAS") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: bdas_zero.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
SECTION("Force MBD") {
|
||||||
|
CAPTURE(length);
|
||||||
|
size_t i = 0;
|
||||||
|
for (auto &c: mbd_zero.getComputes()) {
|
||||||
|
CAPTURE(c);
|
||||||
|
CHECK(c.getDifference() <= delta*c.getAgg().getMean());
|
||||||
|
++i;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/*
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
SECTION("Euler") {
|
SECTION("Euler") {
|
||||||
|
|
||||||
Calculation euler(Integratoren2d_forceless::Set1_Euler, {{Compute::Type::msd, 1},
|
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED);
|
|
||||||
{
|
{
|
||||||
euler.run(10000);
|
|
||||||
for (auto &c: euler.getComputes()) { CHECK(c.getDifference() <= 0.005); }
|
|
||||||
}
|
}
|
||||||
|
|
||||||
Calculation euler_force(Integratoren2d_force::Set1_Euler, {
|
Calculation euler_force(
|
||||||
{Compute::Type::msd, 1},
|
Integratoren2d_force::Set1_Euler,
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED, force, torque);
|
{{Compute::Type::msd, 1}, {Compute::Type::oaf, 1}}, 0.01, SEED,
|
||||||
|
force, torque);
|
||||||
{
|
{
|
||||||
euler_force.run(10000);
|
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()[0].getAgg().getMean() ==
|
||||||
CHECK(euler.getComputes()[1].getAgg().getMean() == euler_force.getComputes()[1].getAgg().getMean());
|
euler_force.getComputes()[0].getAgg().getMean());
|
||||||
}SECTION("Heun") {
|
CHECK(euler.getComputes()[1].getAgg().getMean() ==
|
||||||
|
euler_force.getComputes()[1].getAgg().getMean());
|
||||||
Calculation heun(Integratoren2d_forceless::Set2_Heun, {{Compute::Type::msd, 1},
|
}
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED);
|
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);
|
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},
|
Calculation heun_force(Integratoren2d_force::Set2_Heun,
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED, force,
|
{{Compute::Type::msd, 1},
|
||||||
torque);
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED, force, torque);
|
||||||
{
|
{
|
||||||
heun_force.run(10000);
|
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()[0].getAgg().getMean() ==
|
||||||
CHECK(heun.getComputes()[1].getAgg().getMean() == heun_force.getComputes()[1].getAgg().getMean());
|
heun_force.getComputes()[0].getAgg().getMean());
|
||||||
}SECTION("Exact") {
|
CHECK(heun.getComputes()[1].getAgg().getMean() ==
|
||||||
|
heun_force.getComputes()[1].getAgg().getMean());
|
||||||
Calculation exact(Integratoren2d_forceless::Set3_Exact, {{Compute::Type::msd, 1},
|
}
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED);
|
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);
|
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},
|
Calculation exact_force(Integratoren2d_force::Set3_Exact,
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED, force,
|
{{Compute::Type::msd, 1},
|
||||||
torque);
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED, force, torque);
|
||||||
{
|
{
|
||||||
exact_force.run(10000);
|
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()[0].getAgg().getMean() ==
|
||||||
CHECK(exact.getComputes()[1].getAgg().getMean() == exact_force.getComputes()[1].getAgg().getMean());
|
Approx(exact_force.getComputes()[0].getAgg().getMean())
|
||||||
}SECTION("BDAS") {
|
.epsilon(1e-10));
|
||||||
|
CHECK(exact.getComputes()[1].getAgg().getMean() ==
|
||||||
Calculation bdas(Integratoren2d_forceless::Set4_BDAS, {{Compute::Type::msd, 1},
|
exact_force.getComputes()[1].getAgg().getMean());
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED);
|
}
|
||||||
|
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);
|
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},
|
Calculation bdas_force(Integratoren2d_force::Set4_BDAS,
|
||||||
{Compute::Type::oaf, 1}}, 0.01, SEED, force,
|
{{Compute::Type::msd, 1},
|
||||||
torque);
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED, force, torque);
|
||||||
{
|
{
|
||||||
bdas_force.run(10000);
|
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() == bdas_force.getComputes()[0].getAgg().getMean());
|
CHECK(bdas.getComputes()[0].getAgg().getMean() ==
|
||||||
CHECK(bdas.getComputes()[1].getAgg().getMean() == bdas_force.getComputes()[1].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));
|
||||||
|
}
|
||||||
|
|
||||||
|
SECTION("MBD") {
|
||||||
|
Calculation mbd(Integratoren2d_forceless::Set5_MBD,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED);
|
||||||
|
{
|
||||||
|
mbd.run(10000);
|
||||||
|
for (auto &c : mbd.getComputes()) {
|
||||||
|
CHECK(c.getDifference() <= 0.005);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
Calculation mbd_force(Integratoren2d_force::Set5_MBD,
|
||||||
|
{{Compute::Type::msd, 1},
|
||||||
|
{Compute::Type::oaf, 1},
|
||||||
|
{Compute::Type::empxx, 1},
|
||||||
|
{Compute::Type::empyy, 1}},
|
||||||
|
0.01, SEED, force, torque);
|
||||||
|
{
|
||||||
|
mbd_force.run(10000);
|
||||||
|
for (auto &c : mbd_force.getComputes()) {
|
||||||
|
CHECK(c.getDifference() <= 0.005);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
CHECK(mbd.getComputes()[0].getAgg().getMean() ==
|
||||||
|
Approx(mbd_force.getComputes()[0].getAgg().getMean())
|
||||||
|
.epsilon(10e-10));
|
||||||
|
CHECK(mbd.getComputes()[1].getAgg().getMean() ==
|
||||||
|
Approx(mbd_force.getComputes()[1].getAgg().getMean())
|
||||||
|
.epsilon(10e-10));
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
*/
|
@ -3,14 +3,15 @@
|
|||||||
//
|
//
|
||||||
|
|
||||||
#include <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
|
|
||||||
#include "LiveAgg.hpp"
|
#include "LiveAgg.hpp"
|
||||||
|
|
||||||
TEST_CASE("LiveAgg") {
|
TEST_CASE("LiveAgg") {
|
||||||
auto ag = LiveAgg();
|
auto ag = LiveAgg();
|
||||||
SECTION("Mean of same values"){
|
SECTION("Mean of same values") {
|
||||||
ag.feed(1.0);
|
ag.feed(1.0);
|
||||||
ag.feed(1.0);
|
ag.feed(1.0);
|
||||||
REQUIRE(ag.getMean()==1.0);
|
REQUIRE(ag.getMean() == 1.0);
|
||||||
REQUIRE(ag.getNumPoints()==2);
|
REQUIRE(ag.getNumPoints() == 2);
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
@ -2,10 +2,11 @@
|
|||||||
// Created by jholder on 21.10.21.
|
// Created by jholder on 21.10.21.
|
||||||
//
|
//
|
||||||
|
|
||||||
#include "Rod2d.hpp"
|
|
||||||
#include <catch2/catch.hpp>
|
|
||||||
#include <Eigen/Dense>
|
#include <Eigen/Dense>
|
||||||
|
#include <catch2/catch.hpp>
|
||||||
|
|
||||||
|
#include "Rod2d.hpp"
|
||||||
|
#include "vec_trafos.h"
|
||||||
TEST_CASE("Rods") {
|
TEST_CASE("Rods") {
|
||||||
Rod2d sphere(1);
|
Rod2d sphere(1);
|
||||||
Rod2d rod(2);
|
Rod2d rod(2);
|
||||||
@ -14,34 +15,19 @@ TEST_CASE("Rods") {
|
|||||||
auto newPos = rod.getPos();
|
auto newPos = rod.getPos();
|
||||||
REQUIRE(newPos[0] == 1);
|
REQUIRE(newPos[0] == 1);
|
||||||
REQUIRE(newPos[1] == -1);
|
REQUIRE(newPos[1] == -1);
|
||||||
REQUIRE(sphere.getDiff_Sqrt()(0,0) == 1.0);
|
REQUIRE(sphere.getDiff_Sqrt()(0, 0) == 1.0);
|
||||||
REQUIRE(sphere.getDiff_Sqrt()(0,1) == 0.0);
|
REQUIRE(sphere.getDiff_Sqrt()(0, 1) == 0.0);
|
||||||
REQUIRE(sphere.getDiff_Sqrt()(1,0) == 0.0);
|
REQUIRE(sphere.getDiff_Sqrt()(1, 0) == 0.0);
|
||||||
REQUIRE(sphere.getDiff_Sqrt()(1,1) == 1.0);
|
REQUIRE(sphere.getDiff_Sqrt()(1, 1) == 1.0);
|
||||||
}SECTION("Checking rotation") {
|
}
|
||||||
SECTION("Forwards == e 90"){
|
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());
|
||||||
rod.setE({0, 1}); // Lab system is 90 degree rotated
|
auto test = Eigen::Vector2d({1, 0});
|
||||||
auto test = Eigen::Vector2d ({1,0});
|
|
||||||
auto erg = (rod.getE_Base_matrix()*test).normalized();
|
|
||||||
auto e = rod.getE();
|
auto e = rod.getE();
|
||||||
|
auto erg = (rotation_Matrix(e) * test).normalized();
|
||||||
REQUIRE(erg.isApprox(e));
|
REQUIRE(erg.isApprox(e));
|
||||||
|
|
||||||
}
|
|
||||||
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();
|
|
||||||
auto e = rod.getE();
|
|
||||||
REQUIRE(erg.isApprox(e));
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
}
|
}
|
||||||
|
|
||||||
}
|
}
|
@ -3,25 +3,26 @@
|
|||||||
//
|
//
|
||||||
|
|
||||||
#include <catch2/catch.hpp>
|
#include <catch2/catch.hpp>
|
||||||
|
|
||||||
#include "Simulation.h"
|
#include "Simulation.h"
|
||||||
|
|
||||||
TEST_CASE("Simulation") {
|
TEST_CASE("Simulation") {
|
||||||
auto sim = Simulation(1,0);
|
auto sim = Simulation(1, 0);
|
||||||
REQUIRE(sim.getSTD()==sqrt(1*2.0));
|
REQUIRE(sim.getSTD() == sqrt(1 * 2.0));
|
||||||
SECTION("Symmetrie"){
|
SECTION("Symmetrie") {
|
||||||
constexpr int num = 10000000;
|
constexpr int num = 10000000;
|
||||||
double sum = 0.0;
|
double sum = 0.0;
|
||||||
for (int i = 0; i < num; ++i) {
|
for (int i = 0; i < num; ++i) {
|
||||||
sum += sim.getNorm(1);
|
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;
|
constexpr int num = 10000000;
|
||||||
double sum = 0.0;
|
double sum = 0.0;
|
||||||
for (int i = 0; i < num; ++i) {
|
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));
|
||||||
}
|
}
|
||||||
}
|
}
|
Loading…
Reference in New Issue
Block a user