diff --git a/src/dieterichstateupdater.hh b/src/dieterichstateupdater.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ee8480c473e4f1f32802a97e51a2d9b21ef8d839
--- /dev/null
+++ b/src/dieterichstateupdater.hh
@@ -0,0 +1,60 @@
+#ifndef DIETERICH_STATE_UPDATER_HH
+#define DIETERICH_STATE_UPDATER_HH
+
+#include "compute_state_dieterich_euler.hh"
+#include "stateupdater.hh"
+
+template <class SingletonVectorType, class VectorType>
+class DieterichStateUpdater
+ : public StateUpdater<SingletonVectorType, VectorType> {
+public:
+ DieterichStateUpdater(SingletonVectorType _alpha_initial,
+ Dune::BitSetVector<1> const &_nodes, double _L);
+
+ virtual void nextTimeStep();
+ virtual void setup(double _tau);
+ virtual void solve(VectorType const &velocity_field);
+ virtual void extractState(SingletonVectorType &state);
+
+private:
+ SingletonVectorType alpha_old;
+ SingletonVectorType alpha;
+ Dune::BitSetVector<1> const &nodes;
+ double const L;
+ double tau;
+};
+
+template <class SingletonVectorType, class VectorType>
+DieterichStateUpdater<SingletonVectorType, VectorType>::DieterichStateUpdater(
+ SingletonVectorType _alpha_initial, Dune::BitSetVector<1> const &_nodes,
+ double _L)
+ : alpha(_alpha_initial), nodes(_nodes), L(_L) {}
+
+template <class SingletonVectorType, class VectorType>
+void DieterichStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
+ alpha_old = alpha;
+}
+
+template <class SingletonVectorType, class VectorType>
+void DieterichStateUpdater<SingletonVectorType, VectorType>::setup(
+ double _tau) {
+ tau = _tau;
+}
+
+template <class SingletonVectorType, class VectorType>
+void DieterichStateUpdater<SingletonVectorType, VectorType>::solve(
+ VectorType const &velocity_field) {
+ for (size_t i = 0; i < nodes.size(); ++i)
+ if (nodes[i][0]) {
+ double const V = velocity_field[i].two_norm();
+ alpha[i] = state_update_dieterich_euler(tau, V / L, alpha_old[i]);
+ }
+}
+
+template <class SingletonVectorType, class VectorType>
+void DieterichStateUpdater<SingletonVectorType, VectorType>::extractState(
+ SingletonVectorType &state) {
+ state = alpha;
+}
+
+#endif
diff --git a/src/one-body-sample.cc b/src/one-body-sample.cc
index 36827736004729f89784a3129f477e6e7fa84c2f..8be52082d744cc99943999bb3e5e75b3c4bc020f 100644
--- a/src/one-body-sample.cc
+++ b/src/one-body-sample.cc
@@ -65,8 +65,6 @@
#include <dune/tectonic/myconvexproblem.hh>
#include "assemblers.hh"
-#include "compute_state_dieterich_euler.hh"
-#include "compute_state_ruina.hh"
#include "mysolver.hh"
#include "vtk.hh"
@@ -77,6 +75,9 @@
#include "timestepping.hh"
+#include "ruinastateupdater.hh"
+#include "dieterichstateupdater.hh"
+
int const dim = DIM;
template <class GridView, class GridCorner>
@@ -284,9 +285,9 @@ int main(int argc, char *argv[]) {
VectorType udd_initial(finestSize);
udd_initial = 0.0;
- SingletonVectorType alpha_old(finestSize);
- alpha_old = parset.get<double>("boundary.friction.state.initial");
- SingletonVectorType alpha(alpha_old);
+ SingletonVectorType alpha_initial(finestSize);
+ alpha_initial = parset.get<double>("boundary.friction.state.initial");
+ SingletonVectorType alpha(alpha_initial);
SingletonVectorType vonMisesStress;
@@ -355,11 +356,28 @@ int main(int argc, char *argv[]) {
break;
}
+ Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType>>
+ stateUpdater;
+ switch (parset.get<Config::state_model>("boundary.friction.state.model")) {
+ case Config::Dieterich:
+ stateUpdater = Dune::make_shared<
+ DieterichStateUpdater<SingletonVectorType, VectorType>>(
+ alpha_initial, frictionalNodes, L);
+ break;
+ case Config::Ruina:
+ stateUpdater = Dune::make_shared<
+ RuinaStateUpdater<SingletonVectorType, VectorType>>(
+ alpha_initial, frictionalNodes, L);
+ break;
+ }
+
auto const state_fpi_max =
parset.get<size_t>("solver.tnnmg.fixed_point_iterations");
for (size_t run = 1; run <= timesteps; ++run) {
double const time = tau * run;
{
+ stateUpdater->nextTimeStep();
+
VectorType ell(finestSize);
assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
leafView, p1Basis, neumannNodes, ell, neumannFunction, time);
@@ -370,6 +388,7 @@ int main(int argc, char *argv[]) {
VectorType problem_rhs(finestSize);
VectorType problem_iterate(finestSize);
+ stateUpdater->setup(tau);
timeSteppingScheme->setup(ell, tau, time, problem_rhs, problem_iterate,
problem_A);
@@ -396,23 +415,9 @@ int main(int argc, char *argv[]) {
timeSteppingScheme->extractDisplacement(u);
timeSteppingScheme->extractVelocity(ud);
- // Update the state
- for (size_t i = 0; i < frictionalNodes.size(); ++i) {
- if (frictionalNodes[i][0]) {
- double const V = ud[i].two_norm();
-
- switch (parset.get<Config::state_model>(
- "boundary.friction.state.model")) {
- case Config::Dieterich:
- alpha[i] =
- state_update_dieterich_euler(tau, V / L, alpha_old[i]);
- break;
- case Config::Ruina:
- alpha[i] = state_update_ruina(tau, V * tau / L, alpha_old[i]);
- break;
- }
- }
- }
+ stateUpdater->solve(ud);
+ stateUpdater->extractState(alpha);
+
if (parset.get<bool>("printProgress")) {
std::cerr << '.';
std::cerr.flush();
@@ -433,8 +438,6 @@ int main(int argc, char *argv[]) {
std::cerr << std::endl;
}
- alpha_old = alpha;
-
{ // Write all kinds of data
state_writer << alpha[specialNode][0] << " " << std::endl;
diff --git a/src/ruinastateupdater.hh b/src/ruinastateupdater.hh
new file mode 100644
index 0000000000000000000000000000000000000000..083d97eedaeb2d86085a193f7b03b331ea8459a4
--- /dev/null
+++ b/src/ruinastateupdater.hh
@@ -0,0 +1,58 @@
+#ifndef RUINA_STATE_UPDATER_HH
+#define RUINA_STATE_UPDATER_HH
+
+#include "compute_state_ruina.hh"
+#include "stateupdater.hh"
+
+template <class SingletonVectorType, class VectorType>
+class RuinaStateUpdater : public StateUpdater<SingletonVectorType, VectorType> {
+public:
+ RuinaStateUpdater(SingletonVectorType _alpha_initial,
+ Dune::BitSetVector<1> const &_nodes, double _L);
+
+ virtual void nextTimeStep();
+ virtual void setup(double _tau);
+ virtual void solve(VectorType const &velocity_field);
+ virtual void extractState(SingletonVectorType &state);
+
+private:
+ SingletonVectorType alpha_old;
+ SingletonVectorType alpha;
+ Dune::BitSetVector<1> const &nodes;
+ double const L;
+ double tau;
+};
+
+template <class SingletonVectorType, class VectorType>
+RuinaStateUpdater<SingletonVectorType, VectorType>::RuinaStateUpdater(
+ SingletonVectorType _alpha_initial, Dune::BitSetVector<1> const &_nodes,
+ double _L)
+ : alpha(_alpha_initial), nodes(_nodes), L(_L) {}
+
+template <class SingletonVectorType, class VectorType>
+void RuinaStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
+ alpha_old = alpha;
+}
+
+template <class SingletonVectorType, class VectorType>
+void RuinaStateUpdater<SingletonVectorType, VectorType>::setup(double _tau) {
+ tau = _tau;
+}
+
+template <class SingletonVectorType, class VectorType>
+void RuinaStateUpdater<SingletonVectorType, VectorType>::solve(
+ VectorType const &velocity_field) {
+ for (size_t i = 0; i < nodes.size(); ++i)
+ if (nodes[i][0]) {
+ double const V = velocity_field[i].two_norm();
+ alpha[i] = state_update_ruina(tau, V * tau / L, alpha_old[i]);
+ }
+}
+
+template <class SingletonVectorType, class VectorType>
+void RuinaStateUpdater<SingletonVectorType, VectorType>::extractState(
+ SingletonVectorType &state) {
+ state = alpha;
+}
+
+#endif
diff --git a/src/stateupdater.hh b/src/stateupdater.hh
new file mode 100644
index 0000000000000000000000000000000000000000..8c5f412dada1512e8cf6958e845128c22bfa73e9
--- /dev/null
+++ b/src/stateupdater.hh
@@ -0,0 +1,12 @@
+#ifndef STATE_UPDATER_HH
+#define STATE_UPDATER_HH
+
+template <class SingletonVectorType, class VectorType> class StateUpdater {
+public:
+ virtual void nextTimeStep() = 0;
+ virtual void setup(double _tau) = 0;
+ virtual void solve(VectorType const &velocity_field) = 0;
+ virtual void extractState(SingletonVectorType &state) = 0;
+};
+
+#endif