src/one-body-problem-3D.cfg

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter= 2e-2  # [m]
+
+[boundary.friction.weakening]
+patchType       = Trapezoidal
+
+[timeSteps]
+refinementTolerance = 1e-5
+
+[u0.solver]
+tolerance         = 1e-6
+
+[a0.solver]
+tolerance         = 1e-6
+
+[v.solver]
+tolerance         = 1e-6
+
+[v.fpi]
+tolerance         = 1e-5
+
+[solver.tnnmg.linear]
+tolerance          = 1e-10

src/one-body-problem-data/bc.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+
+class VelocityDirichletCondition
+    : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const {
+    // Assumed to vanish at time zero
+    double const finalVelocity = -5e-5;
+    y = (relativeTime <= 0.1)
+            ? finalVelocity * (1.0 - std::cos(relativeTime * M_PI / 0.1)) / 2.0
+            : finalVelocity;
+  }
+};
+
+class NeumannCondition : public Dune::VirtualFunction<double, double> {
+  void evaluate(double const &relativeTime, double &y) const { y = 0.0; }
+};
+#endif

src/one-body-problem-data/geometry.tex

+\documentclass[tikz]{minimal}
+
+\usepackage{tikz}
+\usetikzlibrary{calc}
+\usetikzlibrary{decorations.pathreplacing}
+
+\usepackage{siunitx}
+
+\begin{document}
+\pgfmathsetmacro{\rightleg}{0.27}
+\pgfmathsetmacro{\leftleg}{1.00}
+\pgfmathsetmacro{\leftangle}{atan(\rightleg/\leftleg)}
+\begin{tikzpicture}[scale=12, rotate=\leftangle]
+  \pgfmathsetmacro{\mysin}{sin(\leftangle)}
+  \pgfmathsetmacro{\mycos}{cos(\leftangle)}
+  \pgfmathsetmacro{\viscoheight}{0.06}
+  \pgfmathsetmacro{\Zx}{0.35}
+  \pgfmathsetmacro{\weaklen}{0.20}
+
+  \coordinate (A) at (0,0);
+  \node at (A) [left] {A};
+  \coordinate (B) at (\leftleg,-\rightleg);
+  \node at (B) [right] {B};
+  \coordinate (C) at (\leftleg,0);
+  \node at (C) [right] {C};
+
+  \draw (A) -- (B) -- (C) -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
+
+  \coordinate (Z) at (\Zx,0);
+  \node at (Z) [above] {Z};
+  \coordinate (Y) at ($(\Zx,-\Zx/\leftleg * \rightleg)$);
+  \node at (Y) [below] {Y};
+  \coordinate (X) at ($(Y) + (-\weaklen*\mycos,\weaklen*\mysin)$);
+  \node at (X) [below] {X};
+  \path let \p1 = (X) in coordinate (U) at ($(\x1, 0)$);
+  \node at (U) [above] {U};
+
+  \path (A) -- node [above=.25cm, sloped] {$\overline{AZ} = \SI{35}{cm}$} (Z);
+
+  \draw[color=red, thick] (X) -- node [below=.25cm] {$\overline{XY}=\SI{20}{cm}$} (Y);
+  \draw[dashed] (Y) -- (Z);
+  \draw[dashed] (U) -- (X);
+
+  \coordinate (K) at ($(B) + (-\leftleg * \viscoheight / \rightleg,\viscoheight)$);
+  \node at (K) [below] {K};
+  \coordinate (M) at ($(B) + (0, \viscoheight)$);
+  \node at (M) [right] {M};
+  \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
+
+  \path[fill=blue] (K) -- (B) -- node [right=.75cm] {$\overline{BM}=\SI{6}{cm}$} (M) -- cycle;
+
+  \coordinate (G) at ($(A) ! 0.5 ! (X)$);
+  \node at (G) [below] {G};
+  \coordinate (H) at ($(X) ! 0.5 ! (Y)$);
+  \node at (H) [below] {H};
+  \coordinate (J) at ($(Y) ! 0.5 ! (B)$);
+  \node at (J) [below] {J};
+
+  \coordinate (I) at ($(Y) + (G)$);
+  \node at (I) [below] {I};
+
+  \node[align=left] at (0.5,-0.225) {
+    $Z$: coast line\\
+    $\overline{XY}$: velocity weakening zone\\
+    $BKM$: visco-elastic domain};
+\end{tikzpicture}
+
+\end{document}

src/one-body-problem-data/midpoint.hh

+#ifndef SRC_MIDPOINT_HH
+#define SRC_MIDPOINT_HH
+
+#include <dune/solvers/common/arithmetic.hh>
+
+template <class Vector> Vector midPoint(Vector const &x, Vector const &y) {
+  Vector ret(0);
+  Arithmetic::addProduct(ret, 0.5, x);
+  Arithmetic::addProduct(ret, 0.5, y);
+  return ret;
+}
+#endif

src/mybody.hh → src/one-body-problem-data/mybody.hh

-#ifndef MY_BODY_HH
-#define MY_BODY_HH
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
 
 #include <dune/common/fvector.hh>
 
@@ -9,30 +9,36 @@
 #include <dune/tectonic/gravity.hh>
 
 #include "mygeometry.hh"
+#include "segmented-function.hh"
 
 template <int dimension> class MyBody : public Body<dimension> {
   using typename Body<dimension>::ScalarFunction;
   using typename Body<dimension>::VectorField;
-  using MyConstantFunction = ConstantFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
 
 public:
-  MyBody(Dune::ParameterTree const &parset, MyGeometry const &mygeometry)
+  MyBody(Dune::ParameterTree const &parset)
       : poissonRatio_(parset.get<double>("body.poissonRatio")),
-        youngModulus_(parset.get<double>("body.youngModulus")),
-        shearViscosity_(parset.get<double>("body.shearViscosity")),
-        bulkViscosity_(parset.get<double>("body.bulkViscosity")),
-        densityField_(parset.get<double>("body.density")),
-        gravityField_(densityField_, mygeometry.zenith,
+        youngModulus_(3.0 * parset.get<double>("body.bulkModulus") *
+                      (1.0 - 2.0 * poissonRatio_)),
+        shearViscosityField_(
+            parset.get<double>("body.elastic.shearViscosity"),
+            parset.get<double>("body.viscoelastic.shearViscosity")),
+        bulkViscosityField_(
+            parset.get<double>("body.elastic.bulkViscosity"),
+            parset.get<double>("body.viscoelastic.bulkViscosity")),
+        densityField_(parset.get<double>("body.elastic.density"),
+                      parset.get<double>("body.viscoelastic.density")),
+        gravityField_(densityField_, MyGeometry::zenith,
                       parset.get<double>("gravity")) {}
 
   double getPoissonRatio() const override { return poissonRatio_; }
-
   double getYoungModulus() const override { return youngModulus_; }
-
-  double getShearViscosity() const override { return shearViscosity_; }
-  double getBulkViscosity() const override { return bulkViscosity_; }
-
+  ScalarFunction const &getShearViscosityField() const override {
+    return shearViscosityField_;
+  }
+  ScalarFunction const &getBulkViscosityField() const override {
+    return bulkViscosityField_;
+  }
   ScalarFunction const &getDensityField() const override {
     return densityField_;
   }
@@ -41,10 +47,9 @@ template <int dimension> class MyBody : public Body<dimension> {
 private:
   double const poissonRatio_;
   double const youngModulus_;
-  double const shearViscosity_;
-  double const bulkViscosity_;
-  MyConstantFunction const densityField_;
+  SegmentedFunction const shearViscosityField_;
+  SegmentedFunction const bulkViscosityField_;
+  SegmentedFunction const densityField_;
   Gravity<dimension> const gravityField_;
 };
-
 #endif

src/one-body-problem-data/mygeometry.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+#include "mygeometry.hh"
+
+void MyGeometry::write() {
+  std::fstream writer("geometry", std::fstream::out);
+  writer << "A = " << A << std::endl;
+  writer << "B = " << B << std::endl;
+  writer << "C = " << C << std::endl;
+  writer << "Y = " << Y << std::endl;
+  writer << "X = " << X << std::endl;
+  writer << "Z = " << Z << std::endl;
+  writer << "U = " << U << std::endl;
+  writer << "K = " << K << std::endl;
+  writer << "M = " << M << std::endl;
+  writer << "G = " << G << std::endl;
+  writer << "H = " << H << std::endl;
+  writer << "J = " << J << std::endl;
+  writer << "I = " << I << std::endl;
+  writer << "zenith = " << zenith << std::endl;
+}
+
+void MyGeometry::render() {
+#ifdef HAVE_CAIROMM
+  std::string const filename = "geometry.png";
+  double const width = 600;
+  double const height = 400;
+  double const widthScale = 400;
+  double const heightScale = 400;
+
+  auto surface =
+      Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, width, height);
+  auto cr = Cairo::Context::create(surface);
+
+  auto const setRGBColor = [&](int colour) {
+    cr->set_source_rgb(((colour & 0xFF0000) >> 16) / 255.0,
+                       ((colour & 0x00FF00) >> 8) / 255.0,
+                       ((colour & 0x0000FF) >> 0) / 255.0);
+  };
+  auto const moveTo = [&](LocalVector2D const &v) { cr->move_to(v[0], -v[1]); };
+  auto const lineTo = [&](LocalVector2D const &v) { cr->line_to(v[0], -v[1]); };
+
+  cr->scale(widthScale, heightScale);
+  cr->translate(0.1, 0.1);
+  cr->set_line_width(0.0025);
+
+  // triangle
+  {
+    moveTo(reference::A);
+    lineTo(reference::B);
+    lineTo(reference::C);
+    cr->close_path();
+    cr->stroke();
+  }
+
+  // dashed lines
+  {
+    cr->save();
+    std::vector<double> dashPattern = { 0.005 };
+    cr->set_dash(dashPattern, 0);
+    moveTo(reference::Z);
+    lineTo(reference::Y);
+    moveTo(reference::U);
+    lineTo(reference::X);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // fill viscoelastic region
+  {
+    cr->save();
+    setRGBColor(0x0097E0);
+    moveTo(reference::B);
+    lineTo(reference::K);
+    lineTo(reference::M);
+    cr->fill();
+    cr->restore();
+  }
+
+  // mark weakening region
+  {
+    cr->save();
+    setRGBColor(0x7AD3FF);
+    cr->set_line_width(0.005);
+    moveTo(reference::X);
+    lineTo(reference::Y);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // mark points
+  {
+    auto const drawCircle = [&](LocalVector2D const &v) {
+      cr->arc(v[0], -v[1], 0.0075, -M_PI, M_PI); // x,y,radius,angle1,angle2
+      cr->fill();
+    };
+
+    cr->save();
+    setRGBColor(0x002F47);
+    drawCircle(reference::A);
+    drawCircle(reference::B);
+    drawCircle(reference::C);
+    drawCircle(reference::Y);
+    drawCircle(reference::X);
+    drawCircle(reference::Z);
+    drawCircle(reference::U);
+    drawCircle(reference::K);
+    drawCircle(reference::M);
+    drawCircle(reference::G);
+    drawCircle(reference::H);
+    drawCircle(reference::J);
+    drawCircle(reference::I);
+    cr->restore();
+  }
+
+  // labels
+  {
+    auto const label = [&](LocalVector2D const &v, std::string l) {
+      moveTo(v);
+      cr->rel_move_to(0.005, -0.02);
+      cr->show_text(l);
+    };
+    auto font = Cairo::ToyFontFace::create(
+        "monospace", Cairo::FONT_SLANT_NORMAL, Cairo::FONT_WEIGHT_NORMAL);
+
+    cr->save();
+    cr->set_font_face(font);
+    cr->set_font_size(0.03);
+
+    label(reference::A, "A");
+    label(reference::B, "B");
+    label(reference::C, "C");
+    label(reference::K, "K");
+    label(reference::M, "M");
+    label(reference::U, "U");
+    label(reference::X, "X");
+    label(reference::Y, "Y");
+    label(reference::Z, "Z");
+    label(reference::G, "G");
+    label(reference::H, "H");
+    label(reference::J, "J");
+    label(reference::I, "I");
+    cr->restore();
+  }
+
+  surface->write_to_png(filename);
+#endif
+}

src/one-body-problem-data/mygeometry.hh

+#ifndef SRC_MYGEOMETRY_HH
+#define SRC_MYGEOMETRY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "midpoint.hh"
+
+namespace MyGeometry {
+namespace {
+  using LocalVector2D = Dune::FieldVector<double, 2>;
+  using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;
+
+  using LocalVector = Dune::FieldVector<double, MY_DIM>;
+}
+
+namespace reference {
+  double const s = 1.0; // scaling factor
+
+  double const rightLeg = 0.27 * s;
+  double const leftLeg = 1.00 * s;
+  double const leftAngle = atan(rightLeg / leftLeg);
+  double const viscoHeight = 0.06 * s; // Height of the viscous bottom layer
+  double const weakLen = 0.20 * s;     // Length of the weak zone
+
+  double const zDistance = 0.35;
+
+  LocalVector2D const A = {0, 0};
+  LocalVector2D const B = {leftLeg, -rightLeg};
+  LocalVector2D const C = {leftLeg, 0};
+
+  LocalVector2D const Z = {zDistance * s, 0};
+  LocalVector2D const Y = {zDistance * s, -zDistance *s / leftLeg *rightLeg};
+  LocalVector2D const X = {Y[0] - weakLen * std::cos(leftAngle),
+                           Y[1] + weakLen *std::sin(leftAngle)};
+
+  LocalVector2D const U = {X[0], 0};
+
+  LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
+                           B[1] + viscoHeight};
+  LocalVector2D const M = {B[0], B[1] + viscoHeight};
+
+  LocalVector2D const G = midPoint(A, X);
+  LocalVector2D const H = midPoint(X, Y);
+  LocalVector2D const J = midPoint(Y, B);
+
+  LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
+
+  LocalVector2D const zenith = {0, 1};
+
+  LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
+                                  {std::sin(leftAngle), std::cos(leftAngle)}};
+}
+
+namespace {
+  LocalVector rotate(LocalVector2D const &x) {
+    LocalVector2D ret2D;
+    reference::rotation.mv(x, ret2D);
+    LocalVector ret(0);
+    ret[0] = ret2D[0];
+    ret[1] = ret2D[1];
+    return ret;
+  }
+}
+
+double const lengthScale = reference::s;
+
+double const depth = 0.60 * lengthScale;
+
+LocalVector const A = rotate(reference::A);
+LocalVector const B = rotate(reference::B);
+LocalVector const C = rotate(reference::C);
+LocalVector const G = rotate(reference::G);
+LocalVector const H = rotate(reference::H);
+LocalVector const I = rotate(reference::I);
+LocalVector const J = rotate(reference::J);
+LocalVector const K = rotate(reference::K);
+LocalVector const M = rotate(reference::M);
+LocalVector const U = rotate(reference::U);
+LocalVector const X = rotate(reference::X);
+LocalVector const Y = rotate(reference::Y);
+LocalVector const Z = rotate(reference::Z);
+
+LocalVector const zenith = rotate(reference::zenith);
+
+void write();
+
+void render();
+}
+#endif

src/one-body-problem-data/myglobalfrictiondata.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+
+#include <dune/common/function.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "patchfunction.hh"
+
+template <class LocalVector>
+class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
+private:
+  using typename GlobalFrictionData<LocalVector::dimension>::VirtualFunction;
+
+public:
+  MyGlobalFrictionData(Dune::ParameterTree const &parset,
+                       ConvexPolyhedron<LocalVector> const &segment)
+      : C_(parset.get<double>("C")),
+        L_(parset.get<double>("L")),
+        V0_(parset.get<double>("V0")),
+        a_(parset.get<double>("strengthening.a"),
+           parset.get<double>("weakening.a"), segment),
+        b_(parset.get<double>("strengthening.b"),
+           parset.get<double>("weakening.b"), segment),
+        mu0_(parset.get<double>("mu0")) {}
+
+  double const &C() const override { return C_; }
+  double const &L() const override { return L_; }
+  double const &V0() const override { return V0_; }
+  VirtualFunction const &a() const override { return a_; }
+  VirtualFunction const &b() const override { return b_; }
+  double const &mu0() const override { return mu0_; }
+
+private:
+  double const C_;
+  double const L_;
+  double const V0_;
+  PatchFunction const a_;
+  PatchFunction const b_;
+  double const mu0_;
+};
+#endif

src/one-body-problem-data/mygrid.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "mygrid.hh"
+#include "midpoint.hh"
+#include "../diameter.hh"
+
+#if MY_DIM == 3
+SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
+#endif
+
+// back-to-front, front-to-back, front-to-back
+void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
+                                        unsigned int W) {
+#if MY_DIM == 3
+  unsigned int const U2 = U + shift_;
+  unsigned int const V2 = V + shift_;
+  unsigned int const W2 = W + shift_;
+
+  simplices_.push_back({ U, V, W, U2 });
+  simplices_.push_back({ V, V2, W2, U2 });
+  simplices_.push_back({ W, W2, U2, V });
+#else
+  simplices_.push_back({ U, V, W });
+#endif
+}
+
+// back-to-front, back-to-front, front-to-back
+void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
+                                        unsigned int W) {
+#if MY_DIM == 3
+  unsigned int const U2 = U + shift_;
+  unsigned int const V2 = V + shift_;
+  unsigned int const W2 = W + shift_;
+
+  simplices_.push_back({ U, V, W, U2 });
+  simplices_.push_back({ V, V2, W, U2 });
+  simplices_.push_back({ V2, W, U2, W2 });
+#else
+  simplices_.push_back({ U, V, W });
+#endif
+}
+
+auto SimplexManager::getSimplices() -> SimplexList const & {
+  return simplices_;
+}
+
+template <class Grid> GridConstructor<Grid>::GridConstructor() {
+  auto const &A = MyGeometry::A;
+  auto const &B = MyGeometry::B;
+  auto const &C = MyGeometry::C;
+
+  unsigned int const vc = 3;
+
+#if MY_DIM == 3
+  Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
+#else
+  Dune::FieldMatrix<double, vc, MY_DIM> vertices;
+#endif
+  for (size_t i = 0; i < 2; ++i) {
+#if MY_DIM == 3
+    size_t numXYplanes = 2;
+#else
+    size_t numXYplanes = 1;
+#endif
+    size_t k = 0;
+    for (size_t j = 1; j <= numXYplanes; ++j) {
+      vertices[k++][i] = A[i];
+      vertices[k++][i] = B[i];
+      vertices[k++][i] = C[i];
+      assert(k == j * vc);
+    }
+  }
+
+#if MY_DIM == 3
+  for (size_t k = 0; k < vc; ++k) {
+    vertices[k][2] = -MyGeometry::depth / 2.0;
+    vertices[k + vc][2] = MyGeometry::depth / 2.0;
+  }
+#endif
+
+  for (size_t i = 0; i < vertices.N(); ++i)
+    gridFactory.insertVertex(vertices[i]);
+
+  Dune::GeometryType cell;
+#if MY_DIM == 3
+  cell.makeTetrahedron();
+#else
+  cell.makeTriangle();
+#endif
+
+#if MY_DIM == 3
+  SimplexManager sm(vc);
+#else
+  SimplexManager sm;
+#endif
+  sm.addFromVerticesFFB(1, 2, 0);
+  auto const &simplices = sm.getSimplices();
+
+  // sanity-check choices of simplices
+  for (size_t i = 0; i < simplices.size(); ++i) {
+    Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
+    for (size_t j = 0; j < MY_DIM; ++j)
+      check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+    assert(check.determinant() > 0);
+    gridFactory.insertElement(cell, simplices[i]);
+  }
+}
+
+template <class Grid> std::shared_ptr<Grid> GridConstructor<Grid>::getGrid() {
+  return std::shared_ptr<Grid>(gridFactory.createGrid());
+}
+
+template <class Grid>
+template <class GridView>
+MyFaces<GridView> GridConstructor<Grid>::constructFaces(
+    GridView const &gridView) {
+  return MyFaces<GridView>(gridView);
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+                                    Vector const &c) {
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+                                Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 * MyGeometry::lengthScale > x[0] or
+      x[0] > minmax0.second + 1e-14 * MyGeometry::lengthScale)
+    return false;
+  if (minmax1.first - 1e-14 * MyGeometry::lengthScale > x[1] or
+      x[1] > minmax1.second + 1e-14 * MyGeometry::lengthScale)
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView)
+    :
+#if MY_DIM == 3
+      lower(gridView),
+      right(gridView),
+      upper(gridView),
+      front(gridView),
+      back(gridView)
+#else
+      lower(gridView),
+      right(gridView),
+      upper(gridView)
+#endif
+{
+  assert(isClose(MyGeometry::A[1], 0));
+  assert(isClose(MyGeometry::B[1], 0));
+  lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return isClose(0, in.geometry().center()[1]);
+  });
+#if MY_DIM == 3
+  front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return isClose(MyGeometry::depth / 2.0, in.geometry().center()[2]);
+  });
+
+  back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return isClose(-MyGeometry::depth / 2.0, in.geometry().center()[2]);
+  });
+#endif
+
+  upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return xyBetween(MyGeometry::A, MyGeometry::C, in.geometry().center());
+  });
+
+  right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+    return xyBetween(MyGeometry::B, MyGeometry::C, in.geometry().center());
+  });
+}
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter) {
+  return (distance / 0.0125 / MyGeometry::lengthScale + 1.0) * smallestDiameter;
+}
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter) {
+  bool needRefine = true;
+  while (true) {
+    needRefine = false;
+    for (auto &&e : elements(grid.leafGridView())) {
+      auto const geometry = e.geometry();
+
+      auto const weakeningRegionDistance =
+          distance(weakPatch, geometry, 1e-6 * MyGeometry::lengthScale);
+      auto const admissibleDiameter =
+          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter);
+
+      if (diameter(geometry) <= admissibleDiameter)
+        continue;
+
+      needRefine = true;
+      grid.mark(1, e);
+    }
+    if (!needRefine)
+      break;
+
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+  }
+}
+
+#include "mygrid_tmpl.cc"

src/one-body-problem-data/mygrid.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include "mygeometry.hh"
+
+template <class GridView> struct MyFaces {
+  BoundaryPatch<GridView> lower;
+  BoundaryPatch<GridView> right;
+  BoundaryPatch<GridView> upper;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> front;
+  BoundaryPatch<GridView> back;
+#endif
+
+  MyFaces(GridView const &gridView);
+
+private:
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14 * MyGeometry::lengthScale;
+  };
+
+  bool isClose2(double a, double b) {
+    return std::abs(a - b) <
+           1e-14 * MyGeometry::lengthScale * MyGeometry::lengthScale;
+  };
+
+  template <class Vector>
+  bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+  template <class Vector>
+  bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+  template <class Vector>
+  bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
+};
+
+class SimplexManager {
+public:
+  using SimplexList = std::vector<std::vector<unsigned int>>;
+
+#if MY_DIM == 3
+  SimplexManager(unsigned int shift);
+#endif
+
+  void addFromVerticesFBB(unsigned int U, unsigned int V, unsigned int W);
+  void addFromVerticesFFB(unsigned int U, unsigned int V, unsigned int W);
+
+  SimplexList const &getSimplices();
+
+private:
+  SimplexList simplices_;
+
+#if MY_DIM == 3
+  unsigned int const shift_;
+#endif
+};
+
+template <class Grid> class GridConstructor {
+public:
+  GridConstructor();
+
+  std::shared_ptr<Grid> getGrid();
+
+  template <class GridView>
+  MyFaces<GridView> constructFaces(GridView const &gridView);
+
+private:
+  Dune::GridFactory<Grid> gridFactory;
+};
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter);
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter);
+
+#endif

src/one-body-problem-data/mygrid_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+template class GridConstructor<Grid>;
+
+template struct MyFaces<GridView>;
+
+template MyFaces<GridView> GridConstructor<Grid>::constructFaces(
+    GridView const &gridView);
+
+template void refine<Grid, LocalVector>(
+    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+    double smallestDiameter);

src/one-body-problem-data/patchfunction.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+class PatchFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  using Polyhedron = ConvexPolyhedron<Dune::FieldVector<double, MY_DIM>>;
+
+  double const v1_;
+  double const v2_;
+  Polyhedron const &segment_;
+
+public:
+  PatchFunction(double v1, double v2, Polyhedron const &segment)
+      : v1_(v1), v2_(v2), segment_(segment) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+                                                                      : v1_;
+  }
+};
+
+#endif

src/one-body-problem-data/segmented-function.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include "mygeometry.hh"
+
+class SegmentedFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
+                 Dune::FieldVector<double, MY_DIM> const &y,
+                 Dune::FieldVector<double, MY_DIM> const &z) const {
+    return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
+  };
+  bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
+    return liesBelow(MyGeometry::K, MyGeometry::M, z);
+  };
+
+  double const _v1;
+  double const _v2;
+
+public:
+  SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = insideRegion2(x) ? _v2 : _v1;
+  }
+};
+#endif

src/one-body-problem-data/weakpatch.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
+
+template <class LocalVector>
+ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
+  ConvexPolyhedron<LocalVector> weakPatch;
+#if MY_DIM == 3
+  weakPatch.vertices.resize(4);
+  weakPatch.vertices[0] = weakPatch.vertices[2] = MyGeometry::X;
+  weakPatch.vertices[1] = weakPatch.vertices[3] = MyGeometry::Y;
+  for (size_t k = 0; k < 2; ++k) {
+    weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
+    weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
+  }
+  switch (parset.get<Config::PatchType>("patchType")) {
+    case Config::Rectangular:
+      break;
+    case Config::Trapezoidal:
+      weakPatch.vertices[1][0] += 0.05 * MyGeometry::lengthScale;
+      weakPatch.vertices[3][0] -= 0.05 * MyGeometry::lengthScale;
+      break;
+    default:
+      assert(false);
+  }
+#else
+  weakPatch.vertices.resize(2);
+  weakPatch.vertices[0] = MyGeometry::X;
+  weakPatch.vertices[1] = MyGeometry::Y;
+#endif
+  return weakPatch;
+};
+#endif

src/one-body-problem.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/formatstring.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/solvers/solver.hh>
+#include <dune/tnnmg/problem-classes/convexproblem.hh>
+
+#include <dune/tectonic/geocoordinate.hh>
+#include <dune/tectonic/myblockproblem.hh>
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include "assemblers.hh"
+#include "diameter.hh"
+#include "enumparser.hh"
+#include "enums.hh"
+#include "gridselector.hh"
+#include "hdf5-writer.hh"
+#include "hdf5/restart-io.hh"
+#include "matrices.hh"
+#include "program_state.hh"
+#include "one-body-problem-data/bc.hh"
+#include "one-body-problem-data/mybody.hh"
+#include "one-body-problem-data/mygeometry.hh"
+#include "one-body-problem-data/myglobalfrictiondata.hh"
+#include "one-body-problem-data/mygrid.hh"
+#include "one-body-problem-data/weakpatch.hh"
+#include "spatial-solving/solverfactory.hh"
+#include "time-stepping/adaptivetimestepper.hh"
+#include "time-stepping/rate.hh"
+#include "time-stepping/state.hh"
+#include "time-stepping/updaters.hh"
+#include "vtk.hh"
+
+size_t const dims = MY_DIM;
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+    auto const parset = getParameters(argc, argv);
+
+    MyGeometry::render();
+    MyGeometry::write();
+
+    using GridView = Grid::LeafGridView;
+    using MyAssembler = MyAssembler<GridView, dims>;
+    using Matrix = MyAssembler::Matrix;
+    using Vector = MyAssembler::Vector;
+    using LocalVector = Vector::block_type;
+    using ScalarMatrix = MyAssembler::ScalarMatrix;
+    using ScalarVector = MyAssembler::ScalarVector;
+
+    auto const weakPatch =
+        getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
+
+    // {{{ Set up grid
+    GridConstructor<Grid> gridConstructor;
+    auto grid = gridConstructor.getGrid();
+
+    refine(*grid, weakPatch,
+           parset.get<double>("boundary.friction.smallestDiameter"));
+
+    double minDiameter = std::numeric_limits<double>::infinity();
+    double maxDiameter = 0.0;
+    for (auto &&e : elements(grid->leafGridView())) {
+      auto const geometry = e.geometry();
+      auto const diam = diameter(geometry);
+      minDiameter = std::min(minDiameter, diam);
+      maxDiameter = std::max(maxDiameter, diam);
+    }
+    std::cout << "min diameter: " << minDiameter << std::endl;
+    std::cout << "max diameter: " << maxDiameter << std::endl;
+
+    auto const leafView = grid->leafGridView();
+    auto const leafVertexCount = leafView.size(dims);
+
+    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
+
+    auto myFaces = gridConstructor.constructFaces(leafView);
+
+    BoundaryPatch<GridView> const neumannBoundary(leafView);
+    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
+    BoundaryPatch<GridView> const &surface = myFaces.upper;
+
+    // Dirichlet Boundary
+    Dune::BitSetVector<dims> noNodes(leafVertexCount);
+    Dune::BitSetVector<dims> dirichletNodes(leafVertexCount);
+    for (size_t i = 0; i < leafVertexCount; ++i) {
+      if (myFaces.right.containsVertex(i))
+        dirichletNodes[i][0] = true;
+
+      if (myFaces.lower.containsVertex(i))
+        dirichletNodes[i][1] = true;
+
+#if MY_DIM == 3
+      if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i))
+        dirichletNodes[i][2] = true;
+#endif
+    }
+
+    // Set up functions for time-dependent boundary conditions
+    using Function = Dune::VirtualFunction<double, double>;
+    Function const &velocityDirichletFunction = VelocityDirichletCondition();
+    Function const &neumannFunction = NeumannCondition();
+
+    MyAssembler const myAssembler(leafView);
+
+    MyBody<dims> const body(parset);
+
+    Matrices<Matrix> matrices;
+    myAssembler.assembleElasticity(body.getYoungModulus(),
+                                   body.getPoissonRatio(), matrices.elasticity);
+    myAssembler.assembleViscosity(body.getShearViscosityField(),
+                                  body.getBulkViscosityField(),
+                                  matrices.damping);
+    myAssembler.assembleMass(body.getDensityField(), matrices.mass);
+
+    ScalarMatrix relativeFrictionalBoundaryMass;
+    myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
+                                               relativeFrictionalBoundaryMass);
+    relativeFrictionalBoundaryMass /= frictionalBoundary.area();
+    EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
+        relativeFrictionalBoundaryMass);
+
+    // Assemble forces
+    Vector gravityFunctional;
+    myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
+
+    // Problem formulation: right-hand side
+    std::function<void(double, Vector &)> computeExternalForces =
+        [&](double _relativeTime, Vector &_ell) {
+          myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
+                                      _relativeTime);
+          _ell += gravityFunctional;
+        };
+
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(leafVertexCount);
+    auto const firstRestart = parset.get<size_t>("io.restarts.first");
+    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
+    auto const writeRestarts = parset.get<bool>("io.restarts.write");
+    auto const writeData = parset.get<bool>("io.data.write");
+    bool const handleRestarts = writeRestarts or firstRestart > 0;
+
+    auto dataFile =
+        writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+    auto restartFile = handleRestarts
+                           ? std::make_unique<HDF5::File>(
+                                 "restarts.h5",
+                                 writeRestarts ? HDF5::Access::READWRITE
+                                               : HDF5::Access::READONLY)
+                           : nullptr;
+    auto restartIO = handleRestarts
+                         ? std::make_unique<RestartIO<MyProgramState>>(
+                               *restartFile, leafVertexCount)
+                         : nullptr;
+
+    if (firstRestart > 0) // automatically adjusts the time and timestep
+      restartIO->read(firstRestart, programState);
+    else
+      programState.setupInitialConditions(parset, computeExternalForces,
+                                          matrices, myAssembler, dirichletNodes,
+                                          noNodes, frictionalBoundary, body);
+
+    MyGlobalFrictionData<LocalVector> frictionInfo(
+        parset.sub("boundary.friction"), weakPatch);
+    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
+        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
+        frictionalBoundary, frictionInfo, programState.weightedNormalStress);
+    myGlobalFriction->updateAlpha(programState.alpha);
+
+    Vector vertexCoordinates(leafVertexCount);
+    {
+      Dune::MultipleCodimMultipleGeomTypeMapper<
+          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
+      for (auto &&v : vertices(leafView))
+        vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
+    }
+
+    using MyVertexBasis = typename MyAssembler::VertexBasis;
+    auto dataWriter =
+        writeData ? std::make_unique<
+                        HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
+                        *dataFile, vertexCoordinates, myAssembler.vertexBasis,
+                        surface, frictionalBoundary, weakPatch)
+                  : nullptr;
+    MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
+        myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
+    IterationRegister iterationCount;
+    auto const report = [&](bool initial = false) {
+      if (writeData) {
+        dataWriter->reportSolution(programState, *myGlobalFriction);
+        if (!initial)
+          dataWriter->reportIterations(programState, iterationCount);
+        dataFile->flush();
+      }
+
+      if (writeRestarts and !initial and
+          programState.timeStep % restartSpacing == 0) {
+        restartIO->write(programState);
+        restartFile->flush();
+      }
+
+      if (parset.get<bool>("io.printProgress"))
+        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
+                  << ", time = " << std::setw(12) << programState.relativeTime
+                  << ", tau = " << std::setw(12) << programState.relativeTau
+                  << std::endl;
+
+      if (parset.get<bool>("io.vtk.write")) {
+        ScalarVector stress;
+        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+                                           body.getPoissonRatio(),
+                                           programState.u, stress);
+        vtkWriter.write(programState.timeStep, programState.u, programState.v,
+                        programState.alpha, stress);
+      }
+    };
+    report(true);
+
+    // Set up TNNMG solver
+    using NonlinearFactory = SolverFactory<
+        dims,
+        MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+        Grid>;
+    NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
+
+    using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
+                               StateUpdater<ScalarVector, Vector>>;
+    MyUpdater current(
+        initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
+                        velocityDirichletFunction, dirichletNodes, matrices,
+                        programState.u, programState.v, programState.a),
+        initStateUpdater<ScalarVector, Vector>(
+            parset.get<Config::stateModel>("boundary.friction.stateModel"),
+            programState.alpha, *frictionalBoundary.getVertices(),
+            parset.get<double>("boundary.friction.L"),
+            parset.get<double>("boundary.friction.V0")));
+
+    auto const refinementTolerance =
+        parset.get<double>("timeSteps.refinementTolerance");
+    auto const mustRefine = [&](MyUpdater &coarseUpdater,
+                                MyUpdater &fineUpdater) {
+      ScalarVector coarseAlpha;
+      coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+      ScalarVector fineAlpha;
+      fineUpdater.state_->extractAlpha(fineAlpha);
+
+      return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
+    };
+
+    std::signal(SIGXCPU, handleSignal);
+    std::signal(SIGINT, handleSignal);
+    std::signal(SIGTERM, handleSignal);
+
+    AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
+                        EnergyNorm<ScalarMatrix, ScalarVector>>
+        adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
+                            programState.relativeTime, programState.relativeTau,
+                            computeExternalForces, stateEnergyNorm, mustRefine);
+    while (!adaptiveTimeStepper.reachedEnd()) {
+      programState.timeStep++;
+
+      iterationCount = adaptiveTimeStepper.advance();
+
+      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
+      current.rate_->extractDisplacement(programState.u);
+      current.rate_->extractVelocity(programState.v);
+      current.rate_->extractAcceleration(programState.a);
+      current.state_->extractAlpha(programState.alpha);
+
+      report();
+
+      if (terminationRequested) {
+        std::cerr << "Terminating prematurely" << std::endl;
+        break;
+      }
+    }
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  }
+}

src/one-body-sample.parset → src/one-body-problem.cfg

 # -*- mode:conf -*-
-gravity         = 9.81
+gravity         = 9.81  # [m/s^2]
 
 [io]
-printProgress = false
-writeVTK      = false
+data.write      = true
+printProgress   = false
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = true
+vtk.write       = false
 
 [problem]
-finalTime = 15
+finalTime       = 1000  # [s]
 
 [body]
-youngModulus = 5e7
-poissonRatio = 0.3
-density = 5000
-shearViscosity = 0
-bulkViscosity  = 0
+bulkModulus     = 0.5e5 # [Pa]
+poissonRatio    = 0.3   # [1]
+[body.elastic]
+density         = 900   # [kg/m^3]
+shearViscosity  = 1e3   # [Pas]
+bulkViscosity   = 1e3   # [Pas]
+[body.viscoelastic]
+density         = 1000  # [kg/m^3]
+shearViscosity  = 1e4   # [Pas]
+bulkViscosity   = 1e4   # [Pas]
 
 [boundary.friction]
-C            = 0.0
-mu0          = 0.6
-a            = 0.010
-b            = 0.015
-V0           = 1e-6
-L            = 1e-5
-initialState = 4.54e-05 # = exp(-10) = theta, so that alpha = -10
-stateModel   = Dieterich
+C               = 10    # [Pa]
+mu0             = 0.7   # [ ]
+V0              = 5e-5  # [m/s]
+L               = 2.25e-5 # [m]
+initialAlpha    = 0     # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Regularised
+[boundary.friction.weakening]
+a               = 0.002 # [ ]
+b               = 0.017 # [ ]
+[boundary.friction.strengthening]
+a               = 0.020 # [ ]
+b               = 0.005 # [ ]
 
 [timeSteps]
-number = 10000
 scheme = newmark
 
-[grid]
-refinements = 4
-
 [u0.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [a0.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.fpi]
-tolerance         = 1e-10
 maximumIterations = 10000
-relaxation        = 0.5
-requiredReduction = 0.5
+lambda            = 0.5
 
 [solver.tnnmg.linear]
 maxiumumIterations = 100000
-tolerance          = 1e-10
 pre                = 3
 cycle              = 1  # 1 = V, 2 = W, etc.
 post               = 3
@@ -64,6 +68,3 @@ post               = 3
 pre   = 1
 multi = 5 # number of multigrid steps
 post  = 0
-
-[localsolver]
-steps = 1

src/one-body-sample.cc

-#include <Python.h>
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifndef HAVE_PYTHON
-#error Python is required
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#ifndef srcdir
-#error srcdir unset
-#endif
-
-#ifndef DIM
-#error DIM unset
-#endif
-
-#if !HAVE_ALUGRID
-#error ALUGRID is required
-#endif
-
-#include <cmath>
-#include <exception>
-#include <fstream>
-#include <iostream>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/exceptions.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
-#include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/utility/structuredgridfactory.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wdeprecated-declarations"
-#include <dune/fufem/boundarypatch.hh>
-#pragma clang diagnostic pop
-
-#include <dune/fufem/dunepython.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/sharedpointermap.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/norms/sumnorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/globalnonlinearity.hh>
-
-#include "assemblers.hh"
-#include "enum_parser.cc"
-#include "enum_scheme.cc"
-#include "enum_state_model.cc"
-#include "enum_verbosity.cc"
-#include "enums.hh"
-#include "friction_writer.hh"
-#include "mybody.hh"
-#include "mygeometry.hh"
-#include "myglobalfrictiondata.hh"
-#include "mygrid.hh"
-#include "solverfactory.hh"
-#include "state.hh"
-#include "timestepping.hh"
-#include "vtk.hh"
-
-size_t const dims = DIM;
-
-void initPython() {
-  Python::start();
-
-  Python::run("import sys");
-  Python::run("sys.path.append('" srcdir "')");
-}
-
-int main(int argc, char *argv[]) {
-  try {
-    Dune::ParameterTree parset;
-    Dune::ParameterTreeParser::readINITree(srcdir "/one-body-sample.parset",
-                                           parset);
-    Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-
-    MyGeometry const myGeometry;
-
-    // {{{ Set up grid
-    using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
-    auto grid = constructGrid<Grid>(myGeometry); // FIXME
-
-    auto const refinements = parset.get<size_t>("grid.refinements");
-    grid->globalRefine(refinements);
-    size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
-
-    using GridView = Grid::LeafGridView;
-    GridView const leafView = grid->leafView();
-    // }}}
-
-    // Set up myFaces
-    MyFaces<GridView, MyGeometry> myFaces(leafView, myGeometry);
-
-    // Neumann boundary
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
-
-    // Frictional Boundary
-    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
-    frictionalBoundary.getVertices(frictionalNodes);
-
-    // Dirichlet Boundary
-    Dune::BitSetVector<dims> noNodes(fineVertexCount);
-    Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
-    for (size_t i = 0; i < fineVertexCount; ++i) {
-      if (myFaces.lower.containsVertex(i))
-        dirichletNodes[i][1] = true;
-
-      if (myFaces.upper.containsVertex(i))
-        dirichletNodes[i] = true;
-    }
-
-    // Set up functions for time-dependent boundary conditions
-    using Function = Dune::VirtualFunction<double, double>;
-    using FunctionMap = SharedPointerMap<std::string, Function>;
-    FunctionMap functions;
-    {
-      initPython();
-      Python::import("one-body-sample")
-          .get("Functions")
-          .toC<typename FunctionMap::Base>(functions);
-    }
-    auto const &velocityDirichletFunction =
-                   functions.get("velocityDirichletCondition"),
-               &neumannFunction = functions.get("neumannCondition");
-
-    using MyAssembler = MyAssembler<GridView, dims>;
-    using Matrix = MyAssembler::Matrix;
-    using LocalMatrix = Matrix::block_type;
-    using Vector = MyAssembler::Vector;
-    using LocalVector = Vector::block_type;
-    using ScalarMatrix = MyAssembler::ScalarMatrix;
-    using ScalarVector = MyAssembler::ScalarVector;
-
-    MyAssembler myAssembler(leafView);
-
-    MyBody<dims> const body(parset, myGeometry);
-
-    Matrix A, C, M;
-    myAssembler.assembleElasticity(body.getYoungModulus(),
-                                   body.getPoissonRatio(), A);
-    myAssembler.assembleViscosity(body.getShearViscosity(),
-                                  body.getBulkViscosity(), C);
-    myAssembler.assembleMass(body.getDensityField(), M);
-    EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
-    // Q: Does it make sense to weigh them in this manner?
-    SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
-
-    ScalarMatrix frictionalBoundaryMass;
-    myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
-                                               frictionalBoundaryMass);
-    EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
-        frictionalBoundaryMass);
-
-    // Assemble forces
-    Vector gravityFunctional;
-    myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
-
-    // Problem formulation: right-hand side
-    auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
-      myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
-                                  _relativeTime);
-      _ell += gravityFunctional;
-    };
-    Vector ell(fineVertexCount);
-    computeExternalForces(0.0, ell);
-
-    ScalarVector normalStress(fineVertexCount);
-    normalStress = (myGeometry.A[1] - myGeometry.C[1]) *
-                   parset.get<double>("gravity") *
-                   parset.get<double>("body.density");
-    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"),
-                                            myGeometry);
-
-    // {{{ Initial conditions
-    ScalarVector alpha_initial(fineVertexCount);
-    alpha_initial =
-        std::log(parset.get<double>("boundary.friction.initialState"));
-    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionInfo, normalStress);
-    myGlobalNonlinearity->updateLogState(alpha_initial);
-
-    using LinearFactory = SolverFactory<
-        dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                  ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-        Grid>;
-    ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
-    // TODO: clean up once generic lambdas arrive
-    auto const solveLinearProblem = [&](
-        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-        Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
-        Dune::ParameterTree const &_localParset) {
-      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
-                            refinements, *grid, _dirichletNodes);
-
-      typename LinearFactory::ConvexProblem convexProblem(
-          1.0, _matrix, zeroNonlinearity, _rhs, _x);
-      typename LinearFactory::BlockProblem problem(parset, convexProblem);
-
-      auto multigridStep = factory.getSolver();
-      multigridStep->setProblem(_x, problem);
-      LoopSolver<Vector> solver(
-          multigridStep, _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), &_norm,
-          _localParset.get<Solver::VerbosityMode>("verbosity"),
-          false); // absolute error
-
-      solver.preprocess();
-      solver.solve();
-    };
-
-    // Solve the stationary problem
-    Vector u_initial(fineVertexCount);
-    u_initial = 0.0;
-    solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
-                       parset.sub("u0.solver"));
-    Vector v_initial(fineVertexCount);
-    v_initial = 0.0;
-    {
-      // Prescribe a homogeneous velocity field in the x-direction
-      // This way, the initial condition for v and the Dirichlet
-      // condition match up at t=0
-      double v_initial_const;
-      velocityDirichletFunction.evaluate(0.0, v_initial_const);
-      for (auto &x : v_initial)
-        x[0] = v_initial_const;
-    }
-    Vector a_initial(fineVertexCount);
-    a_initial = 0.0;
-    {
-      // We solve Ma = ell - [Au + Cv + Psi(v)]
-      Vector accelerationRHS(fineVertexCount);
-      {
-        accelerationRHS = 0.0;
-        Arithmetic::addProduct(accelerationRHS, A, u_initial);
-        Arithmetic::addProduct(accelerationRHS, C, v_initial);
-        // NOTE: We assume differentiability of Psi at v0 here!
-        myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
-        accelerationRHS *= -1.0;
-        accelerationRHS += ell;
-      }
-      solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
-                         parset.sub("a0.solver"));
-    }
-    // }}}
-
-    Vector vertexCoordinates(fineVertexCount);
-    {
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
-      for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
-        auto const geometry = it->geometry();
-        assert(geometry.corners() == 1);
-        vertexCoordinates[vertexMapper.map(*it)] = geometry.corner(0);
-      }
-    }
-    FrictionWriter<ScalarVector, Vector> writer(vertexCoordinates,
-                                                frictionalNodes);
-    auto const reportFriction = [&](Vector const &_u, Vector const &_v,
-                                    ScalarVector const &_alpha) {
-      ScalarVector c;
-      myGlobalNonlinearity->coefficientOfFriction(_v, c);
-      writer.writeKinetics(_u, _v);
-      writer.writeOther(c, _alpha);
-    };
-    reportFriction(u_initial, v_initial, alpha_initial);
-
-    MyVTKWriter<typename MyAssembler::VertexBasis,
-                typename MyAssembler::CellBasis> const
-    vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
-    // Set up TNNMG solver
-    using NonlinearFactory =
-        SolverFactory<dims, MyBlockProblem<ConvexProblem<
-                                GlobalNonlinearity<Matrix, Vector>, Matrix>>,
-                      Grid>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
-                             dirichletNodes);
-    auto multigridStep = factory.getSolver();
-
-    std::fstream iterationWriter("iterations", std::fstream::out),
-        relaxationWriter("relaxation", std::fstream::out);
-
-    auto timeSteppingScheme =
-        initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
-                        velocityDirichletFunction, dirichletNodes, M, A, C,
-                        u_initial, v_initial, a_initial);
-    auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
-        parset.get<Config::stateModel>("boundary.friction.stateModel"),
-        alpha_initial, frictionalNodes,
-        parset.get<double>("boundary.friction.L"));
-
-    Vector v = v_initial;
-    ScalarVector alpha(fineVertexCount);
-
-    auto const timeSteps = parset.get<size_t>("timeSteps.number"),
-               maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
-               maximumIterations =
-                   parset.get<size_t>("v.solver.maximumIterations");
-    auto const tau = parset.get<double>("problem.finalTime") / timeSteps,
-               tolerance = parset.get<double>("v.solver.tolerance"),
-               fixedPointTolerance = parset.get<double>("v.fpi.tolerance"),
-               relaxation = parset.get<double>("v.fpi.relaxation"),
-               requiredReduction =
-                   parset.get<double>("v.fpi.requiredReduction");
-    auto const printProgress = parset.get<bool>("io.printProgress");
-    auto const verbosity =
-        parset.get<Solver::VerbosityMode>("v.solver.verbosity");
-    for (size_t timeStep = 1; timeStep <= timeSteps; ++timeStep) {
-      if (printProgress)
-        std::cout << std::setw(7) << timeStep << " " << std::flush;
-
-      stateUpdater->nextTimeStep();
-      timeSteppingScheme->nextTimeStep();
-
-      auto const relativeTime = double(timeStep) / double(timeSteps);
-      computeExternalForces(relativeTime, ell);
-
-      Matrix velocityMatrix;
-      Vector velocityRHS(fineVertexCount);
-      Vector velocityIterate(fineVertexCount);
-
-      stateUpdater->setup(tau);
-      timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,
-                                velocityIterate, velocityMatrix);
-
-      LoopSolver<Vector> velocityProblemSolver(multigridStep, maximumIterations,
-                                               tolerance, &AMNorm, verbosity,
-                                               false); // absolute error
-
-      size_t iterationCounter;
-      auto solveVelocityProblem = [&](Vector &_velocityIterate,
-                                      ScalarVector const &_alpha) {
-        myGlobalNonlinearity->updateLogState(_alpha);
-
-        // NIT: Do we really need to pass u here?
-        typename NonlinearFactory::ConvexProblem const convexProblem(
-            1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
-            _velocityIterate);
-        typename NonlinearFactory::BlockProblem velocityProblem(parset,
-                                                                convexProblem);
-        multigridStep->setProblem(_velocityIterate, velocityProblem);
-
-        velocityProblemSolver.preprocess();
-        velocityProblemSolver.solve();
-        iterationCounter = velocityProblemSolver.getResult().iterations;
-      };
-
-      // Since the velocity explodes in the quasistatic case, use the
-      // displacement as a convergence criterion
-      // Q: is this reasonable?
-      Vector u;
-      Vector u_saved;
-      ScalarVector alpha_saved;
-      double lastStateCorrection;
-      for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
-        stateUpdater->solve(v);
-        stateUpdater->extractLogState(alpha);
-
-        if (stateFPI == 1)
-          relaxationWriter << "N ";
-        else {
-          double const stateCorrection =
-              stateEnergyNorm.diff(alpha, alpha_saved);
-          if (stateFPI <= 2 // lastStateCorrection is only set for stateFPI > 2
-              or stateCorrection < requiredReduction * lastStateCorrection)
-            relaxationWriter << "N ";
-          else {
-            alpha *= (1.0 - relaxation);
-            Arithmetic::addProduct(alpha, relaxation, alpha_saved);
-            relaxationWriter << "Y ";
-          }
-          lastStateCorrection = stateCorrection;
-        }
-
-        solveVelocityProblem(velocityIterate, alpha);
-        timeSteppingScheme->postProcess(velocityIterate);
-        timeSteppingScheme->extractDisplacement(u);
-        timeSteppingScheme->extractVelocity(v);
-
-        iterationWriter << iterationCounter << " ";
-        if (printProgress)
-          std::cout << '.' << std::flush;
-
-        if (stateFPI > 1) {
-          double const velocityCorrection = AMNorm.diff(u_saved, u);
-          if (velocityCorrection < fixedPointTolerance)
-            break;
-        }
-        if (stateFPI == maximumStateFPI)
-          DUNE_THROW(Dune::Exception, "FPI failed to converge");
-
-        alpha_saved = alpha;
-        u_saved = u;
-      }
-      if (printProgress)
-        std::cout << std::endl;
-
-      reportFriction(u, v, alpha);
-      iterationWriter << std::endl;
-      relaxationWriter << std::endl;
-
-      if (parset.get<bool>("io.writeVTK")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(), u, stress);
-        vtkWriter.write(timeStep - 1, u, v, alpha, stress);
-      }
-    }
-    iterationWriter.close();
-    relaxationWriter.close();
-
-    Python::stop();
-  }
-  catch (Dune::Exception &e) {
-    Dune::derr << "Dune reported error: " << e << std::endl;
-  }
-  catch (std::exception &e) {
-    std::cerr << "Standard exception: " << e.what() << std::endl;
-  }
-}

src/one-body-sample.py

-class neumannCondition:
-    def __call__(self, relativeTime):
-        return 0
-
-class velocityDirichletCondition:
-    def __call__(self, relativeTime):
-        if relativeTime <= 0.25:
-            factor = 4*relativeTime;
-        else:
-            factor = 1
-        return 2e-4 * factor
-
-Functions = {
-    'neumannCondition' : neumannCondition(),
-    'velocityDirichletCondition' : velocityDirichletCondition()
-}

src/program_state.hh

+#ifndef SRC_PROGRAM_STATE_HH
+#define SRC_PROGRAM_STATE_HH
+
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
+#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
+
+#include <dune/tectonic/body.hh>
+
+#include "assemblers.hh"
+#include "matrices.hh"
+#include "spatial-solving/solverfactory.hh"
+
+template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
+public:
+  using Vector = VectorTEMPLATE;
+  using ScalarVector = ScalarVectorTEMPLATE;
+
+  ProgramState(int leafVertexCount)
+      : u(leafVertexCount),
+        v(leafVertexCount),
+        a(leafVertexCount),
+        alpha(leafVertexCount),
+        weightedNormalStress(leafVertexCount) {}
+
+  // Set up initial conditions
+  template <class Matrix, class GridView>
+  void setupInitialConditions(
+      Dune::ParameterTree const &parset,
+      std::function<void(double, Vector &)> externalForces,
+      Matrices<Matrix> const matrices,
+      MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
+      Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
+      Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      Body<Vector::block_type::dimension> const &body) {
+
+    using LocalVector = typename Vector::block_type;
+    using LocalMatrix = typename Matrix::block_type;
+    auto constexpr dims = LocalVector::dimension;
+
+    // Solving a linear problem with a multigrid solver
+    auto const solveLinearProblem = [&](
+        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
+        Vector const &_rhs, Vector &_x,
+        Dune::ParameterTree const &_localParset) {
+
+      using LinearFactory = SolverFactory<
+          dims, BlockNonlinearTNNMGProblem<ConvexProblem<
+                    ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+          typename GridView::Grid>;
+      ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
+
+      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
+                            myAssembler.gridView.grid(), _dirichletNodes);
+
+      typename LinearFactory::ConvexProblem convexProblem(
+          1.0, _matrix, zeroNonlinearity, _rhs, _x);
+      typename LinearFactory::BlockProblem problem(parset, convexProblem);
+
+      auto multigridStep = factory.getStep();
+      multigridStep->setProblem(_x, problem);
+      EnergyNorm<Matrix, Vector> const norm(_matrix);
+      LoopSolver<Vector> solver(
+          multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
+          _localParset.get<double>("tolerance"), &norm,
+          _localParset.get<Solver::VerbosityMode>("verbosity"),
+          false); // absolute error
+
+      solver.preprocess();
+      solver.solve();
+    };
+
+    timeStep = 0;
+    relativeTime = 0.0;
+    relativeTau = 1e-6;
+
+    Vector ell0(u.size());
+    externalForces(relativeTime, ell0);
+
+    // Initial velocity
+    v = 0.0;
+
+    // Initial displacement: Start from a situation of minimal stress,
+    // which is automatically attained in the case [v = 0 = a].
+    // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
+    solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u,
+                       parset.sub("u0.solver"));
+
+    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+    Vector accelerationRHS = ell0;
+    Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
+    solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
+                       parset.sub("a0.solver"));
+
+    // Initial state
+    alpha = parset.get<double>("boundary.friction.initialAlpha");
+
+    // Initial normal stress
+    myAssembler.assembleWeightedNormalStress(
+        frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
+        body.getPoissonRatio(), u);
+  }
+
+public:
+  Vector u;
+  Vector v;
+  Vector a;
+  ScalarVector alpha;
+  ScalarVector weightedNormalStress;
+  double relativeTime;
+  double relativeTau;
+  size_t timeStep;
+};
+
+#endif

src/spatial-solving/fixedpointiterator.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/common/exceptions.hh>
+
+#include <dune/solvers/common/arithmetic.hh>
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../enums.hh"
+#include "../enumparser.hh"
+
+#include "fixedpointiterator.hh"
+
+void FixedPointIterationCounter::operator+=(
+    FixedPointIterationCounter const &other) {
+  iterations += other.iterations;
+  multigridIterations += other.multigridIterations;
+}
+
+template <class Factory, class Updaters, class ErrorNorm>
+FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
+    Factory &factory, Dune::ParameterTree const &parset,
+    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
+    : step_(factory.getStep()),
+      parset_(parset),
+      globalFriction_(globalFriction),
+      fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
+      fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
+      lambda_(parset.get<double>("v.fpi.lambda")),
+      velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
+      velocityTolerance_(parset.get<double>("v.solver.tolerance")),
+      verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
+      errorNorm_(errorNorm) {}
+
+template <class Factory, class Updaters, class ErrorNorm>
+FixedPointIterationCounter
+FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
+    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
+    Vector &velocityIterate) {
+  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
+  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
+                                           velocityTolerance_, &energyNorm,
+                                           verbosity_, false); // absolute error
+
+  size_t fixedPointIteration;
+  size_t multigridIterations = 0;
+  ScalarVector alpha;
+  updaters.state_->extractAlpha(alpha);
+  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
+       ++fixedPointIteration) {
+    // solve a velocity problem
+    globalFriction_->updateAlpha(alpha);
+    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
+                                velocityRHS, velocityIterate);
+    BlockProblem velocityProblem(parset_, convexProblem);
+    step_->setProblem(velocityIterate, velocityProblem);
+    velocityProblemSolver.preprocess();
+    velocityProblemSolver.solve();
+
+    multigridIterations += velocityProblemSolver.getResult().iterations;
+
+    Vector v_m;
+    updaters.rate_->extractOldVelocity(v_m);
+    v_m *= 1.0 - lambda_;
+    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
+
+    // solve a state problem
+    updaters.state_->solve(v_m);
+    ScalarVector newAlpha;
+    updaters.state_->extractAlpha(newAlpha);
+
+    if (lambda_ < 1e-12 or
+        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
+      fixedPointIteration++;
+      break;
+    }
+    alpha = newAlpha;
+  }
+  if (fixedPointIteration == fixedPointMaxIterations_)
+    DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+  updaters.rate_->postProcess(velocityIterate);
+
+  // Cannot use return { fixedPointIteration, multigridIterations };
+  // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
+  FixedPointIterationCounter ret;
+  ret.iterations = fixedPointIteration;
+  ret.multigridIterations = multigridIterations;
+  return ret;
+}
+
+std::ostream &operator<<(std::ostream &stream,
+                         FixedPointIterationCounter const &fpic) {
+  return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
+                << ")";
+}
+
+#include "fixedpointiterator_tmpl.cc"