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>
-#include <dune/common/shared_ptr.hh>
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 #include <dune/istl/matrixindexset.hh>
 
-#include "frictionpotential.hh"
-#include "localfriction.hh"
+#include <dune/solvers/common/interval.hh>
 
-namespace Dune {
-template <class Matrix, class Vector> class GlobalNonlinearity {
+#include <dune/tectonic/localfriction.hh>
+
+template <class Matrix, class Vector> class GlobalFriction {
 protected:
-  using SingletonVector = BlockVector<FieldVector<double, 1>>;
+  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
 
 public:
   using IndexSet = Dune::MatrixIndexSet;
@@ -23,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;
   }
@@ -37,28 +35,26 @@ template <class Matrix, class Vector> class GlobalNonlinearity {
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  virtual shared_ptr<LocalFriction<block_size>> 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 &,
-                         Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
 
-  void directionalSubDiff(Vector const &u, Vector const &v,
-                          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) {
-      Interval<double> D;
-      auto const res = restriction(i);
-      res->directionalSubDiff(u[i], v[i], D);
+      Dune::Solvers::Interval<double> D;
+      restriction(i).directionalSubDiff(u[i], v[i], D);
       subdifferential[0] += D[0];
       subdifferential[1] += D[1];
     }
@@ -70,18 +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;
   }
 
-  virtual void updateLogState(SingletonVector 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/nullptr.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 "localfriction.hh"
-#include "frictionpotential.hh"
-
-namespace Dune {
-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>::SingletonVector;
-
-public:
-  GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
-                          SingletonVector const &nodalIntegrals,
-                          FrictionData const &fd)
-      : restrictions(nodalIntegrals.size()) {
-    auto trivialNonlinearity =
-        make_shared<Friction>(make_shared<TrivialFunction>());
-    for (size_t i = 0; i < restrictions.size(); ++i) {
-      restrictions[i] =
-          frictionalNodes[i][0]
-              ? make_shared<Friction>(
-                    make_shared<FrictionPotential>(nodalIntegrals[i], fd))
-              : trivialNonlinearity;
-    }
-  }
-
-  void updateLogState(SingletonVector 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.
-  */
-  shared_ptr<Friction> restriction(size_t i) const override {
-    return restrictions[i];
-  }
-
-private:
-  std::vector<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>
 
 #include <dune/common/fvector.hh>
 #include <dune/common/fmatrix.hh>
-#include <dune/common/shared_ptr.hh>
 
 #include <dune/fufem/arithmetic.hh>
-#include <dune/fufem/interval.hh>
+#include <dune/solvers/common/interval.hh>
 
-#include "frictionpotential.hh"
+#include <dune/tectonic/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;
+template <size_t dimension> class LocalFriction {
+public:
+  virtual ~LocalFriction() {}
 
-  size_t const xsize = x.size();
-  assert(xsize == y.size());
+  using VectorType = Dune::FieldVector<double, dimension>;
+  using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
-  double sum = 0;
-  for (size_t i = 0; i < xsize; ++i)
-    sum += x[i] / xnorm * y[i];
+  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;
 
-  return sum;
-}
+  void virtual addHessian(VectorType const &x, MatrixType &A) const = 0;
 
-namespace Dune {
-template <size_t dimension> class LocalFriction {
-public:
-  using VectorType = FieldVector<double, dimension>;
-  using MatrixType = FieldMatrix<double, dimension, dimension>;
+  void virtual addGradient(VectorType const &x, VectorType &y) const = 0;
+
+  void virtual directionalDomain(
+      VectorType const &, VectorType const &,
+      Dune::Solvers::Interval<double> &dom) const = 0;
+};
 
-  explicit LocalFriction(shared_ptr<FrictionPotentialWrapper> func)
-      : func(func) {}
+template <size_t dimension, class ScalarFriction>
+class WrappedScalarFriction : public LocalFriction<dimension> {
+  using VectorType = typename LocalFriction<dimension>::VectorType;
+  using MatrixType = typename LocalFriction<dimension>::MatrixType;
 
-  double operator()(VectorType const &x) const {
-    return func->evaluate(x.two_norm());
-  }
+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 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,
-                          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] = 0.0; // WARNING: we assume that the outer
-                         // function approaches zero superlinearly
+      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:
@@ -81,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;
@@ -106,23 +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 &,
-                         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:
-  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,43 +9,33 @@
 #include <dune/fufem/interval.hh>
 #include <dune/tnnmg/problem-classes/bisection.hh>
 
-#include "mydirectionalconvexfunction.hh"
+#include <dune/tectonic/mydirectionalconvexfunction.hh>
 
-namespace Dune {
+// Warning: this exploits the property v*x = 0
 template <class Functional>
-void descentMinimisation(Functional const &J,
-                         typename Functional::SmallVector &x,
-                         typename Functional::SmallVector const &v,
-                         Bisection const &bisection) {
-  using SmallVector = typename Functional::SmallVector;
-  using LocalNonlinearity = typename Functional::Nonlinearity;
-
-  MyDirectionalConvexFunction<LocalNonlinearity> const JRest(
-      computeDirectionalA(J.A, v), computeDirectionalb(J.A, J.b, x, v), *J.phi,
-      x, v);
-  // }}}
-
+double lineSearch(Functional const &J,
+                  typename Functional::LocalVector const &x,
+                  typename Functional::LocalVector const &v,
+                  Bisection const &bisection) {
+  MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
+      J.alpha * v.two_norm2(), J.b * v, J.phi, x, v);
   int count;
-  double const stepsize = bisection.minimize(JRest, 0.0, 0.0, count);
-
-  Arithmetic::addProduct(x, stepsize, v);
+  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::SmallVector &x,
-              size_t steps, Bisection const &bisection) {
-  using SmallVector = typename Functional::SmallVector;
-
-  for (size_t step = 0; step < steps; ++step) {
-    SmallVector v;
-    J.gradient(x, v);
-    if (v.two_norm() == 0.0)
-      return;
-
-    v *= -1;
-
-    descentMinimisation(J, x, v, bisection);
-  }
-}
+void minimise(Functional const &J, typename Functional::LocalVector &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>
-
-// Just for debugging
-#include "dune/solvers/computeenergy.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"
-
-/* Just for debugging */
-template <class Matrix, class Vector>
-double computeEnergy(Matrix const &A, Vector const &x, Vector const &b,
-                     Dune::GlobalNonlinearity<Matrix, Vector> const &phi) {
-  return computeEnergy(A, x, b) + phi(x);
-}
+#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>>;
@@ -56,10 +54,16 @@ template <class ConvexProblem> class MyBlockProblem {
     Dune::BitSetVector<Linearization::block_size> truncation;
   };
 
-  MyBlockProblem(Dune::ParameterTree const &parset,
-                 ConvexProblem const &problem)
-      : parset(parset), problem(problem) {
-    bisection = Bisection();
+  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 {
@@ -73,6 +77,11 @@ template <class ConvexProblem> class MyBlockProblem {
     return s;
   }
 
+  double computeEnergy(const VectorType &v) const {
+    return 0.0; // FIXME
+    // return ::computeEnergy(problem_.A, v, problem_.f) + problem_.phi(v);
+  }
+
   void projectCoarseCorrection(VectorType const &u,
                                typename Linearization::VectorType const &v,
                                VectorType &projected_v,
@@ -94,20 +103,15 @@ template <class ConvexProblem> class MyBlockProblem {
 
     v /= vnorm; // Rescale for numerical stability
 
-    MyDirectionalConvexFunction<
-        Dune::GlobalNonlinearity<MatrixType, VectorType>> const
-    psi(computeDirectionalA(problem.A, v),
-        computeDirectionalb(problem.A, problem.f, u, v), problem.phi, u, v);
-
-    Interval<double> D;
+    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;
-    Bisection bisection;
-    return bisection.minimize(psi, vnorm, 0.0, bisectionsteps) / vnorm; // TODO
+    Bisection const globalBisection; // NOTE: defaults
+    return globalBisection.minimize(psi, vnorm, 0.0, bisectionsteps) / vnorm;
   }
 
   void assembleTruncate(VectorType const &u, Linearization &linearization,
@@ -119,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;
       }
@@ -130,73 +134,70 @@ 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, bisection, problem);
+    return IterateObject(localBisection, problem_, maxEigenvalues_);
   }
 
 private:
-  Dune::ParameterTree const &parset;
+  std::vector<double> maxEigenvalues_;
 
   // problem data
-  ConvexProblem const &problem;
+  using Base::problem_;
 
-  // commonly used minimization stuff
-  Bisection bisection;
+  Bisection const localBisection;
 
   mutable std::ostringstream outStream;
 };
@@ -211,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 */
@@ -225,67 +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) {
     {
-      size_t ic =
-          block_size; // Special value that indicates nothing should be ignored
-      switch (ignore.count()) {
-        case 0: // Full problem
-          break;
-        case 1:
-          for (ic = 0; ic < block_size; ++ic)
-            if (ignore[ic])
-              break;
-          break;
-        case block_size: // Ignore the whole node
-          return;
-        default:
-          assert(false);
+      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();
+        Arithmetic::subtractProduct(localb, *it, u[j]); // also the diagonal!
       }
+      Arithmetic::addProduct(localb, maxEigenvalues_[m], u[m]);
 
-      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) {
-        size_t const j = it.index();
-        if (j == m)
-          localA = &(*it); // localA = A[m][m]
+      // We minimise over an affine subspace
+      for (size_t j = 0; j < block_size; ++j)
+        if (ignore[j])
+          localb[j] = 0;
         else
-          Arithmetic::subtractProduct(localb, *it, u[j]);
-      }
-      assert(localA != nullptr);
+          ui[j] = 0;
 
-      auto const phi = problem.phi.restriction(m);
-      Dune::EllipticEnergy<block_size> localJ(*localA, localb, phi, ic);
-      Dune::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;
-
-  // commonly used minimization stuff
-  Bisection bisection;
-
+  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/fufem/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);
-}
+#include <dune/solvers/common/interval.hh>
 
 template <class Nonlinearity> class MyDirectionalConvexFunction {
 public:
@@ -41,15 +22,16 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
 
   double linearPart() const { return b; }
 
-  void subDiff(double x, Interval<double> &D) const {
+  void subDiff(double x, Dune::Solvers::Interval<double> &D) const {
     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(Interval<double> &domain) const {
+  void domain(Dune::Solvers::Interval<double> &domain) const {
     domain[0] = this->dom[0];
     domain[1] = this->dom[1];
   }
@@ -62,7 +44,25 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
   Vector const &u;
   Vector const &v;
 
-  Interval<double> dom;
+  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 \
-	one-body-sample-3D
-
-2d: one-body-sample-2D
-3d: one-body-sample-3D
-
-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
-
-## 3D
-one_body_sample_3D_SOURCES = $(SOURCES)
-
-one_body_sample_3D_CPPFLAGS = \
-	$(AM_CPPFLAGS) -Dsrcdir=\"$(abs_srcdir)\" -DDIM=3
-
-# Some are for clang, others are for gcc
-AM_CXXFLAGS = \
-	-Wall \
-	-Wextra \
-	-Wno-unused-parameter \
-	-Wno-overloaded-virtual
-
-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

@@ -2,56 +2,181 @@
 #include "config.h"
 #endif
 
+#include <dune/istl/scaledidmatrix.hh>
+
+#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/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
-#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/functiontools/p0p1interpolation.hh>
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
-#include <dune/tectonic/globalruinanonlinearity.hh>
+#include <dune/tectonic/frictionpotential.hh>
+#include <dune/tectonic/globalratestatefriction.hh>
 
 #include "assemblers.hh"
 
-// Assembles Neumann boundary term in f
-template <class GridView, class LocalVector, class Assembler>
-void assembleNeumann(GridView const &gridView, Assembler const &assembler,
-                     Dune::BitSetVector<1> const &neumannNodes,
-                     Dune::BlockVector<LocalVector> &f,
-                     Dune::VirtualFunction<double, double> const &neumann,
-                     double relativeTime) { // constant sample function on
-                                            // neumann boundary
-  BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
+template <class GridView, int dimension>
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
+    : cellBasis(_gridView),
+      vertexBasis(_gridView),
+      gridView(_gridView),
+      cellAssembler(cellBasis, cellBasis),
+      vertexAssembler(vertexBasis, vertexBasis) {}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    ScalarMatrix &frictionalBoundaryMass) const {
+  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
+                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
+      frictionalBoundaryMassAssembler(frictionalBoundary);
+  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
+                                   frictionalBoundaryMass);
+}
+
+template <class GridView, int dimension>
+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) const {
+  StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
+      localStiffness(E, nu);
+  vertexAssembler.assembleOperator(localStiffness, A);
+}
+
+template <class GridView, int dimension>
+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(
+    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+    Vector &f) const {
+  L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
+      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) const {
   LocalVector localNeumann(0);
   neumann.evaluate(relativeTime, localNeumann[0]);
-  ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
-  NeumannBoundaryAssembler<typename GridView::Grid, LocalVector>
-  neumannBoundaryAssembler(fNeumann);
-  assembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
-                                       neumannBoundary);
-}
-
-// Assembles constant 1-function on frictional boundary in nodalIntegrals
-template <class GridView, class LocalVector, class Assembler>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
-                          Dune::BitSetVector<1> const &frictionalNodes) {
-  using Singleton = Dune::FieldVector<double, 1>;
-  BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
-  ConstantFunction<LocalVector, Singleton> const constantOneFunction(1);
-  NeumannBoundaryAssembler<typename GridView::Grid, Singleton>
-  frictionalBoundaryAssembler(constantOneFunction);
-
-  auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
-  assembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
-                                       *nodalIntegrals, frictionalBoundary);
-  return nodalIntegrals;
-}
-
-template <class Matrix, class Vector>
-Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
-    Dune::BitSetVector<1> const &frictionalNodes,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd) {
-  return Dune::make_shared<Dune::GlobalRuinaNonlinearity<Matrix, Vector>>(
-      frictionalNodes, nodalIntegrals, fd);
+  NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
+      std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
+          localNeumann));
+  vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
+                                             neumannBoundary);
+}
+
+template <class GridView, int dimension>
+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;
+  {
+    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 &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) const {
+  auto const gridDisplacement =
+      std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
+          vertexBasis, u);
+  VonMisesStressAssembler<Grid, LocalCellBasis> localStressAssembler(
+      youngModulus, poissonRatio, gridDisplacement);
+  cellAssembler.assembleFunctional(localStressAssembler, stress);
 }
 
 #include "assemblers_tmpl.cc"

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>
-#include <dune/common/fvector.hh>
-#include <dune/common/shared_ptr.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
 #include <dune/fufem/assemblers/assembler.hh>
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wsign-compare"
+#include <dune/fufem/functionspacebases/p0basis.hh>
+#pragma clang diagnostic pop
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
-#include <dune/tectonic/globalnonlinearity.hh>
-
-template <class GridView, class LocalVector, class Assembler>
-void assembleNeumann(GridView const &gridView, Assembler const &assembler,
-                     Dune::BitSetVector<1> const &neumannNodes,
-                     Dune::BlockVector<LocalVector> &f,
-                     Dune::VirtualFunction<double, double> const &neumann,
-                     double relativeTime);
-
-template <class GridView, class LocalVector, class Assembler>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
-                          Dune::BitSetVector<1> const &frictionalNodes);
-
-template <class Matrix, class Vector>
-Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
-    Dune::BitSetVector<1> const &frictionalNodes,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd);
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "enums.hh"
+
+template <class GridView, int dimension> class MyAssembler {
+public:
+  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+  using Matrix =
+      Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
+  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+
+  using CellBasis = P0Basis<GridView, double>;
+  using VertexBasis = P1NodalBasis<GridView, double>;
+
+  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;
+
+  Assembler<CellBasis, CellBasis> cellAssembler;
+  Assembler<VertexBasis, VertexBasis> vertexAssembler;
+
+public:
+  MyAssembler(GridView const &gridView);
+
+  void assembleFrictionalBoundaryMass(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarMatrix &frictionalBoundaryMass) const;
+
+  void assembleMass(Dune::VirtualFunction<
+                        LocalVector, LocalScalarVector> const &densityFunction,
+                    Matrix &M) const;
+
+  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) const;
+
+  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) const;
+
+  void assembleWeightedNormalStress(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      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) const;
+};
 #endif

src/assemblers_tmpl.cc

-#ifndef DIM
-#error DIM unset
+#ifndef MY_DIM
+#error MY_DIM unset
 #endif
 
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
-#include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
-#include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
+#include "explicitgrid.hh"
 
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-using SmallVector = Dune::FieldVector<double, DIM>;
-using SmallMatrix = Dune::FieldMatrix<double, DIM, DIM>;
-using Matrix = Dune::BCRSMatrix<SmallMatrix>;
-using Vector = Dune::BlockVector<SmallVector>;
-
-using Grid = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
-using GridView = Grid::LeafGridView;
-using P1Basis = P1NodalBasis<GridView, double>;
-using MyAssembler = Assembler<P1Basis, P1Basis>;
-
-template void assembleNeumann<GridView, SmallVector, MyAssembler>(
-    GridView const &gridView, MyAssembler const &assembler,
-    Dune::BitSetVector<1> const &neumannNodes,
-    Dune::BlockVector<SmallVector> &f,
-    Dune::VirtualFunction<double, double> const &neumann, double relativeTime);
-
-template Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assembleFrictionWeightsal<GridView, SmallVector, MyAssembler>(
-    GridView const &gridView, MyAssembler const &assembler,
-    Dune::BitSetVector<1> const &frictionalNodes);
-
-template Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>>
-assembleNonlinearity<Matrix, Vector>(
-    Dune::BitSetVector<1> const &frictionalNodes,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    FrictionData const &fd);
+template class MyAssembler<GridView, MY_DIM>;