src/hdf5/iteration-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "iteration-writer.hh"
+
+IterationWriter::IterationWriter(HDF5::Grouplike &file)
+    : group_(file, "iterations"),
+      fpiSubGroup_(group_, "fixedPoint"),
+      mgSubGroup_(group_, "multiGrid"),
+      finalMGIterationWriter_(mgSubGroup_, "final"),
+      finalFPIIterationWriter_(fpiSubGroup_, "final"),
+      totalMGIterationWriter_(mgSubGroup_, "total"),
+      totalFPIIterationWriter_(fpiSubGroup_, "total") {}
+
+void IterationWriter::write(size_t timeStep,
+                            IterationRegister const &iterationCount) {
+  addEntry(finalMGIterationWriter_, timeStep,
+           iterationCount.finalCount.multigridIterations);
+  addEntry(finalFPIIterationWriter_, timeStep,
+           iterationCount.finalCount.iterations);
+  addEntry(totalMGIterationWriter_, timeStep,
+           iterationCount.totalCount.multigridIterations);
+  addEntry(totalFPIIterationWriter_, timeStep,
+           iterationCount.totalCount.iterations);
+}

src/hdf5/iteration-writer.hh

+#ifndef SRC_HDF_ITERATION_WRITER_HH
+#define SRC_HDF_ITERATION_WRITER_HH
+
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+#include "../time-stepping/adaptivetimestepper.hh"
+
+class IterationWriter {
+public:
+  IterationWriter(HDF5::Grouplike &file);
+
+  void write(size_t timeStep, IterationRegister const &iterationCount);
+
+private:
+  HDF5::Group group_;
+  HDF5::Group fpiSubGroup_;
+  HDF5::Group mgSubGroup_;
+
+  HDF5::SequenceIO<0, size_t> finalMGIterationWriter_;
+  HDF5::SequenceIO<0, size_t> finalFPIIterationWriter_;
+  HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
+  HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
+};
+#endif

src/hdf5/patchinfo-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/grid/hierarchic-approximation.hh>
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+#include "patchinfo-writer.hh"
+
+template <class LocalVector, class GridView>
+GridEvaluator<LocalVector, GridView>::GridEvaluator(
+    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
+  double const bufferWidth = 0.05 * MyGeometry::lengthScale;
+  auto const xminmax = std::minmax_element(
+      weakPatch.vertices.begin(), weakPatch.vertices.end(),
+      [](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
+  double const xmin = (*xminmax.first)[0] - bufferWidth;
+  double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;
+
+  double const zmin = -MyGeometry::depth / 2.0;
+  double const zspan = MyGeometry::depth;
+
+  double spatialResolution = 0.01 * MyGeometry::lengthScale;
+  size_t const xsteps = std::round(xspan / spatialResolution);
+  size_t const zsteps = std::round(zspan / spatialResolution);
+
+  xCoordinates.resize(xsteps + 1);
+  for (size_t xi = 0; xi <= xsteps; xi++)
+    xCoordinates[xi] = xmin + xi * xspan / xsteps;
+  zCoordinates.resize(zsteps + 1);
+  for (size_t zi = 0; zi <= zsteps; zi++)
+    zCoordinates[zi] = zmin + zi * zspan / zsteps;
+
+  HierarchicApproximation<typename GridView::Grid, GridView> const
+      hApproximation(gridView.grid(), gridView, 1e-6 * MyGeometry::lengthScale);
+
+  LocalVector global(0);
+  localInfo.resize(xsteps + 1);
+  for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
+    localInfo[xi].resize(zsteps + 1);
+    for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
+      global[0] = xCoordinates[xi];
+      global[2] = zCoordinates[zi];
+      auto const element = hApproximation.findEntity(global);
+      localInfo[xi][zi] =
+          std::make_pair(element, element.geometry().local(global));
+    }
+  }
+}
+
+template <class LocalVector, class GridView>
+template <class Function>
+Dune::Matrix<typename Function::RangeType>
+GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
+  Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
+                                                 zCoordinates.size());
+  for (size_t xi = 0; xi < localInfo.size(); ++xi) {
+    auto const &localInfoX = localInfo[xi];
+    for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
+      auto const &localInfoXZ = localInfoX[zi];
+      f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
+    }
+  }
+  return ret;
+}
+
+template <class ProgramState, class VertexBasis, class GridView>
+PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
+    HDF5::Grouplike &file, VertexBasis const &vertexBasis,
+    Patch const &frictionalBoundary,
+    ConvexPolyhedron<LocalVector> const &weakPatch)
+    : group_(file, "weakPatchGrid"),
+      vertexBasis_(vertexBasis),
+      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
+      weakPatchGridVelocityWriter_(
+          group_, "velocity", gridEvaluator_.xCoordinates.size(),
+          gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
+  HDF5::SingletonWriter<1> weakPatchGridXCoordinateWriter(
+      group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
+  setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);
+
+  HDF5::SingletonWriter<1> weakPatchGridZCoordinateWriter(
+      group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
+  setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
+}
+
+template <class ProgramState, class VertexBasis, class GridView>
+void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
+    ProgramState const &programState) {
+  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
+  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
+  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
+}
+
+#include "patchinfo-writer_tmpl.cc"

src/hdf5/patchinfo-writer.hh

+#ifndef SRC_HDF_PATCHINFO_WRITER_HH
+#define SRC_HDF_PATCHINFO_WRITER_HH
+
+#include <dune/istl/matrix.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+#include "../one-body-problem-data/mygeometry.hh"
+
+template <class LocalVector, class GridView> class GridEvaluator {
+  using Element = typename GridView::Grid::template Codim<0>::Entity;
+
+public:
+  GridEvaluator(ConvexPolyhedron<LocalVector> const &weakPatch,
+                GridView const &gridView);
+
+  template <class Function>
+  Dune::Matrix<typename Function::RangeType> evaluate(Function const &f) const;
+
+  Dune::BlockVector<Dune::FieldVector<double, 1>> xCoordinates;
+  Dune::BlockVector<Dune::FieldVector<double, 1>> zCoordinates;
+
+private:
+  std::vector<std::vector<std::pair<Element, LocalVector>>> localInfo;
+};
+
+template <class ProgramState, class VertexBasis, class GridView>
+class PatchInfoWriter {
+  using Vector = typename ProgramState::Vector;
+  using LocalVector = typename Vector::block_type;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
+                  Patch const &frictionalBoundary,
+                  ConvexPolyhedron<LocalVector> const &weakPatch);
+
+  void write(ProgramState const &programState);
+
+private:
+  HDF5::Group group_;
+
+  VertexBasis const &vertexBasis_;
+
+  GridEvaluator<LocalVector, GridView> const gridEvaluator_;
+  HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
+};
+#endif

src/hdf5/patchinfo-writer_tmpl.cc

+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+using P1Basis = P1NodalBasis<GridView, double>;
+using MyFunction = BasisGridFunction<P1Basis, Vector>;
+
+template class GridEvaluator<LocalVector, GridView>;
+template Dune::Matrix<typename MyFunction::RangeType>
+GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
+
+template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;

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}