Move time stepping of alpha to a class

src/dieterichstateupdater.hh

+#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

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;
 

src/ruinastateupdater.hh

+#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

src/stateupdater.hh

+#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