diff --git a/src/Makefile.am b/src/Makefile.am index 45d9248544e41fc971b571847c8d44a5542bb25f..42c6012a19d6a29665eee84418567c8a4c086a53 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -5,6 +5,7 @@ common_sources = \ boundary_writer.cc \ friction_writer.cc \ solverfactory.cc \ + state.cc \ timestepping.cc \ vtk.cc diff --git a/src/state.cc b/src/state.cc new file mode 100644 index 0000000000000000000000000000000000000000..49567d3f2b112e2ba72af6870f5458e9f9cf6ecf --- /dev/null +++ b/src/state.cc @@ -0,0 +1,25 @@ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif + +#include "state.hh" +#include "state/ageinglawstateupdater.cc" +#include "state/sliplawstateupdater.cc" + +template <class ScalarVector, class Vector> +std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater( + Config::stateModel model, ScalarVector const &alpha_initial, + Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) { + switch (model) { + case Config::AgeingLaw: + return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>( + alpha_initial, frictionalNodes, L, V0); + case Config::SlipLaw: + return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>( + alpha_initial, frictionalNodes, L, V0); + default: + assert(false); + } +} + +#include "state_tmpl.cc" diff --git a/src/state.hh b/src/state.hh index cfa60f94021780f9394d218bf036ad540e910925..5726ef3b02dbd97d1d7a3c12c80be32a41ac77b0 100644 --- a/src/state.hh +++ b/src/state.hh @@ -6,6 +6,7 @@ #include <dune/common/bitsetvector.hh> #include "enums.hh" + #include "state/stateupdater.hh" #include "state/ageinglawstateupdater.hh" #include "state/sliplawstateupdater.hh" @@ -13,16 +14,6 @@ template <class ScalarVector, class Vector> std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater( Config::stateModel model, ScalarVector const &alpha_initial, - Dune::BitSetVector<1> const &frictionalNodes, double L, double V0) { - switch (model) { - case Config::AgeingLaw: - return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>( - alpha_initial, frictionalNodes, L, V0); - case Config::SlipLaw: - return std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>( - alpha_initial, frictionalNodes, L, V0); - default: - assert(false); - } -} + Dune::BitSetVector<1> const &frictionalNodes, double L, double V0); + #endif diff --git a/src/state/ageinglawstateupdater.cc b/src/state/ageinglawstateupdater.cc new file mode 100644 index 0000000000000000000000000000000000000000..1461e29907e0ac65f6f3a90100a5416088686470 --- /dev/null +++ b/src/state/ageinglawstateupdater.cc @@ -0,0 +1,48 @@ +#include "ageinglawstateupdater.hh" + +template <class ScalarVector, class Vector> +AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater( + ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L, + double _V0) + : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} + +template <class ScalarVector, class Vector> +void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() { + alpha_o = alpha; +} + +template <class ScalarVector, class Vector> +void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) { + tau = _tau; +} + +/* + Compute [ 1-\exp(c*x) ] / x under the assumption that x is + non-negative +*/ +double liftSingularity(double c, double x) { + if (x <= 0) + return -c; + else + return -std::expm1(c * x) / x; +} + +template <class ScalarVector, class Vector> +void AgeingLawStateUpdater<ScalarVector, Vector>::solve( + Vector const &velocity_field) { + for (size_t i = 0; i < nodes.size(); ++i) { + if (not toBool(nodes[i])) + continue; + + double const V = velocity_field[i].two_norm(); + double const mtoL = -tau / L; + alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) + + V0 * liftSingularity(mtoL, V)); + } +} + +template <class ScalarVector, class Vector> +void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha( + ScalarVector &_alpha) { + _alpha = alpha; +} diff --git a/src/state/ageinglawstateupdater.hh b/src/state/ageinglawstateupdater.hh index e7176a1efff6a30845ca7835371a204d796f2b89..5b6fec71f1cbb3411a20f7dd4ca1d867d7d87bb5 100644 --- a/src/state/ageinglawstateupdater.hh +++ b/src/state/ageinglawstateupdater.hh @@ -24,51 +24,4 @@ class AgeingLawStateUpdater : public StateUpdater<ScalarVector, Vector> { double const V0; double tau; }; - -template <class ScalarVector, class Vector> -AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater( - ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L, - double _V0) - : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} - -template <class ScalarVector, class Vector> -void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() { - alpha_o = alpha; -} - -template <class ScalarVector, class Vector> -void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) { - tau = _tau; -} - -/* - Compute [ 1-\exp(c*x) ] / x under the assumption that x is - non-negative -*/ -double liftSingularity(double c, double x) { - if (x <= 0) - return -c; - else - return -std::expm1(c * x) / x; -} - -template <class ScalarVector, class Vector> -void AgeingLawStateUpdater<ScalarVector, Vector>::solve( - Vector const &velocity_field) { - for (size_t i = 0; i < nodes.size(); ++i) { - if (not toBool(nodes[i])) - continue; - - double const V = velocity_field[i].two_norm(); - double const mtoL = -tau / L; - alpha[i] = std::log(std::exp(alpha_o[i] + V * mtoL) + - V0 * liftSingularity(mtoL, V)); - } -} - -template <class ScalarVector, class Vector> -void AgeingLawStateUpdater<ScalarVector, Vector>::extractAlpha( - ScalarVector &_alpha) { - _alpha = alpha; -} #endif diff --git a/src/state/sliplawstateupdater.cc b/src/state/sliplawstateupdater.cc new file mode 100644 index 0000000000000000000000000000000000000000..5397c1d2e4845483dc25843eff816c5978f4b092 --- /dev/null +++ b/src/state/sliplawstateupdater.cc @@ -0,0 +1,37 @@ +#include "sliplawstateupdater.hh" + +template <class ScalarVector, class Vector> +SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater( + ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L, + double _V0) + : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} + +template <class ScalarVector, class Vector> +void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() { + alpha_o = alpha; +} + +template <class ScalarVector, class Vector> +void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) { + tau = _tau; +} + +template <class ScalarVector, class Vector> +void SlipLawStateUpdater<ScalarVector, Vector>::solve( + Vector const &velocity_field) { + for (size_t i = 0; i < nodes.size(); ++i) { + if (not toBool(nodes[i])) + continue; + + double const V = velocity_field[i].two_norm(); + double const mtVoL = -tau * V / L; + alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) + + alpha_o[i] * std::exp(mtVoL); + } +} + +template <class ScalarVector, class Vector> +void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha( + ScalarVector &_alpha) { + _alpha = alpha; +} diff --git a/src/state/sliplawstateupdater.hh b/src/state/sliplawstateupdater.hh index 9bb82210bb3f6769688547d0251b713e0187b3f5..f527c1c40ddd16714a3b79c03a3be59dfde1a1b1 100644 --- a/src/state/sliplawstateupdater.hh +++ b/src/state/sliplawstateupdater.hh @@ -25,40 +25,4 @@ class SlipLawStateUpdater : public StateUpdater<ScalarVector, Vector> { double tau; }; -template <class ScalarVector, class Vector> -SlipLawStateUpdater<ScalarVector, Vector>::SlipLawStateUpdater( - ScalarVector _alpha_initial, Dune::BitSetVector<1> const &_nodes, double _L, - double _V0) - : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {} - -template <class ScalarVector, class Vector> -void SlipLawStateUpdater<ScalarVector, Vector>::nextTimeStep() { - alpha_o = alpha; -} - -template <class ScalarVector, class Vector> -void SlipLawStateUpdater<ScalarVector, Vector>::setup(double _tau) { - tau = _tau; -} - -template <class ScalarVector, class Vector> -void SlipLawStateUpdater<ScalarVector, Vector>::solve( - Vector const &velocity_field) { - for (size_t i = 0; i < nodes.size(); ++i) { - if (not toBool(nodes[i])) - continue; - - double const V = velocity_field[i].two_norm(); - double const mtVoL = -tau * V / L; - alpha[i] = (V <= 0) ? alpha_o[i] : std::expm1(mtVoL) * std::log(V / V0) + - alpha_o[i] * std::exp(mtVoL); - } -} - -template <class ScalarVector, class Vector> -void SlipLawStateUpdater<ScalarVector, Vector>::extractAlpha( - ScalarVector &_alpha) { - _alpha = alpha; -} - #endif diff --git a/src/state_tmpl.cc b/src/state_tmpl.cc new file mode 100644 index 0000000000000000000000000000000000000000..65d7ba9c77e7e70211ebf5c3f79ad75cb75181a3 --- /dev/null +++ b/src/state_tmpl.cc @@ -0,0 +1,7 @@ +#include "explicitvectors.hh" + +template std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater< + ScalarVector, Vector>(Config::stateModel model, + ScalarVector const &alpha_initial, + Dune::BitSetVector<1> const &frictionalNodes, + double L, double V0);