dune/tectonic/geocoordinate.hh

+#ifndef DUNE_TECTONIC_GEOCOORDINATE_HH
+#define DUNE_TECTONIC_GEOCOORDINATE_HH
+
+// tiny helper to make a common piece of code pleasanter to read
+
+template <class Geometry>
+typename Geometry::GlobalCoordinate geoToPoint(Geometry geo) {
+  assert(geo.corners() == 1);
+  return geo.corner(0);
+}
+
+#endif

dune/tectonic/globalnonlinearity.hh → dune/tectonic/globalfriction.hh

-#ifndef DUNE_TECTONIC_GLOBAL_NONLINEARITY_HH
-#define DUNE_TECTONIC_GLOBAL_NONLINEARITY_HH
+#ifndef DUNE_TECTONIC_GLOBALFRICTION_HH
+#define DUNE_TECTONIC_GLOBALFRICTION_HH
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
@@ -9,9 +9,9 @@
 
 #include <dune/solvers/common/interval.hh>
 
-#include "localfriction.hh"
+#include <dune/tectonic/localfriction.hh>
 
-template <class Matrix, class Vector> class GlobalNonlinearity {
+template <class Matrix, class Vector> class GlobalFriction {
 protected:
   using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
@@ -22,13 +22,12 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
   using LocalMatrix = typename Matrix::block_type;
   using LocalVectorType = typename Vector::block_type;
   size_t static const block_size = LocalVectorType::dimension;
-  using Friction = LocalFriction<block_size>;
+  using LocalNonlinearity = LocalFriction<block_size>;
 
   double operator()(Vector const &x) const {
     double tmp = 0;
     for (size_t i = 0; i < x.size(); ++i) {
-      auto const res = restriction(i);
-      tmp += (*res)(x[i]);
+      tmp += restriction(i)(x[i]);
     }
     return tmp;
   }
@@ -36,14 +35,11 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  std::shared_ptr<LocalFriction<block_size>> virtual restriction(size_t i)
-      const = 0;
+  LocalNonlinearity const virtual &restriction(size_t i) const = 0;
 
   void addHessian(Vector const &v, Matrix &hessian) const {
-    for (size_t i = 0; i < v.size(); ++i) {
-      auto const res = restriction(i);
-      res->addHessian(v[i], hessian[i][i]);
-    }
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addHessian(v[i], hessian[i][i]);
   }
 
   void directionalDomain(Vector const &, Vector const &,
@@ -52,14 +48,13 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
     dom[1] = std::numeric_limits<double>::max();
   }
 
-  void directionalSubDiff(Vector const &u, Vector const &v,
-                          Dune::Solvers::Interval<double> &subdifferential)
-      const {
+  void directionalSubDiff(
+      Vector const &u, Vector const &v,
+      Dune::Solvers::Interval<double> &subdifferential) const {
     subdifferential[0] = subdifferential[1] = 0;
     for (size_t i = 0; i < u.size(); ++i) {
       Dune::Solvers::Interval<double> D;
-      auto const res = restriction(i);
-      res->directionalSubDiff(u[i], v[i], D);
+      restriction(i).directionalSubDiff(u[i], v[i], D);
       subdifferential[0] += D[0];
       subdifferential[1] += D[1];
     }
@@ -71,17 +66,21 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
   }
 
   void addGradient(Vector const &v, Vector &gradient) const {
-    for (size_t i = 0; i < v.size(); ++i) {
-      auto const res = restriction(i);
-      res->addGradient(v[i], gradient[i]);
-    }
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addGradient(v[i], gradient[i]);
   }
 
   double regularity(size_t i, typename Vector::block_type const &x) const {
-    auto const res = restriction(i);
-    return res->regularity(x);
+    return restriction(i).regularity(x);
+  }
+
+  ScalarVector coefficientOfFriction(Vector const &x) const {
+    ScalarVector ret(x.size());
+    for (size_t i = 0; i < x.size(); ++i)
+      ret[i] = restriction(i).coefficientOfFriction(x[i]);
+    return ret;
   }
 
-  void virtual updateLogState(ScalarVector const &logState) = 0;
+  void virtual updateAlpha(ScalarVector const &alpha) = 0;
 };
 #endif

dune/tectonic/globalfrictiondata.hh

-#ifndef GLOBAL_FRICTION_DATA_HH
-#define GLOBAL_FRICTION_DATA_HH
+#ifndef DUNE_TECTONIC_GLOBALFRICTIONDATA_HH
+#define DUNE_TECTONIC_GLOBALFRICTIONDATA_HH
 
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
@@ -23,8 +23,9 @@ template <int dimension> class GlobalFrictionData {
   }
 
 protected:
-  using VirtualFunction = Dune::VirtualFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+  using VirtualFunction =
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>>;
 
   double virtual const &C() const = 0;
   double virtual const &L() const = 0;
@@ -33,5 +34,4 @@ template <int dimension> class GlobalFrictionData {
   VirtualFunction virtual const &b() const = 0;
   double virtual const &mu0() const = 0;
 };
-
 #endif

dune/tectonic/globalruinanonlinearity.hh → dune/tectonic/globalratestatefriction.hh

-#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
-#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
+#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
 
 #include <vector>
 
@@ -10,67 +10,62 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
+#include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/index-in-sorted-range.hh>
 
-#include "globalnonlinearity.hh"
-#include "frictionpotential.hh"
-
-template <class Matrix, class Vector, class GridView>
-class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
+template <class Matrix, class Vector, class ScalarFriction, class GridView>
+class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
 public:
-  using GlobalNonlinearity<Matrix, Vector>::block_size;
-  using typename GlobalNonlinearity<Matrix, Vector>::Friction;
+  using GlobalFriction<Matrix, Vector>::block_size;
+  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
 
 private:
-  using typename GlobalNonlinearity<Matrix, Vector>::ScalarVector;
+  using typename GlobalFriction<Matrix, Vector>::ScalarVector;
 
 public:
-  GlobalRuinaNonlinearity(BoundaryPatch<GridView> const &frictionalBoundary,
-                          GridView const &gridView,
+  GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
                           GlobalFrictionData<block_size> const &frictionInfo,
-                          // Note: passing the following two makes no sense
                           ScalarVector const &weights,
-                          ScalarVector const &normalStress)
-      : restrictions(normalStress.size()) {
-    auto trivialNonlinearity =
-        std::make_shared<Friction>(std::make_shared<TrivialFunction>());
-
+                          ScalarVector const &weightedNormalStress)
+      : restrictions(), localToGlobal(), zeroFriction() {
+    auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
         GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
-      auto const i = vertexMapper.map(*it);
-      auto const coordinate = it->geometry().corner(0);
-      if (not frictionalBoundary.containsVertex(i)) {
-        restrictions[i] = trivialNonlinearity;
+      auto const i = vertexMapper.index(*it);
+
+      if (not frictionalBoundary.containsVertex(i))
         continue;
-      }
-      auto const fp = std::make_shared<FrictionPotential>(
-          weights[i], normalStress[i], frictionInfo(coordinate));
-      restrictions[i] = std::make_shared<Friction>(fp);
-    }
-  }
 
-  void updateLogState(ScalarVector const &logState) override {
-    for (size_t i = 0; i < restrictions.size(); ++i)
-      restrictions[i]->updateLogState(logState[i]);
+      localToGlobal.emplace_back(i);
+      restrictions.emplace_back(weights[i], weightedNormalStress[i],
+                                frictionInfo(geoToPoint(it->geometry())));
+    }
+    assert(restrictions.size() == frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == frictionalBoundary.numVertices());
   }
 
-  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
-    assert(x.size() == restrictions.size());
-    coefficient.resize(restrictions.size());
-    for (size_t i = 0; i < restrictions.size(); ++i)
-      coefficient[i] = restrictions[i]->coefficientOfFriction(x[i]);
+  void updateAlpha(ScalarVector const &alpha) override {
+    for (size_t j = 0; j < restrictions.size(); ++j)
+      restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
   }
 
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  std::shared_ptr<Friction> restriction(size_t i) const override {
-    return restrictions[i];
+  LocalNonlinearity const &restriction(size_t i) const override {
+    auto const index = indexInSortedRange(localToGlobal, i);
+    if (index == localToGlobal.size())
+      return zeroFriction;
+    return restrictions[index];
   }
 
 private:
-  std::vector<std::shared_ptr<Friction>> restrictions;
+  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
+  std::vector<size_t> localToGlobal;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
 };
 #endif

dune/tectonic/gravity.hh

-#ifndef GRAVITY_HH
-#define GRAVITY_HH
+#ifndef DUNE_TECTONIC_GRAVITY_HH
+#define DUNE_TECTONIC_GRAVITY_HH
 
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>

dune/tectonic/index-in-sorted-range.hh

+#ifndef DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
+#define DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
+
+#include <algorithm>
+
+// returns v.size() if value does not exist
+template <typename T>
+size_t indexInSortedRange(std::vector<T> const &v, T value) {
+  size_t const specialReturnValue = v.size();
+
+  auto const b = std::begin(v);
+  auto const e = std::end(v);
+  auto const lb = std::lower_bound(b, e, value);
+
+  if (lb == e) // all elements are strictly smaller
+    return specialReturnValue;
+
+  if (value < *lb) // value falls between to elements
+    return specialReturnValue;
+
+  return std::distance(b, lb);
+}
+#endif

dune/tectonic/localfriction.hh

-#ifndef DUNE_TECTONIC_LOCAL_FRICTION_HH
-#define DUNE_TECTONIC_LOCAL_FRICTION_HH
+#ifndef DUNE_TECTONIC_LOCALFRICTION_HH
+#define DUNE_TECTONIC_LOCALFRICTION_HH
 
 #include <cmath>
 #include <limits>
@@ -10,43 +10,63 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-#include "frictionpotential.hh"
+#include <dune/tectonic/frictionpotential.hh>
 
 template <size_t dimension> class LocalFriction {
 public:
+  virtual ~LocalFriction() {}
+
   using VectorType = Dune::FieldVector<double, dimension>;
   using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
-  explicit LocalFriction(std::shared_ptr<FrictionPotentialWrapper> func)
-      : func(func) {}
+  void virtual updateAlpha(double alpha) = 0;
+  double virtual regularity(VectorType const &x) const = 0;
+  double virtual coefficientOfFriction(VectorType const &x) const = 0;
+  void virtual directionalSubDiff(VectorType const &x, VectorType const &v,
+                                  Dune::Solvers::Interval<double> &D) const = 0;
 
-  double operator()(VectorType const &x) const {
-    return func->evaluate(x.two_norm());
-  }
+  void virtual addHessian(VectorType const &x, MatrixType &A) const = 0;
+
+  void virtual addGradient(VectorType const &x, VectorType &y) const = 0;
+
+  void virtual directionalDomain(
+      VectorType const &, VectorType const &,
+      Dune::Solvers::Interval<double> &dom) const = 0;
+};
+
+template <size_t dimension, class ScalarFriction>
+class WrappedScalarFriction : public LocalFriction<dimension> {
+  using VectorType = typename LocalFriction<dimension>::VectorType;
+  using MatrixType = typename LocalFriction<dimension>::MatrixType;
+
+public:
+  template <typename... Args>
+  WrappedScalarFriction(Args... args)
+      : func_(args...) {}
 
-  void updateLogState(double logState) { func->updateLogState(logState); }
+  void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
 
-  double regularity(VectorType const &x) const {
+  double regularity(VectorType const &x) const override {
     double const xnorm = x.two_norm();
     if (xnorm <= 0.0)
       return std::numeric_limits<double>::infinity();
 
-    return func->regularity(xnorm);
+    return func_.regularity(xnorm);
   }
 
-  double coefficientOfFriction(VectorType const &x) const {
-    return func->coefficientOfFriction(x.two_norm());
+  double coefficientOfFriction(VectorType const &x) const override {
+    return func_.coefficientOfFriction(x.two_norm());
   }
 
   // directional subdifferential: at u on the line u + t*v
   // u and v are assumed to be non-zero
   void directionalSubDiff(VectorType const &x, VectorType const &v,
-                          Dune::Solvers::Interval<double> &D) const {
+                          Dune::Solvers::Interval<double> &D) const override {
     double const xnorm = x.two_norm();
     if (xnorm <= 0.0)
-      D[0] = D[1] = func->differential(0.0) * v.two_norm();
+      D[0] = D[1] = func_.differential(0.0) * v.two_norm();
     else
-      D[0] = D[1] = func->differential(xnorm) * (x * v) / xnorm;
+      D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
   }
 
   /** Formula for the derivative:
@@ -65,14 +85,14 @@ template <size_t dimension> class LocalFriction {
       + \frac {H'(|z|)}{|z|} \operatorname{id}
       \f}
   */
-  void addHessian(VectorType const &x, MatrixType &A) const {
+  void addHessian(VectorType const &x, MatrixType &A) const override {
     double const xnorm2 = x.two_norm2();
     double const xnorm = std::sqrt(xnorm2);
     if (xnorm2 <= 0.0)
       return;
 
-    double const H1 = func->differential(xnorm);
-    double const H2 = func->second_deriv(xnorm);
+    double const H1 = func_.differential(xnorm);
+    double const H2 = func_.second_deriv(xnorm);
 
     double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
     double const idweight = H1 / xnorm;
@@ -90,22 +110,22 @@ template <size_t dimension> class LocalFriction {
     }
   }
 
-  void addGradient(VectorType const &x, VectorType &y) const {
+  void addGradient(VectorType const &x, VectorType &y) const override {
     double const xnorm = x.two_norm();
-    if (std::isinf(func->regularity(xnorm)))
+    if (std::isinf(func_.regularity(xnorm)))
       return;
 
     if (xnorm > 0.0)
-      Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x);
+      Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x);
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
-                         Dune::Solvers::Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const override {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
 
 private:
-  std::shared_ptr<FrictionPotentialWrapper> const func;
+  ScalarFriction func_;
 };
 #endif

dune/tectonic/minimisation.hh

-#ifndef MINIMISATION_HH
-#define MINIMISATION_HH
+#ifndef DUNE_TECTONIC_MINIMISATION_HH
+#define DUNE_TECTONIC_MINIMISATION_HH
 
 #include <dune/common/fmatrix.hh>
 #include <dune/common/fvector.hh>
@@ -9,59 +9,33 @@
 #include <dune/fufem/interval.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
-#include "mydirectionalconvexfunction.hh"
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
 
+// Warning: this exploits the property v*x = 0
 template <class Functional>
 double lineSearch(Functional const &J,
                   typename Functional::LocalVector const &x,
                   typename Functional::LocalVector const &v,
                   Bisection const &bisection) {
   MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
-      computeDirectionalA(J.A, v), computeDirectionalb(J.A, J.b, x, v), *J.phi,
-      x, v);
-
+      J.alpha * v.two_norm2(), J.b * v, J.phi, x, v);
   int count;
   return bisection.minimize(JRest, 0.0, 0.0, count);
 }
 
+/** Minimise a quadratic problem, for which both the quadratic and the
+    nonlinear term have gradients which point in the direction of
+    their arguments */
 template <class Functional>
 void minimise(Functional const &J, typename Functional::LocalVector &x,
-              size_t steps, Bisection const &bisection) {
-  using LocalVector = typename Functional::LocalVector;
-  auto const diff = [](LocalVector const &a, LocalVector const &b) {
-    LocalVector tmp = a;
-    tmp -= b;
-    return tmp.two_norm();
-  };
-
-  LocalVector x_initial = x;
-  LocalVector x_o = x;
-
-  for (size_t step = 0; step < steps; ++step) {
-    LocalVector v;
-    J.descentDirection(x, v);
-
-    double const vnorm = v.two_norm();
-    if (vnorm <= 0.0)
-      return;
-    v /= vnorm;
-
-    double const alpha = lineSearch(J, x, v, bisection);
-    Arithmetic::addProduct(x, alpha, v);
-
-    if (alpha < 1e-14) // TODO
-      break;
-
-    double const correction = diff(x, x_o);
-    double const overallCorrection = diff(x, x_initial);
-    if (overallCorrection <= 0.0)
-      return;
-
-    double const correctionQuotient = correction / overallCorrection;
-    if (correctionQuotient < 0.1) // enough descent; TODO
-      break;
-
-    x_o = x;
-  }
+              Bisection const &bisection) {
+  auto v = J.b;
+  double const vnorm = v.two_norm();
+  if (vnorm <= 0.0)
+    return;
+  v /= vnorm;
+
+  double const alpha = lineSearch(J, x, v, bisection);
+  Arithmetic::addProduct(x, alpha, v);
 }
 #endif

dune/tectonic/myblockproblem.hh

-// Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
+#ifndef DUNE_TECTONIC_MYBLOCKPROBLEM_HH
+#define DUNE_TECTONIC_MYBLOCKPROBLEM_HH
 
-#ifndef MY_BLOCK_PROBLEM_HH
-#define MY_BLOCK_PROBLEM_HH
+// Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
 
 #include <dune/common/bitsetvector.hh>
-#include <dune/common/nullptr.hh>
 #include <dune/common/parametertree.hh>
+#include <dune/common/fmatrixev.hh>
 
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
@@ -13,24 +13,24 @@
 #include <dune/tnnmg/problem-classes/bisection.hh>
 #include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
 
-#include "ellipticenergy.hh"
-#include "globalnonlinearity.hh"
-#include "minimisation.hh"
-#include "mydirectionalconvexfunction.hh"
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/minimisation.hh>
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
+#include <dune/tectonic/quadraticenergy.hh>
 
 /** \brief Base class for problems where each block can be solved with a
  * modified gradient method */
 template <class ConvexProblem>
-class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
+class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
 private:
-  typedef BlockNonlinearGSProblem<ConvexProblem> BNGSP;
+  typedef BlockNonlinearGSProblem<ConvexProblem> Base;
 
 public:
-  using typename BNGSP::ConvexProblemType;
-  using typename BNGSP::LocalMatrixType;
-  using typename BNGSP::LocalVectorType;
-  using typename BNGSP::MatrixType;
-  using typename BNGSP::VectorType;
+  using typename Base::ConvexProblemType;
+  using typename Base::LocalMatrixType;
+  using typename Base::LocalVectorType;
+  using typename Base::MatrixType;
+  using typename Base::VectorType;
 
   size_t static const block_size = ConvexProblem::block_size;
   size_t static const coarse_block_size = block_size;
@@ -55,9 +55,16 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
   };
 
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
-      : BNGSP(parset, problem),
-        parset_(parset),
-        localBisection(0.0, 1.0, 1e-12, false) {}
+      : Base(parset, problem),
+        maxEigenvalues_(problem.f.size()),
+        localBisection(0.0, 1.0, 0.0, true, 0.0) {
+    for (size_t i = 0; i < problem.f.size(); ++i) {
+      LocalVectorType eigenvalues;
+      Dune::FMatrixHelp::eigenValues(problem.A[i][i], eigenvalues);
+      maxEigenvalues_[i] =
+          *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+    }
+  }
 
   std::string getOutput(bool header = false) const {
     if (header) {
@@ -96,15 +103,10 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
     v /= vnorm; // Rescale for numerical stability
 
-    MyDirectionalConvexFunction<
-        GlobalNonlinearity<MatrixType, VectorType>> const
-    psi(computeDirectionalA(problem_.A, v),
-        computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u, v);
-
+    auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
     Dune::Solvers::Interval<double> D;
     psi.subDiff(0, D);
-    // NOTE: Numerical instability can actually get us here
-    if (D[1] > 0)
+    if (D[1] > 0) // NOTE: Numerical instability can actually get us here
       return 0;
 
     int bisectionsteps = 0;
@@ -142,8 +144,8 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
     // compute quadratic part of hessian (linearization.A += problem_.A)
     for (size_t i = 0; i < problem_.A.N(); ++i) {
-      auto const end = problem_.A[i].end();
-      for (auto it = problem_.A[i].begin(); it != end; ++it)
+      auto const end = std::end(problem_.A[i]);
+      for (auto it = std::begin(problem_.A[i]); it != end; ++it)
         linearization.A[i][it.index()] += *it;
     }
 
@@ -163,39 +165,37 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
     // apply truncation to stiffness matrix and rhs
     for (size_t row = 0; row < linearization.A.N(); ++row) {
-      auto const col_end = linearization.A[row].end();
-      for (auto col_it = linearization.A[row].begin(); col_it != col_end;
+      auto const col_end = std::end(linearization.A[row]);
+      for (auto col_it = std::begin(linearization.A[row]); col_it != col_end;
            ++col_it) {
         size_t const col = col_it.index();
         for (size_t i = 0; i < col_it->N(); ++i) {
-          auto const blockEnd = (*col_it)[i].end();
-          for (auto blockIt = (*col_it)[i].begin(); blockIt != blockEnd;
+          auto const blockEnd = std::end((*col_it)[i]);
+          for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd;
                ++blockIt)
-            if (linearization.truncation[row][i] ||
+            if (linearization.truncation[row][i] or
                 linearization.truncation[col][blockIt.index()])
               *blockIt = 0.0;
         }
       }
-
       for (size_t j = 0; j < block_size; ++j)
         if (linearization.truncation[row][j])
           linearization.b[row][j] = 0.0;
     }
-
     for (size_t j = 0; j < block_size; ++j)
       outStream << std::setw(9) << linearization.truncation.countmasked(j);
   }
 
   /** \brief Constructs and returns an iterate object */
   IterateObject getIterateObject() {
-    return IterateObject(parset_, localBisection, problem_);
+    return IterateObject(localBisection, problem_, maxEigenvalues_);
   }
 
 private:
-  Dune::ParameterTree const &parset_;
+  std::vector<double> maxEigenvalues_;
 
   // problem data
-  using BNGSP::problem_;
+  using Base::problem_;
 
   Bisection const localBisection;
 
@@ -212,11 +212,11 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
    * \param bisection The class used to do a scalar bisection
    * \param problem The problem including quadratic part and nonlinear part
    */
-  IterateObject(Dune::ParameterTree const &parset, Bisection const &bisection,
-                ConvexProblem const &problem)
+  IterateObject(Bisection const &bisection, ConvexProblem const &problem,
+                std::vector<double> const &maxEigenvalues)
       : problem(problem),
-        bisection_(bisection),
-        localsteps(parset.get<size_t>("localsolver.steps")) {}
+        maxEigenvalues_(maxEigenvalues),
+        bisection_(bisection) {}
 
 public:
   /** \brief Set the current iterate */
@@ -240,36 +240,34 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
       LocalVectorType &ui, size_t m,
       typename Dune::BitSetVector<block_size>::const_reference ignore) {
     {
-      LocalMatrixType const *localA = nullptr;
-      LocalVectorType localb(problem.f[m]);
-
-      auto const end = problem.A[m].end();
-      for (auto it = problem.A[m].begin(); it != end; ++it) {
+      LocalVectorType localb = problem.f[m];
+      auto const end = std::end(problem.A[m]);
+      for (auto it = std::begin(problem.A[m]); it != end; ++it) {
         size_t const j = it.index();
-        if (j == m)
-          localA = &(*it); // localA = A[m][m]
-        else
-          Arithmetic::subtractProduct(localb, *it, u[j]);
+        Arithmetic::subtractProduct(localb, *it, u[j]); // also the diagonal!
       }
-      assert(localA != nullptr);
+      Arithmetic::addProduct(localb, maxEigenvalues_[m], u[m]);
 
-      auto const phi = problem.phi.restriction(m);
-      EllipticEnergy<block_size> localJ(*localA, localb, phi, ignore);
-      minimise(localJ, ui, localsteps, bisection_);
+      // We minimise over an affine subspace
+      for (size_t j = 0; j < block_size; ++j)
+        if (ignore[j])
+          localb[j] = 0;
+        else
+          ui[j] = 0;
+
+      QuadraticEnergy<
+          typename ConvexProblem::NonlinearityType::LocalNonlinearity>
+          localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
+      minimise(localJ, ui, bisection_);
     }
   }
 
 private:
-  // problem data
   ConvexProblem const &problem;
-
+  std::vector<double> maxEigenvalues_;
   Bisection const bisection_;
-
   // state data for smoothing procedure used by:
   // setIterate, updateIterate, solveLocalProblem
   VectorType u;
-
-  size_t const localsteps;
 };
-
 #endif

dune/tectonic/mydirectionalconvexfunction.hh

+#ifndef DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
+#define DUNE_TECTONIC_MYDIRECTIONALCONVEXFUNCTION_HH
+
 // Copied from dune/tnnmg/problem-classes/directionalconvexfunction.hh
 // Allows phi to be const
 
-#ifndef MY_DIRECTIONAL_CONVEX_FUNCTION_HH
-#define MY_DIRECTIONAL_CONVEX_FUNCTION_HH
-
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-/*
-  1/2 <A(u + hv),u + hv> - <b, u + hv>
-  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
-  localA = <Av,v>
-  localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector>
-double computeDirectionalA(Matrix const &A, Vector const &v) {
-  return Arithmetic::Axy(A, v, v);
-}
-
-template <class Matrix, class Vector>
-double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
-                           Vector const &v) {
-  return Arithmetic::bmAxy(A, b, u, v);
-}
-
 template <class Nonlinearity> class MyDirectionalConvexFunction {
 public:
   using Vector = typename Nonlinearity::VectorType;
@@ -45,8 +26,9 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
     Vector uxv = u;
     Arithmetic::addProduct(uxv, x, v);
     phi.directionalSubDiff(uxv, v, D);
-    D[0] += A * x - b;
-    D[1] += A * x - b;
+    auto const Axmb = A * x - b;
+    D[0] += Axmb;
+    D[1] += Axmb;
   }
 
   void domain(Dune::Solvers::Interval<double> &domain) const {
@@ -65,4 +47,22 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
   Dune::Solvers::Interval<double> dom;
 };
 
+/*
+  1/2 <A(u + hv),u + hv> - <b, u + hv>
+  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
+
+  localA = <Av,v>
+  localb = <b - Au, v>
+*/
+
+template <class Matrix, class Vector, class Nonlinearity>
+MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
+                                                   Vector const &b,
+                                                   Vector const &u,
+                                                   Vector const &v,
+                                                   Nonlinearity const &phi) {
+  return MyDirectionalConvexFunction<Nonlinearity>(
+      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
+}
+
 #endif

dune/tectonic/quadraticenergy.hh

+#ifndef DUNE_TECTONIC_QUADRATICENERGY_HH
+#define DUNE_TECTONIC_QUADRATICENERGY_HH
+
+#include <memory>
+
+template <class NonlinearityTEMPLATE> class QuadraticEnergy {
+public:
+  using Nonlinearity = NonlinearityTEMPLATE;
+  using LocalVector = typename Nonlinearity::VectorType;
+
+  QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
+      : alpha(alpha), b(b), phi(phi) {}
+
+  double const alpha;
+  LocalVector const &b;
+  Nonlinearity const &phi;
+};
+#endif

dune/tectonic/tectonic.hh

-#ifndef DUNE_tectonic.hh
-#define DUNE_tectonic .hh
+#ifndef DUNE_TECTONIC_TECTONIC_HH
+#define DUNE_TECTONIC_TECTONIC_HH
 
 // add your classes here
-
-#endif // DUNE_tectonic.hh
+#endif

m4/Makefile.am

-M4FILES = dune-tectonic.m4
-
-aclocaldir = $(datadir)/aclocal
-aclocal_DATA = $(M4FILES)
-
-EXTRA_DIST = $(M4FILES)
-
-include $(top_srcdir)/am/global-rules

m4/dune-tectonic.m4

-dnl -*- autoconf -*-
-# Macros needed to find dune-tectonic and dependent libraries.  They are called by
-# the macros in ${top_src_dir}/dependencies.m4, which is generated by
-# "dunecontrol autogen"
-
-# Additional checks needed to build dune-tectonic
-# This macro should be invoked by every module which depends on dune-tectonic, as
-# well as by dune-tectonic itself
-AC_DEFUN([DUNE_TECTONIC_CHECKS])
-
-# Additional checks needed to find dune-tectonic
-# This macro should be invoked by every module which depends on dune-tectonic, but
-# not by dune-tectonic itself
-AC_DEFUN([DUNE_TECTONIC_CHECK_MODULE],
-[
-  DUNE_CHECK_MODULES([dune-tectonic],[tectonic/tectonic.hh])
-])

src/CMakeLists.txt

+set(SW_SOURCE_FILES
+  assemblers.cc
+  enumparser.cc
+  hdf5/frictionalboundary-writer.cc
+  hdf5/iteration-writer.cc
+  hdf5/patchinfo-writer.cc
+  hdf5/restart-io.cc
+  hdf5/surface-writer.cc
+  hdf5/time-writer.cc
+  one-body-problem-data/mygeometry.cc
+  one-body-problem-data/mygrid.cc
+  one-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  vtk.cc
+)
+set(UGW_SOURCE_FILES
+  assemblers.cc # FIXME
+  one-body-problem-data/mygrid.cc
+  uniform-grid-writer.cc
+  vtk.cc
+)
+
+foreach(_dim 2 3)
+  set(_sw_target one-body-problem-${_dim}D)
+  set(_ugw_target uniform-grid-writer-${_dim}D)
+
+  add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_ugw_target})
+
+  add_dune_hdf5_flags(${_sw_target})
+
+  set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()

src/Makefile.am

-bin_PROGRAMS = sliding-block-2D
-
-common_sources = \
-	assemblers.cc \
-	boundary_writer.cc \
-	friction_writer.cc \
-	solverfactory.cc \
-	timestepping.cc \
-	vtk.cc
-
-sliding_block_2D_SOURCES = $(common_sources) sliding-block.cc
-sliding_block_2D_CPPFLAGS = \
-	$(AM_CPPFLAGS) \
-	-Ddatadir=\"$(abs_srcdir)/sliding-block-data/\" -DDIM=2
-
-# Some are for clang, others are for gcc
-AM_CXXFLAGS = \
-	-Wall \
-	-Wextra \
-	-Wno-unused-parameter \
-	-Wno-overloaded-virtual \
-	-Wno-new-returns-null \
-	-Wno-unknown-warning-option \
-	-Wno-unknown-pragmas
-
-AM_CPPFLAGS = \
-	-DDUNE_FMatrix_WITH_CHECKING \
-	$(DUNE_CPPFLAGS) \
-	$(PYTHON_CPPFLAGS) \
-	$(ALUGRID_CPPFLAGS) \
-	-I$(top_srcdir)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last).
-LDADD = \
-	$(DUNE_LDFLAGS) $(DUNE_LIBS) \
-	$(ALUGRID_LIBS) \
-	$(PYTHON_LIBS)
-AM_LDFLAGS = \
-	$(DUNE_LDFLAGS) \
-	$(ALUGRID_LDFLAGS) \
-	$(PYTHON_LDFLAGS)
-
-include $(top_srcdir)/am/global-rules

src/assemblers.cc

@@ -7,16 +7,20 @@
 #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/computestress.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/functiontools/p0p1interpolation.hh>
 #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
+#include <dune/tectonic/frictionpotential.hh>
+#include <dune/tectonic/globalratestatefriction.hh>
+
 #include "assemblers.hh"
 
 template <class GridView, int dimension>
@@ -30,10 +34,10 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
     BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarMatrix &frictionalBoundaryMass) {
+    ScalarMatrix &frictionalBoundaryMass) const {
   BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                         LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
-  frictionalBoundaryMassAssembler(frictionalBoundary);
+      frictionalBoundaryMassAssembler(frictionalBoundary);
   vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
                                    frictionalBoundaryMass);
 }
@@ -42,22 +46,22 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleMass(
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
         densityFunction,
-    Matrix &M) {
+    Matrix &M) const {
   // NOTE: We treat the weight as a constant function
   QuadratureRuleKey quadKey(dimension, 0);
 
   WeightedMassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
                         Dune::VirtualFunction<LocalVector, LocalScalarVector>,
                         Dune::ScaledIdentityMatrix<double, dimension>>
-  localWeightedMass(gridView.grid(), densityFunction, quadKey);
+      localWeightedMass(gridView.grid(), densityFunction, quadKey);
   vertexAssembler.assembleOperator(localWeightedMass, M);
 }
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
-                                                          Matrix &A) {
+                                                          Matrix &A) const {
   StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
-  localStiffness(E, nu);
+      localStiffness(E, nu);
   vertexAssembler.assembleOperator(localStiffness, A);
 }
 
@@ -65,82 +69,108 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleViscosity(
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
-    Matrix &C) {
+    Matrix &C) const {
   // NOTE: We treat the weights as constant functions
   QuadratureRuleKey shearViscosityKey(dimension, 0);
   QuadratureRuleKey bulkViscosityKey(dimension, 0);
   VariableCoefficientViscosityAssembler<
       Grid, LocalVertexBasis, LocalVertexBasis,
       Dune::VirtualFunction<LocalVector, LocalScalarVector>> const
-  localViscosity(gridView.grid(), shearViscosity, bulkViscosity,
-                 shearViscosityKey, bulkViscosityKey);
+      localViscosity(gridView.grid(), shearViscosity, bulkViscosity,
+                     shearViscosityKey, bulkViscosityKey);
   vertexAssembler.assembleOperator(localViscosity, C);
 }
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleBodyForce(
     Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-    Vector &f) {
+    Vector &f) const {
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
-  gravityFunctionalAssembler(gravityField);
+      gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
 }
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleNeumann(
     BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
-    Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
+    Dune::VirtualFunction<double, double> const &neumann,
+    double relativeTime) const {
   LocalVector localNeumann(0);
   neumann.evaluate(relativeTime, localNeumann[0]);
-  ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
-      fNeumann);
+      std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
+          localNeumann));
   vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
                                              neumannBoundary);
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleNormalStress(
+void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
     BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarVector &normalStress, double youngModulus, double poissonRatio,
-    Vector const &displacement) {
-
-  Vector traction;
-  Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
-      (frictionalBoundary, displacement, traction, youngModulus, poissonRatio);
-
-  std::vector<typename Vector::block_type> normals;
-  frictionalBoundary.getNormals(normals);
-  for (size_t i = 0; i < traction.size(); ++i) {
-    normalStress[i] = normals[i] * traction[i];
-    assert(normalStress[i] <= 0.0);
+    ScalarVector &weightedNormalStress, double youngModulus,
+    double poissonRatio, Vector const &displacement) const {
+
+  BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
+                                                              displacement);
+  Vector traction(cellBasis.size());
+  NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
+      youngModulus, poissonRatio, &displacementFunction, 1);
+  cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
+                                           frictionalBoundary);
+
+  auto const nodalTractionAverage =
+      interpolateP0ToP1(frictionalBoundary, traction);
+
+  ScalarVector weights;
+  {
+    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+        frictionalBoundaryAssembler(
+            std::make_shared<ConstantFunction<
+                LocalVector, typename ScalarVector::block_type>>(1));
+    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
+                                               weights, frictionalBoundary);
   }
+  auto const normals = frictionalBoundary.getNormals();
+  for (size_t i = 0; i < vertexBasis.size(); ++i)
+    weightedNormalStress[i] =
+        std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
 }
 
 template <class GridView, int dimension>
 auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+    Config::FrictionModel frictionModel,
     BoundaryPatch<GridView> const &frictionalBoundary,
     GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &normalStress)
-    -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>> {
-  // Lump negative normal stress (kludge)
+    ScalarVector const &weightedNormalStress) const
+    -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
+  // Lumping of the nonlinearity
   ScalarVector weights;
   {
-    ConstantFunction<LocalVector, typename ScalarVector::block_type> const
-    constantOneFunction(1);
     NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
-    frictionalBoundaryAssembler(constantOneFunction);
+        frictionalBoundaryAssembler(std::make_shared<
+            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+            1));
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
-  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>(
-      frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+  switch (frictionModel) {
+    case Config::Truncated:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, TruncatedRateState, GridView>>(
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
+    case Config::Regularised:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, RegularisedRateState, GridView>>(
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
+    default:
+      assert(false);
+  }
 }
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleVonMisesStress(
     double youngModulus, double poissonRatio, Vector const &u,
-    ScalarVector &stress) {
+    ScalarVector &stress) const {
   auto const gridDisplacement =
       std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
           vertexBasis, u);

src/assemblers.hh

-#ifndef ASSEMBLERS_HH
-#define ASSEMBLERS_HH
+#ifndef SRC_ASSEMBLERS_HH
+#define SRC_ASSEMBLERS_HH
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/function.hh>
@@ -13,8 +13,10 @@
 #pragma clang diagnostic pop
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
+#include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
-#include <dune/tectonic/globalruinanonlinearity.hh>
+
+#include "enums.hh"
 
 template <class GridView, int dimension> class MyAssembler {
 public:
@@ -27,8 +29,9 @@ template <class GridView, int dimension> class MyAssembler {
   using CellBasis = P0Basis<GridView, double>;
   using VertexBasis = P1NodalBasis<GridView, double>;
 
-  CellBasis cellBasis;
-  VertexBasis vertexBasis;
+  CellBasis const cellBasis;
+  VertexBasis const vertexBasis;
+  GridView const &gridView;
 
 private:
   using Grid = typename GridView::Grid;
@@ -38,7 +41,6 @@ template <class GridView, int dimension> class MyAssembler {
   using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  GridView const &gridView;
   Assembler<CellBasis, CellBasis> cellAssembler;
   Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
@@ -47,42 +49,42 @@ template <class GridView, int dimension> class MyAssembler {
 
   void assembleFrictionalBoundaryMass(
       BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass);
+      ScalarMatrix &frictionalBoundaryMass) const;
 
   void assembleMass(Dune::VirtualFunction<
                         LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M);
+                    Matrix &M) const;
 
-  void assembleElasticity(double E, double nu, Matrix &A);
+  void assembleElasticity(double E, double nu, Matrix &A) const;
 
   void assembleViscosity(
       Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
           shearViscosity,
       Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
           bulkViscosity,
-      Matrix &C);
+      Matrix &C) const;
 
   void assembleBodyForce(
       Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f);
+      Vector &f) const;
 
   void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                        Vector &f,
                        Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime);
+                       double relativeTime) const;
 
-  void assembleNormalStress(BoundaryPatch<GridView> const &frictionalBoundary,
-                            ScalarVector &normalStress, double youngModulus,
-                            double poissonRatio, Vector const &displacement);
+  void assembleWeightedNormalStress(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarVector &weightedNormalStress, double youngModulus,
+      double poissonRatio, Vector const &displacement) const;
 
-  std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector, GridView>>
-  assembleFrictionNonlinearity(
+  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
+      Config::FrictionModel frictionModel,
       BoundaryPatch<GridView> const &frictionalBoundary,
       GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &normalStress);
+      ScalarVector const &weightedNormalStress) const;
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress);
+                              Vector const &u, ScalarVector &stress) const;
 };
-
 #endif

src/assemblers_tmpl.cc

-#ifndef DIM
-#error DIM unset
+#ifndef MY_DIM
+#error MY_DIM unset
 #endif
 
 #include "explicitgrid.hh"
 
-template class MyAssembler<GridView, DIM>;
+template class MyAssembler<GridView, MY_DIM>;

src/boundary_writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "boundary_writer.hh"
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::BoundaryWriter(
-    Vector const &vertexCoordinates,
-    Dune::BitSetVector<1> const &_boundaryNodes, Projector projector)
-    : displacementWriter("displacements", std::fstream::out),
-      velocityWriter("velocities", std::fstream::out),
-      boundaryNodes(_boundaryNodes),
-      projector_(projector) {
-  std::fstream vertexCoordinateWriter("coordinates", std::fstream::out);
-  for (size_t i = 0; i < boundaryNodes.size(); ++i)
-    if (boundaryNodes[i][0])
-      vertexCoordinateWriter << vertexCoordinates[i] << std::endl;
-  vertexCoordinateWriter.close();
-}
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::~BoundaryWriter() {
-  displacementWriter.close();
-  velocityWriter.close();
-}
-
-template <class ScalarVector, class Vector>
-void BoundaryWriter<ScalarVector, Vector>::writeKinetics(Vector const &u,
-                                                         Vector const &v) {
-  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
-    if (!boundaryNodes[i][0])
-      continue;
-
-    displacementWriter << projector_(u[i]) << " ";
-    velocityWriter << projector_(v[i]) << " ";
-  }
-
-  displacementWriter << std::endl;
-  velocityWriter << std::endl;
-}
-
-#include "boundary_writer_tmpl.cc"