src/diameter.hh

+#ifndef SRC_DIAMETER_HH
+#define SRC_DIAMETER_HH
+
+template <class Geometry> double diameter(Geometry const &geometry) {
+  auto const numCorners = geometry.corners();
+  std::vector<typename Geometry::GlobalCoordinate> corners(numCorners);
+
+  double diameter = 0.0;
+  for (int i = 0; i < numCorners; ++i) {
+    corners[i] = geometry.corner(i);
+    for (int j = 0; j < i; ++j)
+      diameter = std::max(diameter, (corners[i] - corners[j]).two_norm());
+  }
+  return diameter;
+}
+#endif

src/enum_parser.cc

-// Copyright Carsten Graeser 2012
-
-#include <dune/common/exceptions.hh>
-#include <dune/common/typetraits.hh>
-
-template <class EnumType> struct StringToEnum : public Dune::NotImplemented {};
-
-template <class EnumType>
-typename Dune::enable_if<
-    !Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<EnumType>>::value,
-    std::istream &>::type
-operator>>(std::istream &lhs, EnumType &e) {
-  std::string s;
-  lhs >> s;
-
-  try {
-    e = StringToEnum<EnumType>::convert(s);
-  }
-  catch (typename Dune::Exception) {
-    lhs.setstate(std::ios_base::failbit);
-  }
-  return lhs;
-}

src/enum_scheme.cc

-#include <dune/common/exceptions.hh>
-
-template <> struct StringToEnum<Config::scheme> {
-  static Config::scheme convert(std::string const &s) {
-    if (s == "newmark")
-      return Config::Newmark;
-
-    if (s == "backwardEuler")
-      return Config::BackwardEuler;
-
-    DUNE_THROW(Dune::Exception, "failed to parse enum");
-  }
-};

src/enum_state_model.cc

-#include <dune/common/exceptions.hh>
-
-template <> struct StringToEnum<Config::stateModel> {
-  static Config::stateModel convert(std::string const &s) {
-    if (s == "Dieterich")
-      return Config::Dieterich;
-
-    if (s == "Ruina")
-      return Config::Ruina;
-
-    DUNE_THROW(Dune::Exception, "failed to parse enum");
-  }
-};

src/enum_verbosity.cc

-#include <dune/solvers/common/numproc.hh> // Solver::VerbosityMode
-#include <dune/common/exceptions.hh>
-
-template <> struct StringToEnum<Solver::VerbosityMode> {
-  static Solver::VerbosityMode convert(std::string const &s) {
-    if (s == "full")
-      return Solver::FULL;
-
-    if (s == "reduced")
-      return Solver::REDUCED;
-
-    if (s == "quiet")
-      return Solver::QUIET;
-
-    DUNE_THROW(Dune::Exception, "failed to parse enum");
-  }
-};

src/enumparser.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+// Copyright Carsten Graeser 2012
+
+#include <dune/common/exceptions.hh>
+
+#include "enumparser.hh"
+
+template <class Enum>
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
+    std::istream &>::type
+operator>>(std::istream &lhs, Enum &e) {
+  std::string s;
+  lhs >> s;
+
+  try {
+    e = StringToEnum<Enum>::convert(s);
+  } catch (typename Dune::Exception) {
+    lhs.setstate(std::ios_base::failbit);
+  }
+  return lhs;
+}
+
+Config::FrictionModel StringToEnum<Config::FrictionModel>::convert(
+    std::string const &s) {
+  if (s == "Truncated")
+    return Config::Truncated;
+
+  if (s == "Regularised")
+    return Config::Regularised;
+
+  DUNE_THROW(Dune::Exception, "failed to parse enum");
+}
+
+Config::scheme StringToEnum<Config::scheme>::convert(std::string const &s) {
+  if (s == "newmark")
+    return Config::Newmark;
+
+  if (s == "backwardEuler")
+    return Config::BackwardEuler;
+
+  DUNE_THROW(Dune::Exception, "failed to parse enum");
+}
+
+Config::stateModel StringToEnum<Config::stateModel>::convert(
+    std::string const &s) {
+  if (s == "AgeingLaw")
+    return Config::AgeingLaw;
+
+  if (s == "SlipLaw")
+    return Config::SlipLaw;
+
+  DUNE_THROW(Dune::Exception, "failed to parse enum");
+}
+
+Config::PatchType StringToEnum<Config::PatchType>::convert(
+    std::string const &s) {
+  if (s == "Rectangular")
+    return Config::Rectangular;
+
+  if (s == "Trapezoidal")
+    return Config::Trapezoidal;
+
+  DUNE_THROW(Dune::Exception, "failed to parse enum");
+}
+
+template std::istream &operator>>(std::istream &lhs, Config::FrictionModel &);
+template std::istream &operator>>(std::istream &lhs, Config::stateModel &);
+template std::istream &operator>>(std::istream &lhs, Config::scheme &);
+template std::istream &operator>>(std::istream &lhs, Config::PatchType &);

src/enumparser.hh

+#ifndef SRC_ENUMPARSER_HH
+#define SRC_ENUMPARSER_HH
+
+// Copyright Carsten Graeser 2012
+
+#include <type_traits>
+
+#include <dune/solvers/solvers/solver.hh>
+
+#include "enums.hh"
+
+template <class Enum> struct StringToEnum : public Dune::NotImplemented {};
+
+template <> struct StringToEnum<Config::FrictionModel> {
+  static Config::FrictionModel convert(std::string const &s);
+};
+
+template <> struct StringToEnum<Config::stateModel> {
+  static Config::stateModel convert(std::string const &s);
+};
+
+template <> struct StringToEnum<Config::scheme> {
+  static Config::scheme convert(std::string const &s);
+};
+
+template <> struct StringToEnum<Config::PatchType> {
+  static Config::PatchType convert(std::string const &s);
+};
+
+template <class Enum>
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
+    std::istream &>::type
+operator>>(std::istream &lhs, Enum &e);
+#endif

src/enums.hh

-#ifndef ENUMS_HH
-#define ENUMS_HH
+#ifndef SRC_ENUMS_HH
+#define SRC_ENUMS_HH
 
 struct Config {
-  enum stateModel {
-    Dieterich,
-    Ruina
-  };
-  enum scheme {
-    Newmark,
-    BackwardEuler
-  };
+  enum FrictionModel { Truncated, Regularised };
+  enum stateModel { AgeingLaw, SlipLaw };
+  enum scheme { Newmark, BackwardEuler };
+  enum PatchType { Rectangular, Trapezoidal };
 };
+
 #endif

src/explicitgrid.hh

+#ifndef SRC_EXPLICITGRID_HH
+#define SRC_EXPLICITGRID_HH
+
+#include "gridselector.hh"
+
+using GridView = Grid::LeafGridView;
+#endif

src/explicitvectors.hh

+#ifndef SRC_EXPLICITVECTORS_HH
+#define SRC_EXPLICITVECTORS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+using LocalVector = Dune::FieldVector<double, MY_DIM>;
+using LocalMatrix = Dune::FieldMatrix<double, MY_DIM, MY_DIM>;
+using Vector = Dune::BlockVector<LocalVector>;
+using Matrix = Dune::BCRSMatrix<LocalMatrix>;
+using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+#endif

src/friction_writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "friction_writer.hh"
-
-double computeCOF(FrictionData const &fd, double V, double log_state) {
-  double const mu = fd.mu0 + fd.a * std::log(V / fd.V0) +
-                    fd.b * (log_state + std::log(fd.V0 / fd.L));
-  return (mu > 0) ? mu : 0;
-}
-
-template <class BitVectorType>
-FrictionWriter<BitVectorType>::FrictionWriter(
-    FrictionData const &_fd, BitVectorType const &_boundaryNodes)
-    : coefficientWriter("coefficients", std::fstream::out),
-      displacementWriter("displacements", std::fstream::out),
-      stateWriter("states", std::fstream::out),
-      velocityWriter("velocities", std::fstream::out),
-      fd(_fd),
-      boundaryNodes(_boundaryNodes) {}
-
-template <class BitVectorType>
-FrictionWriter<BitVectorType>::~FrictionWriter() {
-  stateWriter.close();
-  displacementWriter.close();
-  velocityWriter.close();
-  coefficientWriter.close();
-}
-
-template <class BitVectorType>
-template <class SingletonVectorType, class VectorType>
-void FrictionWriter<BitVectorType>::writeInfo(SingletonVectorType const &alpha,
-                                              VectorType const &u,
-                                              VectorType const &v) {
-  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
-    if (!boundaryNodes[i][0])
-      continue;
-
-    coefficientWriter << computeCOF(fd, v[i].two_norm(), alpha[i]) << " ";
-    displacementWriter << u[i][0] << " ";
-    stateWriter << alpha[i][0] << " ";
-    velocityWriter << v[i][0] << " ";
-  }
-
-  stateWriter << std::endl;
-  displacementWriter << std::endl;
-  velocityWriter << std::endl;
-  coefficientWriter << std::endl;
-}
-
-#include "friction_writer_tmpl.cc"

src/friction_writer.hh

-//#include <iostream>
-#include <fstream>
-
-#include <dune/tectonic/frictiondata.hh>
-
-template <class BitVectorType> class FrictionWriter {
-public:
-  FrictionWriter(FrictionData const &_fd, BitVectorType const &_boundaryNodes);
-
-  ~FrictionWriter();
-
-  template <class SingletonVectorType, class VectorType>
-  void writeInfo(SingletonVectorType const &alpha, VectorType const &u,
-                 VectorType const &v);
-
-private:
-  std::fstream coefficientWriter;
-  std::fstream displacementWriter;
-  std::fstream stateWriter;
-  std::fstream velocityWriter;
-  FrictionData const &fd;
-  BitVectorType const &boundaryNodes;
-};

src/friction_writer_tmpl.cc

-#include <dune/common/bitsetvector.hh>
-#include <dune/istl/bvector.hh>
-
-using BitVectorType = Dune::BitSetVector<1>;
-using SingletonVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-using VectorType2 = Dune::BlockVector<Dune::FieldVector<double, 2>>;
-using VectorType3 = Dune::BlockVector<Dune::FieldVector<double, 3>>;
-
-template class FrictionWriter<BitVectorType>;
-
-template void FrictionWriter<BitVectorType>::writeInfo(
-    SingletonVectorType const &alpha, VectorType2 const &u,
-    VectorType2 const &v);
-template void FrictionWriter<BitVectorType>::writeInfo(
-    SingletonVectorType const &alpha, VectorType3 const &u,
-    VectorType3 const &v);

src/gridselector.hh

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#define WANT_ALUGRID 0
+#define WANT_UG 1
+
+#define WANT_GRID WANT_UG
+
+#if WANT_GRID == WANT_ALUGRID
+
+#if !HAVE_ALUGRID
+#error ALUGRID was requested but not found
+#endif
+#include <dune/grid/alugrid.hh>
+using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
+
+#elif WANT_GRID == WANT_UG
+
+#if !HAVE_UG
+#error UG was requested but not found
+#endif
+#include <dune/grid/uggrid.hh>
+using Grid = Dune::UGGrid<MY_DIM>;
+
+#else
+
+#error requested a grid that does not exist
+
+#endif

src/hdf5-writer.hh

+#ifndef SRC_HDF5_WRITER_HH
+#define SRC_HDF5_WRITER_HH
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/hdf5/file.hh>
+
+#include "hdf5/frictionalboundary-writer.hh"
+#include "hdf5/iteration-writer.hh"
+#include "hdf5/patchinfo-writer.hh"
+#include "hdf5/surface-writer.hh"
+#include "hdf5/time-writer.hh"
+
+template <class ProgramState, class VertexBasis, class GridView>
+class HDF5Writer {
+private:
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+  using LocalVector = typename Vector::block_type;
+
+public:
+  HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
+             VertexBasis const &vertexBasis, Patch const &surface,
+             Patch const &frictionalBoundary,
+             ConvexPolyhedron<LocalVector> const &weakPatch)
+      : file_(file),
+        iterationWriter_(file_),
+        timeWriter_(file_),
+#if MY_DIM == 3
+        patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
+#endif
+        surfaceWriter_(file_, vertexCoordinates, surface),
+        frictionalBoundaryWriter_(file_, vertexCoordinates,
+                                  frictionalBoundary) {
+  }
+
+  template <class Friction>
+  void reportSolution(ProgramState const &programState,
+                      // for the friction coefficient
+                      Friction &friction) {
+    timeWriter_.write(programState);
+#if MY_DIM == 3
+    patchInfoWriter_.write(programState);
+#endif
+    surfaceWriter_.write(programState);
+    frictionalBoundaryWriter_.write(programState, friction);
+  }
+
+  void reportIterations(ProgramState const &programState,
+                        IterationRegister const &iterationCount) {
+    iterationWriter_.write(programState.timeStep, iterationCount);
+  }
+
+private:
+  HDF5::Grouplike &file_;
+
+  IterationWriter iterationWriter_;
+  TimeWriter<ProgramState> timeWriter_;
+#if MY_DIM == 3
+  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
+#endif
+  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
+  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
+};
+#endif

src/hdf5/frictionalboundary-writer.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/hdf5/singletonwriter.hh>
+
+#include "frictionalboundary-writer.hh"
+#include "restrict.hh"
+
+template <class ProgramState, class GridView>
+FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
+    HDF5::Grouplike &file, Vector const &vertexCoordinates,
+    Patch const &frictionalBoundary)
+    : group_(file, "frictionalBoundary"),
+      frictionalBoundary_(frictionalBoundary),
+      frictionalBoundaryDisplacementWriter_(group_, "displacement",
+                                            frictionalBoundary.numVertices(),
+                                            Vector::block_type::dimension),
+      frictionalBoundaryVelocityWriter_(group_, "velocity",
+                                        frictionalBoundary.numVertices(),
+                                        Vector::block_type::dimension),
+      frictionalBoundaryStateWriter_(group_, "state",
+                                     frictionalBoundary.numVertices()),
+      frictionalBoundaryCoefficientWriter_(group_, "coefficient",
+                                           frictionalBoundary.numVertices()) {
+  auto const frictionalBoundaryCoordinates =
+      restrictToSurface(vertexCoordinates, frictionalBoundary);
+  HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
+      group_, "coordinates", frictionalBoundaryCoordinates.size(),
+      Vector::block_type::dimension);
+  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
+}
+
+template <class ProgramState, class GridView>
+template <class Friction>
+void FrictionalBoundaryWriter<ProgramState, GridView>::write(
+    ProgramState const &programState, Friction &friction) {
+  auto const frictionalBoundaryDisplacements =
+      restrictToSurface(programState.u, frictionalBoundary_);
+  addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
+           frictionalBoundaryDisplacements);
+
+  auto const frictionalBoundaryVelocities =
+      restrictToSurface(programState.v, frictionalBoundary_);
+  addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
+           frictionalBoundaryVelocities);
+
+  auto const frictionalBoundaryStates =
+      restrictToSurface(programState.alpha, frictionalBoundary_);
+  addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
+           frictionalBoundaryStates);
+
+  friction.updateAlpha(programState.alpha);
+  auto const c = friction.coefficientOfFriction(programState.v);
+  auto const frictionalBoundaryCoefficient =
+      restrictToSurface(c, frictionalBoundary_);
+  addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
+           frictionalBoundaryCoefficient);
+}
+
+#include "frictionalboundary-writer_tmpl.cc"

src/hdf5/frictionalboundary-writer.hh

+#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include <dune/fufem/hdf5/sequenceio.hh>
+
+template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
+  using ScalarVector = typename ProgramState::ScalarVector;
+  using Vector = typename ProgramState::Vector;
+  using Patch = BoundaryPatch<GridView>;
+
+public:
+  FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
+                           Patch const &frictionalBoundary);
+
+  template <class Friction>
+  void write(ProgramState const &programState, Friction &friction);
+
+private:
+  HDF5::Group group_;
+
+  Patch const &frictionalBoundary_;
+
+  HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
+  HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
+  HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
+};
+#endif

src/hdf5/frictionalboundary-writer_tmpl.cc

+#include "../explicitvectors.hh"
+#include "../explicitgrid.hh"
+
+#include "../program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+using MyFriction = GlobalFriction<Matrix, Vector>;
+
+template class FrictionalBoundaryWriter<MyProgramState, GridView>;
+template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
+    MyProgramState const &programState, MyFriction &friction);

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