[Problem] Add a new 2D sand wedge geometry

.gitignore

@@ -26,4 +26,5 @@
 /test-driver
 Makefile
 Makefile.in
+src/sand-wedge
 src/sliding-block-?D

configure.ac

@@ -9,6 +9,16 @@ AC_CONFIG_HEADERS([config.h])
 
 DUNE_CHECK_ALL
 
+AC_ARG_WITH([cairomm], AS_HELP_STRING([--with-cairomm], [Use cairomm]))
+
+AS_IF([test "x$with_cairomm" != "xno"], [
+   PKG_CHECK_MODULES([CAIROMM], [cairomm-1.0],
+                     [AM_CONDITIONAL(CAIROMM, true)
+                      AC_DEFINE([HAVE_CAIROMM], [1], [If CairoMM is available])
+   ])
+])
+
+
 AC_CONFIG_FILES([
   Makefile
   src/Makefile

m4/dune-tectonic.m4

@@ -6,7 +6,9 @@ dnl -*- autoconf -*-
 # 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])
+AC_DEFUN([DUNE_TECTONIC_CHECKS],[
+  AC_REQUIRE([AX_BOOST_BASE])
+])
 
 # Additional checks needed to find dune-tectonic
 # This macro should be invoked by every module which depends on dune-tectonic, but

src/Makefile.am

-bin_PROGRAMS = sliding-block-2D
+bin_PROGRAMS = sand-wedge sliding-block-2D
 
 common_sources = \
 	assemblers.cc \
@@ -8,6 +8,10 @@ common_sources = \
 	timestepping.cc \
 	vtk.cc
 
+sand_wedge_SOURCES = $(common_sources) sand-wedge.cc
+sand_wedge_CPPFLAGS = \
+	$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
+	-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=2
 sliding_block_2D_SOURCES = $(common_sources) sliding-block.cc
 sliding_block_2D_CPPFLAGS = \
 	$(AM_CPPFLAGS) \
@@ -41,4 +45,9 @@ AM_LDFLAGS = \
 	$(ALUGRID_LDFLAGS) \
 	$(PYTHON_LDFLAGS)
 
+if CAIROMM
+AM_CPPFLAGS += $(CAIROMM_CFLAGS)
+LDADD += $(CAIROMM_LIBS)
+endif
+
 include $(top_srcdir)/am/global-rules

src/sand-wedge-data/boundaryconditions.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 -5e-5 * factor
+
+Functions = {
+    'neumannCondition' : neumannCondition(),
+    'velocityDirichletCondition' : velocityDirichletCondition()
+}

src/sand-wedge-data/geometry.tex

+\documentclass[tikz]{minimal}
+
+\usepackage{tikz}
+\usetikzlibrary{calc}
+\usetikzlibrary{decorations.pathreplacing}
+
+
+\usepackage{siunitx}
+
+\begin{document}
+\begin{tikzpicture}[scale=12]
+  \pgfmathsetmacro{\hypot}{sqrt(1+.27^2)};
+  \pgfmathsetmacro{\ssin}{sin(15)};
+  \pgfmathsetmacro{\scos}{cos(15)};
+  \pgfmathsetmacro{\lsin}{sin(75)};
+
+  \coordinate (A) at (0,0);
+  \node at (A) [left] {A};
+  \coordinate (B) at (\hypot,0);
+  \node at (B) [right] {B};
+  \coordinate (C) at (\lsin, .27*\lsin);
+  \node at (C) [above] {C};
+
+  \draw (A)
+  -- (B)
+  -- (C)
+  -- node [above=.5cm, sloped] {$\overline{AC}=\SI{100}{cm}$} (A);
+
+  \coordinate (Y) at ($(0.35/\lsin, 0)$);
+  \node at (Y) [below right] {Y};
+  \coordinate (Z) at ($(0.35*\lsin,0.35*\ssin)$);
+  \node at (Z) [above] {Z};
+  \coordinate (X) at ($(Y) - (0.20,0)$);
+  \node at (X) [below left] {X};
+  \path let \p1 = (X), \p2 = (Y), \p3 = (Z) in
+    coordinate (U) at ($(\x1 * \x3 / \x2, \x1 * \y3 / \x2)$);
+  \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, sloped] {$\overline{XY}=\SI{20}{cm}$} (Y);
+  \draw[dashed] (Y) -- (Z);
+  \draw[dashed] (U) -- (X);
+
+  \coordinate (K) at ($(B) - (0.06 / \ssin,0)$);
+  \node at (K) [below] {K};
+  \coordinate (M) at ($(B) + (-0.06 * \ssin, 0.06 * \lsin)$);
+  \node at (M) [right] {M};
+
+  \path (C) -- node [right=.5cm] {$\overline{CM} = \SI{21}{cm}$} (M);
+
+  \draw[thick] (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};
+
+  % Note: Lengths and position are arbitrarily chosen here
+  \draw[->] (0.4,0.45) -- node [above, sloped] {vertical} +(-0.15*\ssin,0.15*\scos);
+  \draw[->] (0.4,0.45) -- node [above, sloped] {horizontal} +(-0.15*\scos,-0.15*\ssin);
+
+  \node[align=left] at (0.5,-0.125) {
+    $Z$: coast line\\
+    $\overline{XY}$: velocity weakening zone\\
+    $BKM$: visco-elastic domain};
+\end{tikzpicture}
+
+\end{document}

src/sand-wedge-data/mybody.hh

+#ifndef SRC_SAND_WEDGE_DATA_MYBODY_HH
+#define SRC_SAND_WEDGE_DATA_MYBODY_HH
+
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/functions/constantfunction.hh>
+
+#include <dune/tectonic/body.hh>
+#include <dune/tectonic/gravity.hh>
+
+#include "twopiece.hh"
+#include "mygeometry.hh"
+
+template <int dimension> class MyBody : public Body<dimension> {
+  using typename Body<dimension>::ScalarFunction;
+  using typename Body<dimension>::VectorField;
+
+public:
+  MyBody(Dune::ParameterTree const &parset)
+      : poissonRatio_(parset.get<double>("body.poissonRatio")),
+        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_; }
+  ScalarFunction const &getShearViscosityField() const override {
+    return shearViscosityField_;
+  }
+  ScalarFunction const &getBulkViscosityField() const override {
+    return bulkViscosityField_;
+  }
+  ScalarFunction const &getDensityField() const override {
+    return densityField_;
+  }
+  VectorField const &getGravityField() const override { return gravityField_; }
+
+private:
+  double const poissonRatio_;
+  double const youngModulus_;
+  TwoPieceFunction const shearViscosityField_;
+  TwoPieceFunction const bulkViscosityField_;
+  TwoPieceFunction const densityField_;
+  Gravity<dimension> const gravityField_;
+};
+#endif

src/sand-wedge-data/mygeometry.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#include "mygeometry.hh"
+
+double MyGeometry::verticalProjection(LocalVector const &x) {
+  return zenith * x;
+}
+double MyGeometry::horizontalProjection(LocalVector const &x) {
+  return rotatedZenith * x;
+}
+
+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 = [&](LocalVector const &v) { cr->move_to(v[0], -v[1]); };
+  auto const lineTo = [&](LocalVector const &v) { cr->line_to(v[0], -v[1]); };
+
+  cr->scale(widthScale, heightScale);
+  cr->translate(0.1, 0.1);
+  cr->rotate(leftAngle);
+  cr->set_line_width(0.0025);
+
+  // triangle
+  {
+    moveTo(A);
+    lineTo(B);
+    lineTo(C);
+    cr->close_path();
+    cr->stroke();
+  }
+
+  // dashed lines
+  {
+    cr->save();
+    std::vector<double> dashPattern = { 0.005 };
+    cr->set_dash(dashPattern, 0);
+    moveTo(Z);
+    lineTo(Y);
+    moveTo(U);
+    lineTo(X);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // fill viscoelastic region
+  {
+    cr->save();
+    setRGBColor(0x0097E0);
+    moveTo(B);
+    lineTo(K);
+    lineTo(M);
+    cr->fill();
+    cr->restore();
+  }
+
+  // mark weakening region
+  {
+    cr->save();
+    setRGBColor(0x7AD3FF);
+    cr->set_line_width(0.005);
+    moveTo(X);
+    lineTo(Y);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // mark points
+  {
+    auto const drawCircle = [&](LocalVector 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(A);
+    drawCircle(B);
+    drawCircle(C);
+    drawCircle(Y);
+    drawCircle(X);
+    drawCircle(Z);
+    drawCircle(U);
+    drawCircle(K);
+    drawCircle(M);
+    drawCircle(G);
+    drawCircle(H);
+    drawCircle(J);
+    drawCircle(I);
+    cr->restore();
+  }
+
+  // labels
+  {
+    auto const label = [&](LocalVector const &v,
+                           std::string l) { // TODO: implicit?
+      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(A, "A");
+    label(B, "B");
+    label(C, "C");
+    label(K, "K");
+    label(M, "M");
+    label(U, "U");
+    label(X, "X");
+    label(Y, "Y");
+    label(Z, "Z");
+    label(G, "G");
+    label(H, "H");
+    label(J, "J");
+    label(I, "I");
+    cr->restore();
+  }
+
+  surface->write_to_png(filename);
+#endif
+}

src/sand-wedge-data/mygeometry.hh

+#ifndef MY_GEOMETRY_HH
+#define MY_GEOMETRY_HH
+
+#include <dune/common/fassign.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/fufem/arithmetic.hh>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+namespace MyGeometry {
+namespace {
+  using LocalVector = Dune::FieldVector<double, 2>;
+
+  // kludge because fieldvectors have no initialiser_list constructor, see
+  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
+  LocalVector generateVector(double x, double y) {
+    LocalVector tmp;
+    tmp <<= x, y;
+    return tmp;
+  }
+
+  double degreesToRadians(double degrees) {
+    return M_PI * degrees / 180;
+  };
+
+  LocalVector convexCombination(LocalVector const &x, LocalVector const &y) {
+    LocalVector ret;
+    ret = 0.0;
+    Arithmetic::addProduct(ret, 0.5, x);
+    Arithmetic::addProduct(ret, 0.5, y);
+    return ret;
+  }
+
+  double leftLeg = 1.00;
+  double rightLeg = 0.27;
+  double leftAngle = degreesToRadians(15);
+  double rightAngle = degreesToRadians(75);
+
+  double const depth = 0.10;
+}
+
+LocalVector const A = generateVector(0, 0);
+LocalVector const B = generateVector(std::hypot(leftLeg, rightLeg), 0);
+LocalVector const C = generateVector(leftLeg * std::sin(rightAngle),
+                                     leftLeg * std::sin(leftAngle));
+
+LocalVector const Y = generateVector(0.35 / std::sin(rightAngle), 0);
+LocalVector const X = generateVector(Y[0] - 0.20, 0);
+LocalVector const Z =
+    generateVector(0.35 * std::sin(rightAngle), 0.35 * std::sin(leftAngle));
+
+LocalVector const U = generateVector(Z[0] * X[0] / Y[0], Z[1] * X[0] / Y[0]);
+
+LocalVector const K = generateVector(B[0] - 0.06 / std::sin(leftAngle), 0);
+LocalVector const M = generateVector(B[0] - 0.06 * std::sin(leftAngle),
+                                     0.06 * std::sin(rightAngle));
+
+LocalVector const G = convexCombination(A, X);
+LocalVector const H = convexCombination(X, Y);
+LocalVector const J = convexCombination(Y, B);
+
+LocalVector const I = generateVector(Y[0] + G[0], 0);
+
+LocalVector const zenith =
+    generateVector(-std::sin(leftAngle), std::cos(leftAngle));
+
+LocalVector const rotatedZenith =
+    generateVector(-std::cos(leftAngle), -std::sin(leftAngle));
+
+double verticalProjection(LocalVector const &);
+double horizontalProjection(LocalVector const &);
+
+void write();
+
+void render();
+}
+#endif

src/sand-wedge-data/myglobalfrictiondata.hh

+#ifndef SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "mygeometry.hh"
+#include "patchfunction.hh"
+
+template <int dimension>
+class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+private:
+  using typename GlobalFrictionData<dimension>::VirtualFunction;
+
+public:
+  MyGlobalFrictionData(Dune::ParameterTree const &parset)
+      : 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")),
+        b_(parset.get<double>("strengthening.b"),
+           parset.get<double>("weakening.b")),
+        mu0_(parset.get<double>("mu0")) {}
+
+  double virtual const &C() const override { return C_; }
+  double virtual const &L() const override { return L_; }
+  double virtual const &V0() const override { return V0_; }
+  VirtualFunction virtual const &a() const override { return a_; }
+  VirtualFunction virtual const &b() const override { return b_; }
+  double virtual 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/sand-wedge-data/mygrid.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "mygrid.hh"
+
+template <class Grid> std::shared_ptr<Grid> constructGrid() {
+  Dune::GridFactory<Grid> gridFactory;
+
+  Dune::FieldMatrix<double, 3, DIM> vertices;
+  vertices[0] = MyGeometry::A;
+  vertices[1] = MyGeometry::B;
+  vertices[2] = MyGeometry::C;
+  for (size_t i = 0; i < vertices.N(); ++i)
+    gridFactory.insertVertex(vertices[i]);
+
+  Dune::GeometryType cell;
+  cell.makeTriangle();
+
+  std::vector<std::vector<unsigned int>> simplices = { { 0, 1, 2 } };
+
+  // sanity-check choices of simplices
+  for (size_t i = 0; i < simplices.size(); ++i) {
+    Dune::FieldMatrix<double, DIM, DIM> check;
+    for (size_t j = 0; j < DIM; ++j)
+      check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+    assert(check.determinant() > 0);
+    gridFactory.insertElement(cell, simplices[i]);
+  }
+  return std::shared_ptr<Grid>(gridFactory.createGrid());
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView)
+    : lower(gridView), right(gridView) {
+  using Grid = typename GridView::Grid;
+  using LevelGridView = typename Grid::LevelGridView;
+  LevelGridView const rootView = gridView.grid().levelView(0);
+
+  lower.insertFacesByProperty([&](
+      typename Grid::LeafGridView::Intersection const &in) {
+    return isClose(0, in.geometry().center()[1]);
+  });
+
+  BoundaryPatch<LevelGridView> upperRoot(rootView);
+  upperRoot.insertFacesByProperty([&](
+      typename LevelGridView::Intersection const &in) {
+    auto const geometry = in.geometry();
+    auto const range = boost::irange(0, geometry.corners());
+    return boost::algorithm::all_of(range, [&](int i) {
+      auto const &co = geometry.corner(i);
+      return (isClose(co[1], MyGeometry::A[1]) &&
+              isClose(co[0], MyGeometry::A[0])) ||
+             (isClose(co[1], MyGeometry::C[1]) &&
+              isClose(co[0], MyGeometry::C[0]));
+    });
+  });
+  BoundaryPatchProlongator<Grid>::prolong(upperRoot, upper);
+
+  BoundaryPatch<LevelGridView> rightRoot(rootView);
+  rightRoot.insertFacesByProperty([&](
+      typename LevelGridView::Intersection const &in) {
+    auto const geometry = in.geometry();
+    auto const range = boost::irange(0, geometry.corners());
+    return boost::algorithm::all_of(range, [&](int i) {
+      auto const &co = geometry.corner(i);
+      return (isClose(co[1], MyGeometry::C[1]) &&
+              isClose(co[0], MyGeometry::C[0])) ||
+             (isClose(co[1], MyGeometry::B[1]) &&
+              isClose(co[0], MyGeometry::B[0]));
+    });
+  });
+  BoundaryPatchProlongator<Grid>::prolong(rightRoot, right);
+}
+
+#include "mygrid_tmpl.cc"

src/sand-wedge-data/mygrid.hh

+#ifndef SRC_SAND_WEDGE_DATA_MYGRID_HH
+#define SRC_SAND_WEDGE_DATA_MYGRID_HH
+
+#include <boost/range/irange.hpp>
+#include <boost/algorithm/cxx11/all_of.hpp>
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wdeprecated-declarations"
+#include <dune/fufem/boundarypatch.hh>
+#pragma clang diagnostic pop
+#include <dune/fufem/boundarypatchprolongator.hh>
+
+#include "mygeometry.hh"
+
+template <class Grid> std::shared_ptr<Grid> constructGrid();
+
+template <class GridView> struct MyFaces {
+  BoundaryPatch<GridView> lower;
+  BoundaryPatch<GridView> right;
+  BoundaryPatch<GridView> upper;
+  MyFaces(GridView const &gridView);
+
+private:
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14;
+  };
+};
+#endif

src/sand-wedge-data/mygrid_tmpl.cc

+#ifndef DIM
+#error DIM unset
+#endif
+
+#include "explicitgrid.hh"
+
+template std::shared_ptr<Grid> constructGrid();
+
+template struct MyFaces<GridView>;

src/sand-wedge-data/parset.cfg

+# -*- mode:conf -*-
+gravity         = 9.81  # [m/s^2]
+
+[io]
+printProgress   = false
+writeVTK        = false
+
+[problem]
+finalTime       = 1800  # [s]
+
+[body]
+bulkModulus     = 1e5   # [Pa]
+poissonRatio    = 0.3   # [1]      0.2 - 0.3
+[body.elastic]
+density         = 900   # [kg/m^3]
+shearViscosity  = 1e3   # [Pas]    0
+bulkViscosity   = 1e3   # [Pas]    0
+[body.viscoelastic]
+density         = 1000  # [kg/m^3]
+shearViscosity  = 1e4   # [Pas]
+bulkViscosity   = 1e4   # [Pas]
+
+[boundary.friction]
+C               = 10    # [Pa]
+mu0             = 0.7   # [1]
+V0              = 5e-5  # [m/s]
+L               = 2e-5  # [m]      ?
+initialLogState = 0     # [ ]      ?
+stateModel      = Dieterich
+[boundary.friction.weakening]
+a               = 0.015 # [1]      ?
+b               = 0.030 # [1]      ?
+[boundary.friction.strengthening]
+a               = 0.030 # [1]      ?
+b               = 0.015 # [1]      ?
+
+[timeSteps]
+number = 100000
+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
+
+[solver.tnnmg.linear]
+maxiumumIterations = 100000
+tolerance          = 1e-10
+pre                = 3
+cycle              = 1  # 1 = V, 2 = W, etc.
+post               = 3
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0
+
+[localsolver]
+steps = 100

src/sand-wedge-data/patchfunction.hh

+#ifndef SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
+#define SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+class PatchFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, 2>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  bool static isBetween(double x, double lower, double upper) {
+    return lower <= x and x <= upper;
+  }
+
+  bool static isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14;
+  };
+
+  bool insideRegion2(Dune::FieldVector<double, 2> const &z) const {
+    assert(isClose(0, z[1]));
+    return isBetween(z[0], _X[0], _Y[0]);
+  }
+
+  Dune::FieldVector<double, 2> const &_X;
+  Dune::FieldVector<double, 2> const &_Y;
+
+  double const _v1;
+  double const _v2;
+
+public:
+  PatchFunction(double v1, double v2)
+      : _X(MyGeometry::X), _Y(MyGeometry::Y), _v1(v1), _v2(v2) {}
+
+  void evaluate(Dune::FieldVector<double, 2> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = insideRegion2(x) ? _v2 : _v1;
+  }
+};
+#endif

src/sand-wedge-data/twopiece.hh

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

src/sand-wedge.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 datadir
+#error datadir 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 "sand-wedge-data/mybody.hh"
+#include "sand-wedge-data/mygeometry.hh"
+#include "sand-wedge-data/mygeometry.cc" // FIXME
+#include "sand-wedge-data/myglobalfrictiondata.hh"
+#include "sand-wedge-data/mygrid.cc" // FIXME
+#include "sand-wedge-data/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('" datadir "')");
+}
+
+template <class GridView> class SpecialWriter {
+  using LocalVector = Dune::FieldVector<double, dims>;
+  using Element = typename GridView::Grid::template Codim<0>::Entity;
+  using ElementPointer =
+      typename GridView::Grid::template Codim<0>::EntityPointer;
+
+  void writeMagnitude(LocalVector const &v) {
+    writer_ << v[0] << " "; // Because of Dirichlet conditions, this is like a
+                            // signed norm
+  }
+
+  void writeHorizontal(LocalVector const &v) {
+    writer_ << MyGeometry::horizontalProjection(v) << " ";
+  }
+  void writeVertical(LocalVector const &v) {
+    writer_ << MyGeometry::verticalProjection(v) << " ";
+  }
+
+  std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
+      const {
+    auto const element = hsearch_.findEntity(x);
+    return std::pair<ElementPointer, LocalVector>(element,
+                                                  element->geometry().local(x));
+  }
+
+  std::fstream writer_;
+
+  Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
+
+  std::pair<ElementPointer, LocalVector> const G;
+  std::pair<ElementPointer, LocalVector> const H;
+  std::pair<ElementPointer, LocalVector> const J;
+  std::pair<ElementPointer, LocalVector> const I;
+  std::pair<ElementPointer, LocalVector> const U;
+  std::pair<ElementPointer, LocalVector> const Z;
+
+public:
+  SpecialWriter(std::string filename, GridView const &gridView)
+      : writer_(filename, std::fstream::out),
+        hsearch_(gridView.grid(), gridView),
+        G(globalToLocal(MyGeometry::G)),
+        H(globalToLocal(MyGeometry::H)),
+        J(globalToLocal(MyGeometry::J)),
+        I(globalToLocal(MyGeometry::I)),
+        U(globalToLocal(MyGeometry::U)),
+        Z(globalToLocal(MyGeometry::Z)) {
+    writer_ << "|G| |H| |J| |I| Uv Uh Zv Zh" << std::endl;
+  }
+
+  void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
+                 specialField) {
+    LocalVector value;
+
+    specialField.evaluateLocal(*G.first, G.second, value);
+    writeMagnitude(value);
+
+    specialField.evaluateLocal(*H.first, H.second, value);
+    writeMagnitude(value);
+
+    specialField.evaluateLocal(*J.first, J.second, value);
+    writeMagnitude(value);
+
+    specialField.evaluateLocal(*I.first, I.second, value);
+    writeMagnitude(value);
+
+    specialField.evaluateLocal(*U.first, U.second, value);
+    writeVertical(value);
+    writeHorizontal(value);
+
+    specialField.evaluateLocal(*Z.first, Z.second, value);
+    writeVertical(value);
+    writeHorizontal(value);
+
+    writer_ << std::endl;
+  }
+};
+
+int main(int argc, char *argv[]) {
+  try {
+    Dune::ParameterTree parset;
+    Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset);
+    Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+
+    MyGeometry::render();
+    MyGeometry::write();
+
+    // {{{ Set up grid
+    using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
+    auto grid = constructGrid<Grid>();
+
+    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();
+
+    MyFaces<GridView> myFaces(leafView);
+
+    // 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.right.containsVertex(i))
+        dirichletNodes[i][0] = true;
+
+      if (myFaces.lower.containsVertex(i))
+        dirichletNodes[i][1] = 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("boundaryconditions")
+          .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;
+    using LocalScalarVector = ScalarVector::block_type;
+
+    MyAssembler myAssembler(leafView);
+
+    MyBody<dims> const body(parset);
+
+    Matrix A, C, M;
+    myAssembler.assembleElasticity(body.getYoungModulus(),
+                                   body.getPoissonRatio(), A);
+    myAssembler.assembleViscosity(body.getShearViscosityField(),
+                                  body.getBulkViscosityField(), 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);
+
+    // {{{ Initial conditions
+    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"));
+
+    ScalarVector alpha_initial(fineVertexCount);
+    alpha_initial = parset.get<double>("boundary.friction.initialLogState");
+    ScalarVector normalStress(fineVertexCount);
+    myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
+                                     body.getYoungModulus(),
+                                     body.getPoissonRatio(), u_initial);
+    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
+    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
+        frictionalBoundary, frictionInfo, normalStress);
+    myGlobalNonlinearity->updateLogState(alpha_initial);
+
+    Vector v_initial(fineVertexCount);
+    v_initial = 0.0;
+    {
+      double v_initial_const;
+      velocityDirichletFunction.evaluate(0.0, v_initial_const);
+      assert(v_initial_const == 0.0);
+    }
+
+    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 0 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> frictionWriter(
+        vertexCoordinates, frictionalNodes,
+        [](LocalVector const &x) { return x[0]; });
+    auto const reportFriction = [&](Vector const &_u, Vector const &_v,
+                                    ScalarVector const &_alpha) {
+      ScalarVector c;
+      myGlobalNonlinearity->coefficientOfFriction(_v, c);
+      frictionWriter.writeKinetics(_u, _v);
+      frictionWriter.writeOther(c, _alpha);
+    };
+    reportFriction(u_initial, v_initial, alpha_initial);
+
+    MyVTKWriter<typename MyAssembler::VertexBasis,
+                typename MyAssembler::CellBasis> const
+    vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
+
+    if (parset.get<bool>("io.writeVTK")) {
+      ScalarVector stress;
+      myAssembler.assembleVonMisesStress(
+          body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
+      vtkWriter.write(0, u_initial, v_initial, alpha_initial, stress);
+    }
+
+    SpecialWriter<GridView> specialVelocityWriter("specialVelocities",
+                                                  leafView);
+    SpecialWriter<GridView> specialDisplacementWriter("specialDisplacements",
+                                                      leafView);
+
+    // 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;
+    Vector v_m(fineVertexCount);
+    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 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;
+      };
+
+      Vector u;
+      Vector v_saved;
+      ScalarVector alpha_saved;
+      double lastStateCorrection;
+      for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
+        timeSteppingScheme->extractOldVelocity(v_m);
+        v_m *= 0.5;
+        Arithmetic::addProduct(v_m, 0.5, v);
+
+        stateUpdater->solve(v_m);
+        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(v_saved, v);
+          if (velocityCorrection < fixedPointTolerance)
+            break;
+        }
+        if (stateFPI == maximumStateFPI)
+          DUNE_THROW(Dune::Exception, "FPI failed to converge");
+
+        alpha_saved = alpha;
+        v_saved = v;
+      }
+      if (printProgress)
+        std::cout << std::endl;
+
+      reportFriction(u, v, alpha);
+      iterationWriter << std::endl;
+      relaxationWriter << std::endl;
+
+      {
+        BasisGridFunction<typename MyAssembler::VertexBasis, Vector> velocity(
+            myAssembler.vertexBasis, v);
+        BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
+        displacement(myAssembler.vertexBasis, u);
+        specialVelocityWriter.write(velocity);
+        specialDisplacementWriter.write(displacement);
+      }
+
+      if (parset.get<bool>("io.writeVTK")) {
+        ScalarVector stress;
+        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
+                                           body.getPoissonRatio(), u, stress);
+        vtkWriter.write(timeStep, 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;
+  }
+}