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 DUNE_TECTONIC_GLOBALFRICTIONDATA_HH
+#define DUNE_TECTONIC_GLOBALFRICTIONDATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/frictiondata.hh>
+
+template <class DomainType>
+double evaluateScalarFunction(
+    Dune::VirtualFunction<DomainType, Dune::FieldVector<double, 1>> const &f,
+    DomainType const &x) {
+  Dune::FieldVector<double, 1> ret;
+  f.evaluate(x, ret);
+  return ret;
+};
+
+template <int dimension> class GlobalFrictionData {
+public:
+  FrictionData operator()(Dune::FieldVector<double, dimension> const &x) const {
+    return FrictionData(C(), L(), V0(), evaluateScalarFunction(a(), x),
+                        evaluateScalarFunction(b(), x), mu0());
+  }
+
+protected:
+  using VirtualFunction =
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>>;
+
+  double virtual const &C() const = 0;
+  double virtual const &L() const = 0;
+  double virtual const &V0() const = 0;
+  VirtualFunction virtual const &a() const = 0;
+  VirtualFunction virtual const &b() const = 0;
+  double virtual const &mu0() const = 0;
+};
+#endif

dune/tectonic/globalratestatefriction.hh

+#ifndef DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+#define DUNE_TECTONIC_GLOBALRATESTATEFRICTION_HH
+
+#include <vector>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+#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>
+
+template <class Matrix, class Vector, class ScalarFriction, class GridView>
+class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
+public:
+  using GlobalFriction<Matrix, Vector>::block_size;
+  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
+
+private:
+  using typename GlobalFriction<Matrix, Vector>::ScalarVector;
+
+public:
+  GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
+                          GlobalFrictionData<block_size> const &frictionInfo,
+                          ScalarVector const &weights,
+                          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.index(*it);
+
+      if (not frictionalBoundary.containsVertex(i))
+        continue;
+
+      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 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.
+  */
+  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<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
+  std::vector<size_t> localToGlobal;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
+};
+#endif

dune/tectonic/globalruinanonlinearity.hh

-#ifndef DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
-#define DUNE_TECTONIC_GLOBAL_RUINA_NONLINEARITY_HH
-
-#include <vector>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include "globalnonlinearity.hh"
-#include "frictionpotential.hh"
-
-template <class Matrix, class Vector>
-class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
-public:
-  using GlobalNonlinearity<Matrix, Vector>::block_size;
-  using typename GlobalNonlinearity<Matrix, Vector>::Friction;
-
-private:
-  using typename GlobalNonlinearity<Matrix, Vector>::ScalarVector;
-
-public:
-  GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
-                          ScalarVector const &weights, FrictionData const &fd)
-      : restrictions(frictionalNodes.size()) {
-    auto trivialNonlinearity =
-        std::make_shared<Friction>(std::make_shared<TrivialFunction>());
-    for (size_t i = 0; i < restrictions.size(); ++i) {
-      if (not frictionalNodes[i][0]) {
-        restrictions[i] = trivialNonlinearity;
-        continue;
-      }
-      auto const fp = std::make_shared<FrictionPotential>(weights[i], fd);
-      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]);
-  }
-
-  /*
-    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];
-  }
-
-private:
-  std::vector<std::shared_ptr<Friction>> restrictions;
-};
-#endif

dune/tectonic/gravity.hh

+#ifndef DUNE_TECTONIC_GRAVITY_HH
+#define DUNE_TECTONIC_GRAVITY_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+template <int dimension>
+class Gravity
+    : public Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                                   Dune::FieldVector<double, dimension>> {
+public:
+  Gravity(
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>> const &_densityField,
+      Dune::FieldVector<double, dimension> const &_zenith, double _gravity)
+      : densityField(_densityField), zenith(_zenith), gravity(_gravity) {}
+
+  void evaluate(Dune::FieldVector<double, dimension> const &x,
+                Dune::FieldVector<double, dimension> &y) const {
+    y = zenith;
+
+    Dune::FieldVector<double, 1> density;
+    densityField.evaluate(x, density);
+
+    y *= -gravity * density;
+  }
+
+private:
+  Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                        Dune::FieldVector<double, 1>> const &densityField;
+  Dune::FieldVector<double, dimension> const &zenith;
+  double const gravity;
+};
+#endif

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,56 +10,63 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-#include "frictionpotential.hh"
-
-// In order to compute (x * y) / |x|, compute x/|x| first
-template <class Vector>
-double dotFirstNormalised(Vector const &x, Vector const &y) {
-  double const xnorm = x.two_norm();
-  if (xnorm <= 0.0)
-    return 0.0;
-
-  size_t const xsize = x.size();
-  assert(xsize == y.size());
-
-  double sum = 0;
-  for (size_t i = 0; i < xsize; ++i)
-    sum += x[i] / xnorm * y[i];
-
-  return sum;
-}
+#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 updateLogState(double logState) { func->updateLogState(logState); }
+  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...) {}
 
-  double regularity(VectorType const &x) const {
+  void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
+
+  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 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) * dotFirstNormalised(x, v);
+      D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
   }
 
   /** Formula for the derivative:
@@ -78,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;
@@ -103,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,37 +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;
-
-  for (size_t step = 0; step < steps; ++step) {
-    LocalVector v;
-    J.gradient(x, v);
-
-    double const vnorm = v.two_norm();
-    if (vnorm <= 0.0)
-      return;
-    v *= -1;
-
-    double const alpha = lineSearch(J, x, v, bisection);
-    Arithmetic::addProduct(x, alpha, v);
-  }
+              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>
 #include <dune/solvers/computeenergy.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineargsproblem.hh>
 
-#include "globalnonlinearity.hh"
-#include "minimisation.hh"
-#include "mydirectionalconvexfunction.hh"
-#include "ellipticenergy.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 {
+template <class ConvexProblem>
+class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
+private:
+  typedef BlockNonlinearGSProblem<ConvexProblem> Base;
+
 public:
-  using ConvexProblemType = ConvexProblem;
-  using VectorType = typename ConvexProblem::VectorType;
-  using MatrixType = typename ConvexProblem::MatrixType;
-  using LocalVector = typename ConvexProblem::LocalVectorType;
-  using LocalMatrix = typename ConvexProblem::LocalMatrixType;
+  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;
 
-  /** \brief Solves one local system using a modified gradient method */
+  /** \brief Solves one local system */
   class IterateObject;
 
   struct Linearization {
     size_t static const block_size = coarse_block_size;
 
-    using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrix;
+    using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrixType;
     using MatrixType = Dune::BCRSMatrix<typename Linearization::LocalMatrix>;
     using VectorType =
         Dune::BlockVector<Dune::FieldVector<double, Linearization::block_size>>;
@@ -48,10 +54,17 @@ template <class ConvexProblem> class MyBlockProblem {
     Dune::BitSetVector<Linearization::block_size> truncation;
   };
 
-  MyBlockProblem(Dune::ParameterTree const &parset,
-                 ConvexProblem const &problem)
-      : parset(parset), problem(problem), localBisection() // NOTE: defaults
-  {}
+  MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
+      : 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) {
@@ -90,15 +103,10 @@ template <class ConvexProblem> class MyBlockProblem {
 
     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;
@@ -115,7 +123,7 @@ template <class ConvexProblem> class MyBlockProblem {
     linearization.truncation.resize(u.size());
     linearization.truncation.unsetAll();
     for (size_t i = 0; i < u.size(); ++i) {
-      if (problem.phi.regularity(i, u[i]) > 1e8) { // TODO: Make customisable
+      if (problem_.phi.regularity(i, u[i]) > 1e8) { // TODO: Make customisable
         linearization.truncation[i] = true;
         continue;
       }
@@ -126,70 +134,68 @@ template <class ConvexProblem> class MyBlockProblem {
     }
 
     // construct sparsity pattern for linearization
-    Dune::MatrixIndexSet indices(problem.A.N(), problem.A.M());
-    indices.import(problem.A);
-    problem.phi.addHessianIndices(indices);
+    Dune::MatrixIndexSet indices(problem_.A.N(), problem_.A.M());
+    indices.import(problem_.A);
+    problem_.phi.addHessianIndices(indices);
 
     // construct matrix from pattern and initialize it
     indices.exportIdx(linearization.A);
     linearization.A = 0.0;
 
-    // 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)
+    // compute quadratic part of hessian (linearization.A += problem_.A)
+    for (size_t i = 0; i < problem_.A.N(); ++i) {
+      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;
     }
 
     // compute nonlinearity part of hessian
-    problem.phi.addHessian(u, linearization.A);
+    problem_.phi.addHessian(u, linearization.A);
 
     // compute quadratic part of gradient
     linearization.b.resize(u.size());
-    problem.A.mv(u, linearization.b);
-    linearization.b -= problem.f;
+    problem_.A.mv(u, linearization.b);
+    linearization.b -= problem_.f;
 
     // compute nonlinearity part of gradient
-    problem.phi.addGradient(u, linearization.b);
+    problem_.phi.addGradient(u, linearization.b);
 
     // -grad is needed for Newton step
     linearization.b *= -1.0;
 
     // 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] or linearization
-                    .truncation[col][blockIt.index()])
+            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
-  ConvexProblem const &problem;
+  using Base::problem_;
 
   Bisection const localBisection;
 
@@ -206,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 */
@@ -220,50 +226,48 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
   }
 
   /** \brief Update the i-th block of the current iterate */
-  void updateIterate(LocalVector const &ui, size_t i) {
+  void updateIterate(LocalVectorType const &ui, size_t i) {
     u[i] = ui;
     return;
   }
 
-  /** \brief Minimise a local problem using a modified gradient method
+  /** \brief Minimise a local problem
    * \param[out] ui The solution
    * \param m Block number
    * \param ignore Set of degrees of freedom to leave untouched
    */
   void solveLocalProblem(
-      LocalVector &ui, size_t m,
+      LocalVectorType &ui, size_t m,
       typename Dune::BitSetVector<block_size>::const_reference ignore) {
     {
-      LocalMatrix const *localA = nullptr;
-      LocalVector 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]);
+
+      // 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;
 
-      auto const phi = problem.phi.restriction(m);
-      EllipticEnergy<block_size> localJ(*localA, localb, phi, ignore);
-      minimise(localJ, ui, localsteps, bisection_);
+      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 = \
-	one-body-sample-2D
-
-SOURCES = \
-	assemblers.cc \
-	friction_writer.cc \
-	state/compute_state_dieterich_euler.cc \
-	state/compute_state_ruina.cc \
-	solverfactory.cc \
-	one-body-sample.cc \
-	timestepping.cc \
-	vtk.cc
-
-## 2D
-one_body_sample_2D_SOURCES = $(SOURCES)
-
-one_body_sample_2D_CPPFLAGS = \
-	$(AM_CPPFLAGS) -Dsrcdir=\"$(abs_srcdir)\" -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

@@ -6,14 +6,20 @@
 
 #include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
 #include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
-#include <dune/fufem/assemblers/localassemblers/massassembler.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/viscosityassembler.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/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"
 
@@ -28,95 +34,143 @@ 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);
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleMass(double density, Matrix &M) {
-  MassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
-                Dune::ScaledIdentityMatrix<double, dimension>> const localMass;
-  vertexAssembler.assembleOperator(localMass, M);
-  M *= density;
+void MyAssembler<GridView, dimension>::assembleMass(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+        densityFunction,
+    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);
+  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);
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleViscosity(double shearViscosity,
-                                                         double bulkViscosity,
-                                                         Matrix &C) {
-  ViscosityAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
-  localViscosity(shearViscosity, bulkViscosity);
+void MyAssembler<GridView, dimension>::assembleViscosity(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+    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);
   vertexAssembler.assembleOperator(localViscosity, C);
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleBodyForce(double gravity,
-                                                         double density,
-                                                         LocalVector zenith,
-                                                         Vector &f) {
-  LocalVector weightedGravitationalDirection = zenith;
-  weightedGravitationalDirection *= density * gravity;
-  weightedGravitationalDirection *= -1;
-
-  ConstantFunction<LocalVector, LocalVector> const gravityFunction(
-      weightedGravitationalDirection);
+void MyAssembler<GridView, dimension>::assembleBodyForce(
+    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+    Vector &f) const {
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
-  gravityFunctionalAssembler(gravityFunction);
+      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>
-auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
-    BoundaryPatch<GridView> const &frictionalBoundary, FrictionData const &fd)
-    -> std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>> {
-  // Lump negative normal stress (kludge)
+void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    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;
   {
-    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);
   }
-  for (size_t i = 0; i < weights.size(); ++i)
-    assert(weights[i] >= 0.0);
+  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];
+}
 
-  Dune::BitSetVector<1> frictionalNodes;
-  frictionalBoundary.getVertices(frictionalNodes);
-  return std::make_shared<GlobalRuinaNonlinearity<Matrix, Vector>>(
-      frictionalNodes, weights, fd);
+template <class GridView, int dimension>
+auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+    Config::FrictionModel frictionModel,
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    GlobalFrictionData<dimension> const &frictionInfo,
+    ScalarVector const &weightedNormalStress) const
+    -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
+  // Lumping of the nonlinearity
+  ScalarVector weights;
+  {
+    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+        frictionalBoundaryAssembler(std::make_shared<
+            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+            1));
+    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
+                                               weights, frictionalBoundary);
+  }
+  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,7 +13,10 @@
 #pragma clang diagnostic pop
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
-#include <dune/tectonic/globalruinanonlinearity.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "enums.hh"
 
 template <class GridView, int dimension> class MyAssembler {
 public:
@@ -26,17 +29,18 @@ 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;
   using LocalVector = typename Vector::block_type;
+  using LocalScalarVector = typename ScalarVector::block_type;
 
   using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  GridView const &gridView;
   Assembler<CellBasis, CellBasis> cellAssembler;
   Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
@@ -45,30 +49,42 @@ template <class GridView, int dimension> class MyAssembler {
 
   void assembleFrictionalBoundaryMass(
       BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass);
+      ScalarMatrix &frictionalBoundaryMass) const;
 
-  void assembleMass(double density, Matrix &M);
+  void assembleMass(Dune::VirtualFunction<
+                        LocalVector, LocalScalarVector> const &densityFunction,
+                    Matrix &M) const;
 
-  void assembleElasticity(double E, double nu, Matrix &A);
+  void assembleElasticity(double E, double nu, Matrix &A) const;
 
-  void assembleViscosity(double shearViscosity, double bulkViscosity,
-                         Matrix &C);
+  void assembleViscosity(
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          shearViscosity,
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          bulkViscosity,
+      Matrix &C) const;
 
-  void assembleBodyForce(double gravity, double density, LocalVector zenith,
-                         Vector &f);
+  void assembleBodyForce(
+      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+      Vector &f) const;
 
   void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                        Vector &f,
                        Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime);
+                       double relativeTime) const;
 
-  std::shared_ptr<GlobalRuinaNonlinearity<Matrix, Vector>>
-  assembleFrictionNonlinearity(
+  void assembleWeightedNormalStress(
       BoundaryPatch<GridView> const &frictionalBoundary,
-      FrictionData const &fd);
+      ScalarVector &weightedNormalStress, double youngModulus,
+      double poissonRatio, Vector const &displacement) const;
+
+  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
+      Config::FrictionModel frictionModel,
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      GlobalFrictionData<dimension> const &frictionInfo,
+      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>;