dune/tectonic/time-stepping/coupledtimestepper.hh

+#ifndef SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
+#define SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
+
+#include <functional>
+#include <memory>
+
+#include <dune/common/parametertree.hh>
+
+#include "../spatial-solving/fixedpointiterator.hh"
+
+template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+class CoupledTimeStepper {
+  using Vector = typename Factory::Vector;
+  using Matrix = typename Factory::Matrix;
+  using IgnoreVector = typename Factory::BitVector;
+  using FixedPointIterator = FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+
+public:
+  using GlobalFriction = typename FixedPointIterator::GlobalFriction;
+  using BitVector = typename FixedPointIterator::BitVector;
+  using ExternalForces = std::vector<std::unique_ptr<const std::function<void(double, Vector &)>>>;
+
+public:
+  CoupledTimeStepper(double finalTime,
+                     Dune::ParameterTree const &parset,
+                     const NBodyAssembler& nBodyAssembler,
+                     const IgnoreVector& ignoreNodes,
+                     GlobalFriction& globalFriction,
+                     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
+                     Updaters updaters,
+                     const ErrorNorms& errorNorms,
+                     ExternalForces& externalForces);
+
+  template <class LinearSolver>
+  FixedPointIterationCounter step(std::shared_ptr<LinearSolver>& linearSolver, double relativeTime, double relativeTau);
+
+private:
+  double finalTime_;
+  Dune::ParameterTree const &parset_;
+  const NBodyAssembler& nBodyAssembler_;
+  const IgnoreVector& ignoreNodes_;
+
+  GlobalFriction& globalFriction_;
+  const std::vector<const BitVector*>& bodywiseNonmortarBoundaries_;
+
+  Updaters updaters_;
+  ExternalForces& externalForces_;
+  const ErrorNorms& errorNorms_;
+};
+
+#endif

dune/tectonic/time-stepping/coupledtimestepper_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+
+#include "../data-structures/network/contactnetwork.hh"
+#include "../data-structures/friction/globalfriction.hh"
+
+#include "rate/rateupdater.hh"
+#include "state/stateupdater.hh"
+#include "updaters.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+
+template class CoupledTimeStepper<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
+template FixedPointIterationCounter CoupledTimeStepper<MySolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::step<LinearSolver>(std::shared_ptr<LinearSolver>&, double, double);
+
+using NoFriction = ZeroNonlinearity;
+using NoFrictionFunctional = Functional<Matrix&, Vector&, NoFriction&, Vector&, Vector&, double>;
+using NoFrictionSolverFactory = SolverFactory<NoFrictionFunctional, BitVector>;
+
+template class CoupledTimeStepper<NoFrictionSolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>;
+template FixedPointIterationCounter CoupledTimeStepper<NoFrictionSolverFactory, MyNBodyAssembler, MyUpdaters, ErrorNorms>::step<LinearSolver>(std::shared_ptr<LinearSolver>&, double, double);

src/time-stepping/rate.cc → dune/tectonic/time-stepping/rate.cc

@@ -6,23 +6,26 @@
 #include "rate/backward_euler.hh"
 #include "rate/newmark.hh"
 
-template <class Vector, class Matrix, class Function, int dimension>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
-initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial) {
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto initRateUpdater(Config::scheme scheme,
+                const BoundaryFunctions& velocityDirichletFunctions,
+                const BoundaryNodes& velocityDirichletNodes,
+                const Matrices<Matrix,2>& matrices,
+                const std::vector<Vector>& u_initial,
+                const std::vector<Vector>& v_initial,
+                const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
   switch (scheme) {
     case Config::Newmark:
-      return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
+      return std::make_shared<Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     case Config::BackwardEuler:
       return std::make_shared<
-          BackwardEuler<Vector, Matrix, Function, dimension>>(
+          BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     default:
       assert(false);
   }

dune/tectonic/time-stepping/rate.hh

+#ifndef SRC_TIME_STEPPING_RATE_HH
+#define SRC_TIME_STEPPING_RATE_HH
+
+#include <memory>
+
+#include "../data-structures/enums.hh"
+#include "rate/rateupdater.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto initRateUpdater(Config::scheme scheme,
+                const BoundaryFunctions& velocityDirichletFunctions,
+                const BoundaryNodes& velocityDirichletNodes,
+                const Matrices<Matrix,2>& matrices,
+                const std::vector<Vector>& u_initial,
+                const std::vector<Vector>& v_initial,
+                const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>;
+#endif

dune/tectonic/time-stepping/rate/CMakeLists.txt

+add_custom_target(tectonic_dune_time-stepping_rate SOURCES
+  backward_euler.hh
+  backward_euler.cc
+  newmark.hh
+  newmark.cc
+  rateupdater.hh
+  rateupdater.cc
+)
+
+#install headers
+install(FILES
+  backward_euler.hh
+  newmark.hh
+  rateupdater.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/time-stepping/rate/backward_euler.cc

+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include "backward_euler.hh"
+
+#include "../../utils/debugutils.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::BackwardEuler(
+        const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+        const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+        const BoundaryNodes& _dirichletNodes,
+        const BoundaryFunctions& _dirichletFunctions) :
+    RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+        _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
+        _dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
+        const std::vector<Vector>& ell,
+        double _tau,
+        double relativeTime,
+        std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+        std::vector<Matrix>& AM) {
+
+    for (size_t i=0; i<this->u.size(); i++) {
+        this->tau = _tau;
+
+        /* We start out with the formulation
+
+             M a + C v + A u = ell
+
+           Backward Euler means
+
+             a1 = 1.0/tau ( v1 - v0 )
+             u1 = tau v1 + u0
+
+           in summary, we get at time t=1
+
+             M [1.0/tau ( v1 - v0 )] + C v1
+             + A [tau v1 + u0] = ell
+
+           or
+
+             1.0/tau M v1 + C v1 + tau A v1
+             = [1.0/tau M + C + tau A] v1
+             = ell + 1.0/tau M v0 - A u0
+        */
+
+        // set up LHS (for fixed tau, we'd only really have to do this once)
+        Matrix& LHS = AM[i];
+        {
+        Dune::MatrixIndexSet indices(
+                    this->matrices.elasticity[i]->N(),
+                    this->matrices.elasticity[i]->M());
+        indices.import(*this->matrices.elasticity[i]);
+        indices.import(*this->matrices.mass[i]);
+        indices.import(*this->matrices.damping[i]);
+        indices.exportIdx(LHS);
+        }
+        LHS = 0.0;
+        Dune::MatrixVector::addProduct(LHS, 1.0 / this->tau, *this->matrices.mass[i]);
+        Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+        Dune::MatrixVector::addProduct(LHS, this->tau, *this->matrices.elasticity[i]);
+
+        // set up RHS
+        {
+        Vector& rhss = rhs[i];
+        rhss = ell[i];
+        Dune::MatrixVector::addProduct(rhss, 1.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+        Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+        }
+
+    }
+
+    iterate = this->v_o;
+
+    for (size_t i=0; i<iterate.size(); i++) {
+        auto& bodyIterate = iterate[i];
+
+        const auto& bodyDirichletNodes = this->dirichletNodes[i];
+        const auto& bodyDirichletFunctions = this->dirichletFunctions[i];
+        for (size_t bc=0; bc<bodyDirichletNodes.size(); ++bc) {
+            const auto& bcDirichletNodes = *bodyDirichletNodes[bc];
+            size_t dim = bcDirichletNodes[0].size();
+
+            double dirichletValue;
+            bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+
+            for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
+                for (size_t j=0; j<dim; ++j) {
+                    if (bcDirichletNodes[k][j]) {
+                        iterate[k][j] = dirichletValue;
+                    }
+                }
+            }
+        }
+    }
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) {
+    // v1 = 1/tau ( u1 - u0 )
+    v1Obstacles.resize(u0.size());
+    v1Obstacles = 0;
+
+    Dune::MatrixVector::addProduct(v1Obstacles, 1.0 / this->tau, uObstacles);
+    Dune::MatrixVector::subtractProduct(v1Obstacles, 1.0 / this->tau, u0);
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::postProcess(const std::vector<Vector>& iterate) {
+    this->postProcessCalled = true;
+
+    this->v = iterate;
+
+    this->u = this->u_o;
+
+    for (size_t i=0; i<this->u.size(); i++) {
+        Dune::MatrixVector::addProduct(this->u[i], this->tau, this->v[i]);
+
+        Vector& ai = this->a[i];
+        ai = this->v[i];
+        ai -= this->v_o[i];
+        ai /= this->tau;
+    }
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::clone() const
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
+    return std::make_shared<BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(*this);
+}

dune/tectonic/time-stepping/rate/backward_euler.hh

+#ifndef SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
+#define SRC_TIME_STEPPING_RATE_BACKWARD_EULER_HH
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class BackwardEuler : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
+public:
+    BackwardEuler(
+            const Matrices<Matrix,2> &_matrices,
+            const std::vector<Vector> &_u_initial,
+            const std::vector<Vector> &_v_initial,
+            const std::vector<Vector> &_a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
+
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) override;
+
+    virtual void velocityObstacles(const Vector&, const Vector&, const Vector&, Vector&) override;
+
+    virtual void postProcess(const std::vector<Vector>&) override;
+
+    virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
+};
+#endif

dune/tectonic/time-stepping/rate/newmark.cc

+#include <dune/matrix-vector/axpy.hh>
+#include <dune/istl/matrixindexset.hh>
+
+#include "newmark.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::Newmark(
+        const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+        const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+        const BoundaryNodes& _dirichletNodes,
+        const BoundaryFunctions& _dirichletFunctions) :
+    RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+        _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
+        _dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::setup(
+        const std::vector<Vector>& ell,
+        double _tau,
+        double relativeTime,
+        std::vector<Vector>& rhs, std::vector<Vector>& iterate,
+        std::vector<Matrix>& AM) {
+
+    const size_t bodyCount = this->u.size();
+
+    rhs.resize(bodyCount);
+    iterate.resize(bodyCount);
+    AM.resize(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+        this->tau = _tau;
+
+        /*  We start out with the formulation
+
+            M a + C v + A u = ell
+
+            Newmark means
+
+            a1 = 2/tau ( v1 - v0 ) - a0
+            u1 = tau/2 ( v1 + v0 ) + u0
+
+            in summary, we get at time t=1
+
+            M [2/tau ( u1 - u0 ) - a0] + C v1
+            + A [tau/2 ( v1 + v0 ) + u0] = ell
+
+            or
+
+            2/tau M v1 + C v1 + tau/2 A v1
+            = [2/tau M + C + tau/2 A] v1
+            = ell + 2/tau M v0 + M a0
+            - tau/2 A v0 - A u0
+        */
+
+        // set up LHS (for fixed tau, we'd only really have to do this once)
+        Matrix& LHS = AM[i];
+        {
+        Dune::MatrixIndexSet indices(
+                    this->matrices.elasticity[i]->N(),
+                    this->matrices.elasticity[i]->M());
+        indices.import(*this->matrices.elasticity[i]);
+        indices.import(*this->matrices.mass[i]);
+        indices.import(*this->matrices.damping[i]);
+        indices.exportIdx(LHS);
+        }
+        LHS = 0.0;
+        Dune::MatrixVector::addProduct(LHS, 2.0 / this->tau, *this->matrices.mass[i]);
+        Dune::MatrixVector::addProduct(LHS, 1.0, *this->matrices.damping[i]);
+        Dune::MatrixVector::addProduct(LHS, this->tau / 2.0, *this->matrices.elasticity[i]);
+
+        // set up RHS
+        {
+        Vector& rhss = rhs[i];
+        rhss = ell[i];
+        Dune::MatrixVector::addProduct(rhss, 2.0 / this->tau, *this->matrices.mass[i],
+                           this->v_o[i]);
+        Dune::MatrixVector::addProduct(rhss, *this->matrices.mass[i], this->a_o[i]);
+        Dune::MatrixVector::subtractProduct(rhss, this->tau / 2.0, *this->matrices.elasticity[i],
+                           this->v_o[i]);
+        Dune::MatrixVector::subtractProduct(rhss, *this->matrices.elasticity[i], this->u_o[i]);
+        }
+    }
+
+    iterate = this->v_o;
+
+    for (size_t i=0; i<bodyCount; i++) {
+        auto& bodyIterate = iterate[i];
+
+        const auto& bodyDirichletNodes = this->dirichletNodes[i];
+        const auto& bodyDirichletFunctions = this->dirichletFunctions[i];
+        for (size_t bc=0; bc<bodyDirichletNodes.size(); ++bc) {
+            const auto& bcDirichletNodes = *bodyDirichletNodes[bc];
+            size_t dim = bcDirichletNodes[0].size();
+
+            double dirichletValue;
+            bodyDirichletFunctions[bc]->evaluate(relativeTime, dirichletValue);
+
+            //std::cout << "dirichletValue: " << dirichletValue << std::endl;
+
+                  
+            for (size_t k=0; k<bcDirichletNodes.size() ; ++k) {
+                for (size_t j=0; j<dim; ++j) {
+                    if (bcDirichletNodes[k][j]) {
+                              //std::cout << k << " " << j << std::endl;
+                        bodyIterate[k][j] = dirichletValue;
+                    }
+                }
+            }
+        }
+    }
+    //print(iterate, "iterate:");
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) {
+    // v1 = 2/tau ( u1 - u0 ) - v0
+    v1Obstacles.resize(v0.size());
+    v1Obstacles = 0;
+
+    Dune::MatrixVector::addProduct(v1Obstacles, 2.0 / this->tau, uObstacles);
+    Dune::MatrixVector::subtractProduct(v1Obstacles, 2.0 / this->tau, u0);
+    v1Obstacles -= v0;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::postProcess(const std::vector<Vector>& iterate) {
+    this->postProcessCalled = true;
+
+    this->v = iterate;
+
+    // u1 = tau/2 ( v1 + v0 ) + u0
+    this->u = this->u_o;
+
+    for (size_t i=0; i<this->u.size(); i++) {
+        Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v[i]);
+        Dune::MatrixVector::addProduct(this->u[i], this->tau / 2.0, this->v_o[i]);
+
+        // a1 = 2/tau ( v1 - v0 ) - a0
+        this->a[i] = 0.0;
+        Dune::MatrixVector::addProduct(this->a[i], 2.0 / this->tau, this->v[i]);
+        Dune::MatrixVector::subtractProduct(this->a[i], 2.0 / this->tau, this->v_o[i]);
+        Dune::MatrixVector::subtractProduct(this->a[i], 1.0, this->a_o[i]);
+    }
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+auto Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::clone() const
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> {
+
+    return std::make_shared<Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>(*this);
+}

dune/tectonic/time-stepping/rate/newmark.hh

+#ifndef SRC_TIME_STEPPING_RATE_NEWMARK_HH
+#define SRC_TIME_STEPPING_RATE_NEWMARK_HH
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class Newmark : public RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes> {
+public:
+    Newmark(
+            const Matrices<Matrix,2>& _matrices,
+            const std::vector<Vector>& _u_initial,
+            const std::vector<Vector>& _v_initial,
+            const std::vector<Vector>& _a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
+
+    virtual void setup(
+            const std::vector<Vector>&,
+            double,
+            double,
+            std::vector<Vector>&,
+            std::vector<Vector>&,
+            std::vector<Matrix>&) override;
+
+    virtual void velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) override;
+
+    virtual void postProcess(const std::vector<Vector>&) override;
+
+    virtual auto clone() const -> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> override;
+};
+#endif

dune/tectonic/time-stepping/rate/rateupdater.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "rateupdater.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::RateUpdater(
+    const Matrices<Matrix,2>& _matrices, const std::vector<Vector>& _u_initial,
+    const std::vector<Vector>& _v_initial, const std::vector<Vector>& _a_initial,
+    const BoundaryNodes& _dirichletNodes,
+    const BoundaryFunctions& _dirichletFunctions)
+    : matrices(_matrices),
+      u(_u_initial),
+      v(_v_initial),
+      a(_a_initial),
+      dirichletNodes(_dirichletNodes),
+      dirichletFunctions(_dirichletFunctions) {}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::nextTimeStep() {
+    u_o = u;
+    v_o = v;
+    a_o = a;
+    postProcessCalled = false;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractDisplacement(std::vector<Vector>& displacements) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+    displacements = u;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractOldDisplacement(std::vector<Vector>& displacements) const {
+    displacements = u_o;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractVelocity(std::vector<Vector>& velocity) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+    velocity = v;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractOldVelocity(std::vector<Vector>& oldVelocity) const {
+    oldVelocity = v_o;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+void RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::extractAcceleration(std::vector<Vector>& acceleration) const {
+    if (!postProcessCalled)
+        DUNE_THROW(Dune::Exception, "It seems you forgot to call postProcess!");
+
+    acceleration = a;
+}
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+const Matrices<Matrix,2>& RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>::getMatrices() const {
+    return matrices;
+}
+
+
+#include "backward_euler.cc"
+#include "newmark.cc"
+#include "rateupdater_tmpl.cc"

dune/tectonic/time-stepping/rate/rateupdater.hh

+#ifndef SRC_TIME_STEPPING_RATE_RATEUPDATER_HH
+#define SRC_TIME_STEPPING_RATE_RATEUPDATER_HH
+
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+
+#include "../../data-structures/matrices.hh"
+
+template <class Vector, class Matrix, class BoundaryFunctions, class BoundaryNodes>
+class RateUpdater {
+protected:
+    RateUpdater(
+            const Matrices<Matrix,2>& _matrices,
+            const std::vector<Vector>& _u_initial,
+            const std::vector<Vector>& _v_initial,
+            const std::vector<Vector>& _a_initial,
+            const BoundaryNodes& _dirichletNodes,
+            const BoundaryFunctions& _dirichletFunctions);
+
+public:
+    void nextTimeStep();
+    void virtual setup(
+            const std::vector<Vector>& ell,
+            double _tau,
+            double relativeTime,
+            std::vector<Vector>& rhs,
+            std::vector<Vector>& iterate,
+            std::vector<Matrix>& AB) = 0;
+
+  void virtual postProcess(const std::vector<Vector>& iterate) = 0;
+  void virtual velocityObstacles(const Vector& u0, const Vector& uObstacles, const Vector& v0, Vector& v1Obstacles) = 0;
+  void extractDisplacement(std::vector<Vector>& displacements) const;
+  void extractOldDisplacement(std::vector<Vector>& displacements) const;
+  void extractVelocity(std::vector<Vector>& velocity) const;
+  void extractOldVelocity(std::vector<Vector>& velocity) const;
+  void extractAcceleration(std::vector<Vector>& acceleration) const;
+
+  const Matrices<Matrix,2>& getMatrices() const;
+
+  std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>> virtual clone() const = 0;
+
+protected:
+  const Matrices<Matrix,2>& matrices;
+  std::vector<Vector> u, v, a;
+  const BoundaryNodes& dirichletNodes;
+  const BoundaryFunctions& dirichletFunctions;
+
+  std::vector<Vector> u_o, v_o, a_o;
+  double tau;
+
+  bool postProcessCalled = true;
+};
+#endif

dune/tectonic/time-stepping/rate/rateupdater_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "../../explicitvectors.hh"
+
+#include "../../data-structures/network/contactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+template class RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class Newmark<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+template class BackwardEuler<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;

dune/tectonic/time-stepping/rate_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+#include "../data-structures/network/contactnetwork.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+
+template
+auto initRateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>(
+    Config::scheme scheme,
+    const BoundaryFunctions&  velocityDirichletFunctions,
+    const BoundaryNodes& velocityDirichletNodes,
+    const Matrices<Matrix, 2>& matrices,
+    const std::vector<Vector>& u_initial,
+    const std::vector<Vector>& v_initial,
+    const std::vector<Vector>& a_initial)
+-> std::shared_ptr<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>>;

dune/tectonic/time-stepping/state.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "state.hh"
+#include "state/ageinglawstateupdater.cc"
+#include "state/sliplawstateupdater.cc"
+
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
+auto initStateUpdater(
+        Config::stateModel model,
+        const std::vector<ScalarVector>& alpha_initial,
+        const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+        const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>> {
+
+  assert(contactCouplings.size() == couplings.size());
+
+  std::vector<size_t> leafVertexCounts(alpha_initial.size(), 0);
+  for (size_t i=0; i<leafVertexCounts.size(); i++) {
+    leafVertexCounts[i] = alpha_initial[i].size();
+  }
+
+  auto stateUpdater = std::make_shared<StateUpdater<ScalarVector, Vector>>(leafVertexCounts);
+
+  switch (model) {
+    case Config::AgeingLaw:
+      for (size_t i=0; i<couplings.size(); i++) {
+        const auto& coupling = couplings[i];
+        const auto nonmortarIdx = coupling->gridIdx_[0];
+
+        auto localUpdater =  std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(
+                    alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
+
+        stateUpdater->addLocalUpdater(localUpdater);
+      }
+      break;
+    case Config::SlipLaw:
+      for (size_t i=0; i<couplings.size(); i++) {
+        const auto& coupling = couplings[i];
+        const auto nonmortarIdx = coupling->gridIdx_[0];
+
+        auto localUpdater =  std::make_shared<SlipLawStateUpdater<ScalarVector, Vector>>(
+                  alpha_initial[nonmortarIdx], *contactCouplings[i]->nonmortarBoundary().getVertices(), coupling->frictionData().L(), coupling->frictionData().V0());
+        localUpdater->setBodyIndex(nonmortarIdx);
+
+        stateUpdater->addLocalUpdater(localUpdater);
+      }
+      break;
+    default:
+      assert(false);
+      break;
+    }
+
+    return stateUpdater;
+}
+
+#include "state_tmpl.cc"

dune/tectonic/time-stepping/state.hh

+#ifndef SRC_TIME_STEPPING_STATE_HH
+#define SRC_TIME_STEPPING_STATE_HH
+
+#include <memory>
+
+#include "../data-structures/enums.hh"
+#include "state/ageinglawstateupdater.hh"
+#include "state/sliplawstateupdater.hh"
+#include "state/stateupdater.hh"
+
+template <class ScalarVector, class Vector, class ContactCoupling, class FrictionCouplingPair>
+auto initStateUpdater(
+        Config::stateModel model,
+        const std::vector<ScalarVector>& alpha_initial,
+        const std::vector<std::shared_ptr<ContactCoupling>>& contactCouplings, // contains nonmortarBoundary
+        const std::vector<std::shared_ptr<FrictionCouplingPair>>& couplings) // contains frictionInfo
+-> std::shared_ptr<StateUpdater<ScalarVector, Vector>>;
+#endif

dune/tectonic/time-stepping/state/CMakeLists.txt

+add_custom_target(tectonic_dune_time-stepping_state SOURCES
+  ageinglawstateupdater.hh
+  ageinglawstateupdater.cc
+  sliplawstateupdater.hh
+  sliplawstateupdater.cc
+  stateupdater.hh
+)
+
+#install headers
+install(FILES
+  ageinglawstateupdater.hh
+  sliplawstateupdater.hh
+  stateupdater.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

src/time-stepping/state/ageinglawstateupdater.cc → dune/tectonic/time-stepping/state/ageinglawstateupdater.cc

 #include <cmath>
 
 #include "ageinglawstateupdater.hh"
-#include "../../tobool.hh"
+#include "../../utils/tobool.hh"
+#include "../../utils/debugutils.hh"
 
 template <class ScalarVector, class Vector>
 AgeingLawStateUpdater<ScalarVector, Vector>::AgeingLawStateUpdater(
-    ScalarVector const &_alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L, double _V0)
-    : alpha(_alpha_initial), nodes(_nodes), L(_L), V0(_V0) {}
+        const ScalarVector& alpha_initial,
+        const BitVector& nodes,
+        const double L,
+        const double V0) :
+    nodes_(nodes),
+    L_(L),
+    V0_(V0) {
+
+    localToGlobal_.resize(nodes_.count());
+    alpha_.resize(localToGlobal_.size());
+
+    size_t localIdx = 0;
+    for (size_t i=0; i<nodes_.size(); i++) {
+        if (not toBool(nodes_[i]))
+            continue;
+
+        localToGlobal_[localIdx] = i;
+        alpha_[localIdx] = alpha_initial[i];
+        localIdx++;
+    }
+}
 
 template <class ScalarVector, class Vector>
 void AgeingLawStateUpdater<ScalarVector, Vector>::nextTimeStep() {
-  alpha_o = alpha;
+  alpha_o_ = alpha_;
 }
 
 template <class ScalarVector, class Vector>
-void AgeingLawStateUpdater<ScalarVector, Vector>::setup(double _tau) {
-  tau = _tau;
+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) {
+auto liftSingularity(
+        double c,
+        double x) {
+
   if (x <= 0)
     return -c;
   else
@@ -31,27 +53,36 @@ double liftSingularity(double c, double 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));
-  }
+void AgeingLawStateUpdater<ScalarVector, Vector>::solve(const Vector& velocity_field) {
+
+    for (size_t i=0; i<alpha_.size(); ++i) {
+        auto tangentVelocity = velocity_field[localToGlobal_[i]];
+        tangentVelocity[0] = 0.0;
+
+        double const V = tangentVelocity.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;
+        ScalarVector& alpha) {
+
+    //std::cout << "alpha size: " << alpha.size() << " nodes_.size() " << nodes_.size() << std::endl;
+    if (alpha.size() != nodes_.size()) {
+        alpha.resize(nodes_.size());
+    }
+
+    for (size_t i=0; i<localToGlobal_.size(); i++) {
+        alpha[localToGlobal_[i]] = alpha_[i];
+    }
 }
 
 template <class ScalarVector, class Vector>
-std::shared_ptr<StateUpdater<ScalarVector, Vector>>
-AgeingLawStateUpdater<ScalarVector, Vector>::clone() const {
+auto AgeingLawStateUpdater<ScalarVector, Vector>::clone() const
+ -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> {
+
   return std::make_shared<AgeingLawStateUpdater<ScalarVector, Vector>>(*this);
 }

dune/tectonic/time-stepping/state/ageinglawstateupdater.hh

+#ifndef SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+#define SRC_TIME_STEPPING_STATE_AGEINGLAWSTATEUPDATER_HH
+
+#include <dune/common/bitsetvector.hh>
+
+#include "stateupdater.hh"
+
+template <class ScalarVector, class Vector>
+class AgeingLawStateUpdater : public LocalStateUpdater<ScalarVector, Vector> {
+private:
+    using BitVector = Dune::BitSetVector<1>;
+
+public:
+    AgeingLawStateUpdater(
+            const ScalarVector& alpha_initial,
+            const BitVector& nodes,
+            const double L,
+            const double V0);
+
+  void nextTimeStep() override;
+  void setup(double tau) override;
+  void solve(const Vector& velocity_field) override;
+  void extractAlpha(ScalarVector&) override;
+
+  auto clone() const -> std::shared_ptr<LocalStateUpdater<ScalarVector, Vector>> override;
+
+private:
+  std::vector<int> localToGlobal_;
+
+  ScalarVector alpha_o_;
+  ScalarVector alpha_;
+  const BitVector& nodes_;
+  const double L_;
+  const double V0_;
+  double tau_;
+};
+#endif

dune/tectonic/time-stepping/state/explicit.aux

+\relax