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/sand-wedge-data/midpoint.hh → src/one-body-problem-data/midpoint.hh

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

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

-#ifndef SRC_SAND_WEDGE_DATA_MYBODY_HH
-#define SRC_SAND_WEDGE_DATA_MYBODY_HH
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYBODY_HH
 
 #include <dune/common/fvector.hh>
 
@@ -8,8 +8,8 @@
 #include <dune/tectonic/body.hh>
 #include <dune/tectonic/gravity.hh>
 
-#include "twopiece.hh"
 #include "mygeometry.hh"
+#include "segmented-function.hh"
 
 template <int dimension> class MyBody : public Body<dimension> {
   using typename Body<dimension>::ScalarFunction;
@@ -47,9 +47,9 @@ template <int dimension> class MyBody : public Body<dimension> {
 private:
   double const poissonRatio_;
   double const youngModulus_;
-  TwoPieceFunction const shearViscosityField_;
-  TwoPieceFunction const bulkViscosityField_;
-  TwoPieceFunction const densityField_;
+  SegmentedFunction const shearViscosityField_;
+  SegmentedFunction const bulkViscosityField_;
+  SegmentedFunction const densityField_;
   Gravity<dimension> const gravityField_;
 };
 #endif

src/sand-wedge-data/mygeometry.cc → src/one-body-problem-data/mygeometry.cc

@@ -12,13 +12,6 @@
 
 #include "mygeometry.hh"
 
-double MyGeometry::verticalProjection(LocalVector const &x) {
-  return zenith * x;
-}
-double MyGeometry::horizontalProjection(LocalVector const &x) {
-  return orthogonalZenith * x;
-}
-
 void MyGeometry::write() {
   std::fstream writer("geometry", std::fstream::out);
   writer << "A = " << A << std::endl;
@@ -132,8 +125,7 @@ void MyGeometry::render() {
 
   // labels
   {
-    auto const label = [&](LocalVector2D const &v,
-                           std::string l) { // TODO: implicit?
+    auto const label = [&](LocalVector2D const &v, std::string l) {
       moveTo(v);
       cr->rel_move_to(0.005, -0.02);
       cr->show_text(l);

src/sand-wedge-data/mygeometry.hh → src/one-body-problem-data/mygeometry.hh

@@ -11,74 +11,60 @@ namespace {
   using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;
 
   using LocalVector = Dune::FieldVector<double, MY_DIM>;
-
-  // kludge because fieldvectors have no initialiser_list constructor, see
-  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
-  LocalVector2D createVector(double x, double y) {
-    LocalVector2D ret{ 0 };
-    ret[0] = x;
-    ret[1] = y;
-    return ret;
-  }
-
-  LocalMatrix2D createMatrix(double a11, double a12, double a21, double a22) {
-    LocalMatrix2D ret{ 0 };
-    ret[0][0] = a11;
-    ret[0][1] = a12;
-    ret[1][0] = a21;
-    ret[1][1] = a22;
-    return ret;
-  }
 }
 
 namespace reference {
-  double const rightLeg = 0.27;
-  double const leftLeg = 1.00;
+  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; // Height of the viscous bottom layer
-  double const weakLen = 0.20;     // Length of the weak zone
+  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 = createVector(0, 0);
-  LocalVector2D const B = createVector(leftLeg, -rightLeg);
-  LocalVector2D const C = createVector(leftLeg, 0);
+  LocalVector2D const A = {0, 0};
+  LocalVector2D const B = {leftLeg, -rightLeg};
+  LocalVector2D const C = {leftLeg, 0};
 
-  LocalVector2D const Z = createVector(0.35, 0);
-  LocalVector2D const Y = createVector(0.35, -0.35 / leftLeg * rightLeg);
-  LocalVector2D const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
-                                       Y[1] + weakLen * std::sin(leftAngle));
+  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 = createVector(X[0], 0);
+  LocalVector2D const U = {X[0], 0};
 
-  LocalVector2D const K =
-      createVector(B[0] - leftLeg * viscoHeight / rightLeg, B[1] + viscoHeight);
-  LocalVector2D const M = createVector(B[0], B[1] + viscoHeight);
+  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 = createVector(Y[0] + G[0], Y[1] + G[1]);
+  LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
 
-  LocalVector2D const zenith = createVector(0, 1);
-  LocalVector2D const orthogonalZenith = createVector(-1, 0);
+  LocalVector2D const zenith = {0, 1};
 
-  LocalMatrix2D const rotation =
-      createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
-                   std::sin(leftAngle), std::cos(leftAngle));
+  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 };
+    LocalVector ret(0);
     ret[0] = ret2D[0];
     ret[1] = ret2D[1];
     return ret;
   }
 }
 
-double const depth = 0.10;
+double const lengthScale = reference::s;
+
+double const depth = 0.60 * lengthScale;
 
 LocalVector const A = rotate(reference::A);
 LocalVector const B = rotate(reference::B);
@@ -95,10 +81,6 @@ LocalVector const Y = rotate(reference::Y);
 LocalVector const Z = rotate(reference::Z);
 
 LocalVector const zenith = rotate(reference::zenith);
-LocalVector const orthogonalZenith = rotate(reference::orthogonalZenith);
-
-double verticalProjection(LocalVector const &);
-double horizontalProjection(LocalVector const &);
 
 void write();
 

src/sand-wedge-data/myglobalfrictiondata.hh → src/one-body-problem-data/myglobalfrictiondata.hh

-#ifndef SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_SAND_WEDGE_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/common/fvector.hh>
 
 #include <dune/tectonic/globalfrictiondata.hh>
 
-#include "mygeometry.hh"
 #include "patchfunction.hh"
 
-template <int dimension>
-class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+template <class LocalVector>
+class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
 private:
-  using typename GlobalFrictionData<dimension>::VirtualFunction;
+  using typename GlobalFrictionData<LocalVector::dimension>::VirtualFunction;
 
 public:
-  MyGlobalFrictionData(Dune::ParameterTree const &parset)
+  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")),
+           parset.get<double>("weakening.a"), segment),
         b_(parset.get<double>("strengthening.b"),
-           parset.get<double>("weakening.b")),
+           parset.get<double>("weakening.b"), segment),
         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_; }
+  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_;

src/sand-wedge-data/mygrid.cc → src/one-body-problem-data/mygrid.cc

@@ -2,8 +2,11 @@
 #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) {}
@@ -41,7 +44,9 @@ void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
 #endif
 }
 
-auto SimplexManager::getSimplices() -> SimplexList const &{ return simplices_; }
+auto SimplexManager::getSimplices() -> SimplexList const & {
+  return simplices_;
+}
 
 template <class Grid> GridConstructor<Grid>::GridConstructor() {
   auto const &A = MyGeometry::A;
@@ -120,7 +125,7 @@ template <class GridView>
 template <class Vector>
 bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
                                     Vector const &c) {
-  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
 }
 
 template <class GridView>
@@ -130,9 +135,11 @@ bool MyFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
   auto const minmax0 = std::minmax(v1[0], v2[0]);
   auto const minmax1 = std::minmax(v1[1], v2[1]);
 
-  if (minmax0.first - 1e-14 > x[0] or x[0] > minmax0.second + 1e-14)
+  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 > x[1] or x[1] > minmax1.second + 1e-14)
+  if (minmax1.first - 1e-14 * MyGeometry::lengthScale > x[1] or
+      x[1] > minmax1.second + 1e-14 * MyGeometry::lengthScale)
     return false;
 
   return true;
@@ -160,7 +167,6 @@ MyFaces<GridView>::MyFaces(GridView const &gridView)
       upper(gridView)
 #endif
 {
-  using Grid = typename GridView::Grid;
   assert(isClose(MyGeometry::A[1], 0));
   assert(isClose(MyGeometry::B[1], 0));
   lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
@@ -185,4 +191,37 @@ MyFaces<GridView>::MyFaces(GridView const &gridView)
   });
 }
 
+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/sand-wedge-data/mygrid.hh → src/one-body-problem-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>
+#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>
 
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wdeprecated-declarations"
 #include <dune/fufem/boundarypatch.hh>
-#pragma clang diagnostic pop
+#include <dune/fufem/geometry/convexpolyhedron.hh>
 
 #include "mygeometry.hh"
 
@@ -28,7 +23,12 @@ template <class GridView> struct MyFaces {
 
 private:
   bool isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14;
+    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>
@@ -74,4 +74,11 @@ template <class Grid> class GridConstructor {
 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/sand-wedge-data/mygrid_tmpl.cc → src/one-body-problem-data/mygrid_tmpl.cc

@@ -2,7 +2,8 @@
 #error MY_DIM unset
 #endif
 
-#include "explicitgrid.hh"
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
 
 template class GridConstructor<Grid>;
 
@@ -10,3 +11,7 @@ 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);