src/hdf5/surface-writer.hh → src/io/hdf5/surface-writer.hh

@@ -12,11 +12,13 @@ template <class ProgramState, class GridView> class SurfaceWriter {
 
 public:
   SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                Patch const &surface);
+                Patch const &surface, size_t const id);
 
   void write(ProgramState const &programState);
 
 private:
+  size_t const id_;
+
   HDF5::Group group_;
 
   Patch const &surface_;

src/io/hdf5/surface-writer_tmpl.cc

+#include "../../explicitvectors.hh"
+#include "../../explicitgrid.hh"
+
+#include "../../data-structures/program_state.hh"
+
+using MyProgramState = ProgramState<Vector, ScalarVector>;
+
+template class SurfaceWriter<MyProgramState, DeformedGrid::LeafGridView>;

src/hdf5/time-writer_tmpl.cc → src/io/hdf5/time-writer_tmpl.cc

-#include "../explicitvectors.hh"
+#include "../../explicitvectors.hh"
 
-#include "../program_state.hh"
+#include "../../data-structures/program_state.hh"
 
 using MyProgramState = ProgramState<Vector, ScalarVector>;
 

src/uniform-grid-writer.cc → src/io/uniform-grid-writer.cc

@@ -11,10 +11,13 @@
 // #include <dune/common/parametertreeparser.hh>
 
 #include "assemblers.hh"
-#include "diameter.hh"
 #include "gridselector.hh"
+
+#include "io/vtk.hh"
+
 #include "one-body-problem-data/mygrid.hh"
-#include "vtk.hh"
+
+#include "utils/diameter.hh"
 
 size_t const dims = MY_DIM;
 size_t const refinements = 5; // FIXME?

src/vtk.cc → src/io/vtk.cc

@@ -9,50 +9,81 @@
 #include "vtk.hh"
 
 template <class VertexBasis, class CellBasis>
-MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(
-    CellBasis const &_cellBasis, VertexBasis const &_vertexBasis,
-    std::string _prefix)
-    : cellBasis(_cellBasis), vertexBasis(_vertexBasis), prefix(_prefix) {}
+BodyVTKWriter<VertexBasis, CellBasis>::BodyVTKWriter(
+    CellBasis const & cellBasis, VertexBasis const & vertexBasis,
+    std::string prefix)
+    : cellBasis_(cellBasis), vertexBasis_(vertexBasis), prefix_(prefix) {}
 
 template <class VertexBasis, class CellBasis>
 template <class Vector, class ScalarVector>
-void MyVTKWriter<VertexBasis, CellBasis>::write(
+void BodyVTKWriter<VertexBasis, CellBasis>::write(
     size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha,
     ScalarVector const &stress) const {
   Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis.getGridView());
+      vertexBasis_.getGridView());
 
   auto const displacementPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis, u, "displacement");
+          vertexBasis_, u, "displacement");
   writer.addVertexData(displacementPointer);
 
   auto const velocityPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, Vector> const>(
-          vertexBasis, v, "velocity");
+          vertexBasis_, v, "velocity");
   writer.addVertexData(velocityPointer);
 
   auto const AlphaPointer =
       std::make_shared<VTKBasisGridFunction<VertexBasis, ScalarVector> const>(
-          vertexBasis, alpha, "Alpha");
+          vertexBasis_, alpha, "Alpha");
   writer.addVertexData(AlphaPointer);
 
   auto const stressPointer =
       std::make_shared<VTKBasisGridFunction<CellBasis, ScalarVector> const>(
-          cellBasis, stress, "stress");
+          cellBasis_, stress, "stress");
   writer.addCellData(stressPointer);
 
-  std::string const filename = prefix + std::to_string(record);
+  std::string const filename = prefix_ + "_" + std::to_string(record);
   writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
 template <class VertexBasis, class CellBasis>
-void MyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
+void BodyVTKWriter<VertexBasis, CellBasis>::writeGrid() const {
   Dune::VTKWriter<typename VertexBasis::GridView> writer(
-      vertexBasis.getGridView());
+      vertexBasis_.getGridView());
 
-  std::string const filename = prefix + "_grid";
+  std::string const filename = prefix_ + "_grid";
   writer.write(filename.c_str(), Dune::VTK::appendedraw);
 }
 
+template <class VertexBasis, class CellBasis>
+MyVTKWriter<VertexBasis, CellBasis>::MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
+  std::string prefix): bodyVTKWriters_(vertexBases.size()), prefix_(prefix) {
+
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i] = new BodyVTKWriter<VertexBasis, CellBasis>(*cellBases[i], *vertexBases[i], prefix_+std::to_string(i));
+  }
+}
+
+template <class VertexBasis, class CellBasis>
+MyVTKWriter<VertexBasis, CellBasis>::~MyVTKWriter() {
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    delete bodyVTKWriters_[i];
+  }
+}
+
+template <class VertexBasis, class CellBasis>
+template <class Vector, class ScalarVector>
+void MyVTKWriter<VertexBasis, CellBasis>::write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+                        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const {
+
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i]->write(record, u[i], v[i], alpha[i], stress[i]);
+  }
+}
+template <class VertexBasis, class CellBasis>
+void MyVTKWriter<VertexBasis, CellBasis>::writeGrids() const {
+  for (size_t i=0; i<bodyVTKWriters_.size(); i++) {
+    bodyVTKWriters_[i]->writeGrid();
+  }
+}
 #include "vtk_tmpl.cc"

src/vtk.hh → src/io/vtk.hh

@@ -3,13 +3,14 @@
 
 #include <string>
 
-template <class VertexBasis, class CellBasis> class MyVTKWriter {
-  CellBasis const &cellBasis;
-  VertexBasis const &vertexBasis;
-  std::string const prefix;
+template <class VertexBasis, class CellBasis> class BodyVTKWriter {
+private:
+  CellBasis const &cellBasis_;
+  VertexBasis const &vertexBasis_;
+  std::string const prefix_;
 
 public:
-  MyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis,
+  BodyVTKWriter(CellBasis const &cellBasis, VertexBasis const &vertexBasis,
               std::string prefix);
 
   template <class Vector, class ScalarVector>
@@ -18,4 +19,22 @@ template <class VertexBasis, class CellBasis> class MyVTKWriter {
 
   void writeGrid() const;
 };
+
+template <class VertexBasis, class CellBasis> class MyVTKWriter {
+private:
+  std::vector<BodyVTKWriter<VertexBasis, CellBasis>* > bodyVTKWriters_;
+  std::string const prefix_;
+
+public:
+  MyVTKWriter(const std::vector<const CellBasis* >& cellBases, const std::vector<const VertexBasis* >& vertexBases,
+              std::string prefix);
+
+  ~MyVTKWriter();
+
+  template <class Vector, class ScalarVector>
+  void write(size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+             const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;
+
+  void writeGrids() const;
+};
 #endif

src/vtk_tmpl.cc → src/io/vtk_tmpl.cc

@@ -2,17 +2,17 @@
 #error MY_DIM unset
 #endif
 
-#include "explicitgrid.hh"
-#include "explicitvectors.hh"
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
 
 #include <dune/fufem/functionspacebases/p0basis.hh>
 #include <dune/fufem/functionspacebases/p1nodalbasis.hh>
 
-using MyP0Basis = P0Basis<GridView, double>;
-using P1Basis = P1NodalBasis<GridView, double>;
+using MyP0Basis = P0Basis<DeformedGrid::LeafGridView, double>;
+using P1Basis = P1NodalBasis<DeformedGrid::LeafGridView, double>;
 
 template class MyVTKWriter<P1Basis, MyP0Basis>;
 
 template void MyVTKWriter<P1Basis, MyP0Basis>::write<Vector, ScalarVector>(
-    size_t record, Vector const &u, Vector const &v, ScalarVector const &alpha,
-    ScalarVector const &stress) const;
+    size_t record, const std::vector<Vector>& u, const std::vector<Vector>& v,
+        const std::vector<ScalarVector>& alpha, const std::vector<ScalarVector>& stress) const;

src/matrices.hh

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

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

+# -*- mode:conf -*-
+[boundary.friction]
+smallestDiameter = 0.5  # 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/multi-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/multi-body-problem-data/bc.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_BC_HH
+
+#include <dune/common/function.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/multi-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/multi-body-problem-data/grid/cube.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cube.hh"
+
+template <int dim>
+Cube<dim>::Cube(bool invariant = false) : invariant_(invariant) {
+    if (invariant_) {
+        nChildren_ = 1;
+    } else {
+        nChildren_ = std::pow(2, dim);
+    }
+    children_.resize(nChildren_, nullptr);
+}
+
+template <int dim>
+Cube<dim>::Cube(const Vector& A, const Vector& B, bool invariant = false) : invariant_(invariant) {
+    setCorners(A, B);
+
+    if (invariant_) {
+        nChildren_ = 1;
+    } else {
+        nChildren_ = std::pow(2, dim);
+    }
+    children_.resize(nChildren_, nullptr);
+}
+
+// generate a combination of unit basis vectors from binary number
+template <int dim>
+void Cube<dim>::parseBinary(int binary, Vector& shift) const {
+
+    shift = 0;
+    for (int d=0; d<dim; d++) {
+        // yields true, if d-th digit of binary is a 1
+        if (binary & std::pow(2, d)) {
+            // use d-th basis vector
+            shift[d] = 1;
+        }
+    }
+}
+
+// constructs child cubes and sets children_
+template <int dim>
+void Cube<dim>::split() {
+    if (invariant_) {
+        children_[0] = this;
+    } else {
+        const int nNewCubes = std::pow(2, dim);
+        newCubes.resize(nNewCubes);
+
+        Vector midPoint = A_;
+        midPoint += 1/2*(B_ - A_);
+
+        const double h = std::pow(2.0, std::round(std::log(midPoint[0]-A_[0])/std::log(2)));
+
+        newCubes[0] = new Cube(A_, midPoint, false);
+        newCubes[nNewCubes-1] = new Cube(midPoint, B_, true);
+
+        for (int i=1; i<nNewCubes-1; i++) {
+            Vector shift;
+            parseBinary(i, shift);
+            shift *= h;
+
+            newCubes[i] = new Cube(A_+shift, midPoint+shift);
+        }
+    } else {
+        children_.resize(1);
+    }
+}
+
+#include "cube_tmpl.cc"

src/multi-body-problem-data/grid/cube.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBE_HH
+
+
+// works in any space dimension
+template <int dim>
+class Cube {
+private:
+    using Vector = Dune::FieldVector<double, dim>;
+
+    // generate a combination of unit basis vectors from binary number
+    void parseBinary(int binary, Vector& shift) const;
+
+public:
+    Cube(bool invariant = false);
+    Cube(const Vector& A, const Vector& B, bool invariant = false);
+
+    void setCorners(const Vector& A, const Vector& B) {
+        A_ = A;
+        B_ = B;
+    }
+
+    void setParent(Cube* parent) {
+        parent_ = parent;
+    }
+
+    const std::vector<std::shared_ptr<Cube>>& children() const {
+        return children_;
+    }
+
+    // constructs child cubes and sets children_
+    void split();
+
+private:
+    using Vector = Dune::FieldVector<double, dim>;
+
+    Vector A_; // two corners that are diagonal define dim-cube in any space dimension
+    Vector B_;
+
+    const bool invariant_; // flag to set if Cube can be split()
+
+    std::shared_ptr<Cube<dim>> parent_ = nullptr;
+    std::vector<std::shared_ptr<Cube<dim>>> children_;
+    int nChildren_;
+};
+
+#endif
+
+

src/multi-body-problem-data/grid/cube_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+template class Cube<MY_DIM>;

src/multi-body-problem-data/grid/cubefaces.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "cubefaces.hh"
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
+                                    Vector const &c) {
+  return isClose2((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBoxed(Vector const &v1, Vector const &v2,
+                                Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 * lengthScale_ > x[0] or
+      x[0] > minmax0.second + 1e-14 * lengthScale_)
+    return false;
+  if (minmax1.first - 1e-14 * lengthScale_ > x[1] or
+      x[1] > minmax1.second + 1e-14 * lengthScale_)
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool CubeFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+CubeFaces<GridView>::CubeFaces(const GridView& gridView, const Cube& cube, double lengthScale)
+      :
+  #if MY_DIM == 3
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+        front(gridView),
+        back(gridView),
+  #else
+        lower(gridView),
+        right(gridView),
+        upper(gridView),
+        left(gridView),
+  #endif
+        cube_(cube),
+        lengthScale_(lengthScale)
+  {
+    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.A, cuboidGeometry.B, in.geometry().center());
+    });
+
+    right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.B, cuboidGeometry.C, in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.D, cuboidGeometry.C, in.geometry().center());
+    });
+
+    left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.A, cuboidGeometry.D, in.geometry().center());
+    });
+
+  #if MY_DIM == 3
+    front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+    });
+
+    back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(-cuboidGeometry.depth / 2.0, in.geometry().center()[2]);
+    });
+  #endif
+  }
+
+#include "cubefaces_tmpl.cc"

src/multi-body-problem-data/grid/cubefaces.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEFACES_HH
+
+#include <dune/fufem/boundarypatch.hh>
+
+#include "cube.hh"
+
+template <class GridView>
+class CubeFaces {
+private:
+    using Cube = Cube<GridView::dimensionworld>;
+
+    bool isClose(double a, double b) {
+      return std::abs(a - b) < 1e-14 * lengthScale_;
+    }
+
+    bool isClose2(double a, double b) {
+      return std::abs(a - b) <
+             1e-14 * lengthScale_ * lengthScale_;
+    }
+
+    template <class Vector>
+    bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
+
+    template <class Vector>
+    bool xyCollinear(Vector const &a, Vector const &b, Vector const &c);
+
+    template <class Vector>
+    bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x);
+
+public:
+    CubeFaces(const GridView& gridView, const Cube<GridView::dimensionworld>& cube, double lengthScale);
+
+    const BoundaryPatch<GridView>& lower() const {
+        return lower_;
+    }
+
+    const BoundaryPatch<GridView>& right() const {
+        return right_;
+    }
+
+    const BoundaryPatch<GridView>& upper() const {
+        return upper_;
+    }
+
+    const BoundaryPatch<GridView>& left() const {
+        return left_;
+    }
+
+#if MY_DIM == 3
+    const BoundaryPatch<GridView>& front() const {
+        return front_;
+    }
+    const BoundaryPatch<GridView>& back() const {
+        return back_;
+    }
+#endif
+
+private:
+  BoundaryPatch<GridView> lower_;
+  BoundaryPatch<GridView> right_;
+  BoundaryPatch<GridView> upper_;
+  BoundaryPatch<GridView> left_;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> front_;
+  BoundaryPatch<GridView> back_;
+#endif
+
+  const Cube& cube_;
+  const double lengthScale_;
+};
+
+#endif
+
+

src/multi-body-problem-data/grid/cubefaces_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+
+template struct CubeFaces<DeformedGrid::LeafGridView>;
+template struct CubeFaces<DeformedGrid::LevelGridView>;