src/hdf5/restart-io.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "restart-io.hh"
+
+template <class ProgramState>
+RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
+    : displacementWriter_(group, "displacement", vertexCount,
+                          Vector::block_type::dimension),
+      velocityWriter_(group, "velocity", vertexCount,
+                      Vector::block_type::dimension),
+      accelerationWriter_(group, "acceleration", vertexCount,
+                          Vector::block_type::dimension),
+      stateWriter_(group, "state", vertexCount),
+      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
+      relativeTimeWriter_(group, "relativeTime"),
+      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
+
+template <class ProgramState>
+void RestartIO<ProgramState>::write(ProgramState const &programState) {
+  addEntry(displacementWriter_, programState.timeStep, programState.u);
+  addEntry(velocityWriter_, programState.timeStep, programState.v);
+  addEntry(accelerationWriter_, programState.timeStep, programState.a);
+  addEntry(stateWriter_, programState.timeStep, programState.alpha);
+  addEntry(weightedNormalStressWriter_, programState.timeStep,
+           programState.weightedNormalStress);
+  addEntry(relativeTimeWriter_, programState.timeStep,
+           programState.relativeTime);
+  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+           programState.relativeTau);
+}
+
+template <class ProgramState>
+void RestartIO<ProgramState>::read(size_t timeStep,
+                                   ProgramState &programState) {
+  programState.timeStep = timeStep;
+  readEntry(displacementWriter_, timeStep, programState.u);
+  readEntry(velocityWriter_, timeStep, programState.v);
+  readEntry(accelerationWriter_, timeStep, programState.a);
+  readEntry(stateWriter_, timeStep, programState.alpha);
+  readEntry(weightedNormalStressWriter_, timeStep,
+            programState.weightedNormalStress);
+  readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
+  readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
+}
+
+#include "restart-io_tmpl.cc"

src/hdf5/restart-io.hh

+#ifndef SRC_HDF_RESTART_HDF_HH
+#define SRC_HDF_RESTART_HDF_HH
+
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+template <class ProgramState> class RestartIO {
+  using ScalarVector = typename ProgramState::ScalarVector;
+  using Vector = typename ProgramState::Vector;
+
+public:
+  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
+
+  void write(ProgramState const &programState);
+
+  void read(size_t timeStep, ProgramState &programState);
+
+private:
+  HDF5::SequenceIO<2> displacementWriter_;
+  HDF5::SequenceIO<2> velocityWriter_;
+  HDF5::SequenceIO<2> accelerationWriter_;
+  HDF5::SequenceIO<1> stateWriter_;
+  HDF5::SequenceIO<1> weightedNormalStressWriter_;
+  HDF5::SequenceIO<0> relativeTimeWriter_;
+  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
+};
+#endif

src/hdf5/restart-io_tmpl.cc

+#include "../explicitvectors.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class RestartIO<MyProgramState>;

src/hdf5/restrict.hh

+#ifndef SRC_HDF5_RESTRICT_HH
+#define SRC_HDF5_RESTRICT_HH
+
+#include <cassert>
+
+#include <dune/common/bitsetvector.hh>
+
+#include "../tobool.hh"
+
+template <class Vector, class Patch>
+Vector restrictToSurface(Vector const &v1, Patch const &patch) {
+  auto const &vertices = *patch.getVertices();
+  assert(vertices.size() == v1.size());
+
+  Vector ret(vertices.count());
+  auto target = ret.begin();
+  for (size_t i = 0; i < v1.size(); ++i)
+    if (toBool(vertices[i]))
+      *(target++) = v1[i];
+  assert(target == ret.end());
+  return ret;
+}
+#endif

src/hdf5/surface-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "restrict.hh"
+#include "surface-writer.hh"
+
+template <class ProgramState, class GridView>
+SurfaceWriter<ProgramState, GridView>::SurfaceWriter(
+    HDF5::Grouplike &file, Vector const &vertexCoordinates, Patch const &surface)
+    : group_(file, "surface"),
+      surface_(surface),
+      surfaceDisplacementWriter_(group_, "displacement", surface.numVertices(),
+                                 Vector::block_type::dimension),
+      surfaceVelocityWriter_(group_, "velocity", surface.numVertices(),
+                             Vector::block_type::dimension) {
+  auto const surfaceCoordinates = restrictToSurface(vertexCoordinates, surface);
+  HDF5::SingletonWriter<2> surfaceCoordinateWriter(group_, "coordinates",
+                                                 surfaceCoordinates.size(),
+                                                 Vector::block_type::dimension);
+  setEntry(surfaceCoordinateWriter, surfaceCoordinates);
+}
+
+template <class ProgramState, class GridView>
+void SurfaceWriter<ProgramState, GridView>::write(
+    ProgramState const &programState) {
+  auto const surfaceDisplacements = restrictToSurface(programState.u, surface_);
+  addEntry(surfaceDisplacementWriter_, programState.timeStep,
+           surfaceDisplacements);
+
+  auto const surfaceVelocities = restrictToSurface(programState.v, surface_);
+  addEntry(surfaceVelocityWriter_, programState.timeStep, surfaceVelocities);
+}
+
+#include "surface-writer_tmpl.cc"

src/hdf5/surface-writer.hh

+#ifndef SRC_HDF_SURFACE_WRITER_HH
+#define SRC_HDF_SURFACE_WRITER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/fufem/hdf5/sequenceio.hh>
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+template <class ProgramState, class GridView> class SurfaceWriter {
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
+                Patch const &surface);
+
+  void write(ProgramState const &programState);
+
+private:
+  HDF5::Group group_;
+
+  Patch const &surface_;
+
+  HDF5::SequenceIO<2> surfaceDisplacementWriter_;
+  HDF5::SequenceIO<2> surfaceVelocityWriter_;
+};
+#endif

src/hdf5/surface-writer_tmpl.cc

+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class SurfaceWriter<MyProgramState, GridView>;

src/hdf5/time-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "time-writer.hh"
+
+template <class ProgramState>
+TimeWriter<ProgramState>::TimeWriter(HDF5::Grouplike &file)
+    : file_(file),
+      relativeTimeWriter_(file_, "relativeTime"),
+      relativeTimeIncrementWriter_(file_, "relativeTimeIncrement") {}
+
+template <class ProgramState>
+void TimeWriter<ProgramState>::write(ProgramState const &programState) {
+  addEntry(relativeTimeWriter_, programState.timeStep,
+           programState.relativeTime);
+  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
+           programState.relativeTau);
+}
+
+#include "time-writer_tmpl.cc"

src/hdf5/time-writer.hh

+#ifndef SRC_HDF_TIME_WRITER_HH
+#define SRC_HDF_TIME_WRITER_HH
+
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+template <class ProgramState> class TimeWriter {
+public:
+  TimeWriter(HDF5::Grouplike &file);
+  void write(ProgramState const &programState);
+
+private:
+  HDF5::Grouplike &file_;
+  HDF5::SequenceIO<0> relativeTimeWriter_;
+  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
+};
+#endif

src/hdf5/time-writer_tmpl.cc

+#include "../explicitvectors.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class TimeWriter<MyProgramState>;

src/matrices.hh

+#ifndef SRC_MATRICES_HH
+#define SRC_MATRICES_HH
+
+template <class Matrix> struct Matrices {
+  Matrix elasticity;
+  Matrix damping;
+  Matrix mass;
+};
+#endif

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

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

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/one-body-problem-data/mybody.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
+#define SRC_ONE_BODY_PROBLEM_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 "mygeometry.hh"
+#include "segmented-function.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_;
+  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