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/midpoint.hh

+#ifndef SRC_MIDPOINT_HH
+#define SRC_MIDPOINT_HH
+
+#include <dune/matrix-vector/axpy.hh>
+
+template <class Vector> Vector midPoint(Vector const &x, Vector const &y) {
+  Vector ret(0);
+  Dune::MatrixVector::addProduct(ret, 0.5, x);
+  Dune::MatrixVector::addProduct(ret, 0.5, y);
+  return ret;
+}
+#endif

src/multi-body-problem-data/mybody.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYBODY_HH
+#define SRC_MULTI_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 "cuboidgeometry.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/multi-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

src/multi-body-problem-data/mygrids.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "mygrids.hh"
+#include "midpoint.hh"
+#include "../diameter.hh"
+
+#if MY_DIM == 3
+SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
+#endif
+
+// back-to-front, front-to-back, front-to-back
+void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
+                                        unsigned int W) {
+#if MY_DIM == 3
+  unsigned int const U2 = U + shift_;
+  unsigned int const V2 = V + shift_;
+  unsigned int const W2 = W + shift_;
+
+  simplices_.push_back({ U, V, W, U2 });
+  simplices_.push_back({ V, V2, W2, U2 });
+  simplices_.push_back({ W, W2, U2, V });
+#else
+  simplices_.push_back({ U, V, W });
+#endif
+}
+
+// back-to-front, back-to-front, front-to-back
+void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
+                                        unsigned int W) {
+#if MY_DIM == 3
+  unsigned int const U2 = U + shift_;
+  unsigned int const V2 = V + shift_;
+  unsigned int const W2 = W + shift_;
+
+  simplices_.push_back({ U, V, W, U2 });
+  simplices_.push_back({ V, V2, W, U2 });
+  simplices_.push_back({ V2, W, U2, W2 });
+#else
+  simplices_.push_back({ U, V, W });
+#endif
+}
+
+auto SimplexManager::getSimplices() -> SimplexList const & {
+  return simplices_;
+}
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_) :
+  cuboidGeometries(cuboidGeometries_)
+{
+  size_t const gridCount = cuboidGeometries.size();
+  grids.resize(gridCount);
+  gridFactories.resize(gridCount);
+
+  for (size_t idx=0; idx<grids.size(); idx++) {
+    const auto& cuboidGeometry = *cuboidGeometries[idx];
+
+    const auto& A = cuboidGeometry.A;
+    const auto& B = cuboidGeometry.B;
+    const auto& C = cuboidGeometry.C;
+    const auto& D = cuboidGeometry.D;
+
+    unsigned int const vc = 4;
+
+#if MY_DIM == 3
+    Dune::FieldMatrix<double, 2 * vc, MY_DIM> vertices;
+#else
+    Dune::FieldMatrix<double, vc, MY_DIM> vertices;
+#endif
+    for (size_t i = 0; i < 2; ++i) {
+#if MY_DIM == 3
+      size_t numXYplanes = 2;
+#else
+      size_t numXYplanes = 1;
+#endif
+      size_t k = 0;
+      for (size_t j = 1; j <= numXYplanes; ++j) {
+        vertices[k++][i] = A[i];
+        vertices[k++][i] = B[i];
+        vertices[k++][i] = C[i];
+        vertices[k++][i] = D[i];
+        assert(k == j * vc);
+      }
+    }
+
+#if MY_DIM == 3
+    for (size_t k = 0; k < vc; ++k) {
+      vertices[k][2] = -cuboidGeometry.depth / 2.0;
+      vertices[k + vc][2] = cuboidGeometry.depth / 2.0;
+    }
+#endif
+
+    for (size_t i = 0; i < vertices.N(); ++i)
+      gridFactories[idx].insertVertex(vertices[i]);
+
+    Dune::GeometryType cell;
+#if MY_DIM == 3
+    cell = Dune::GeometryTypes::tetrahedron;
+#else
+    cell = Dune::GeometryTypes::triangle;
+#endif
+
+#if MY_DIM == 3
+    SimplexManager sm(vc);
+#else
+    SimplexManager sm;
+#endif
+    sm.addFromVerticesFFB(1, 2, 0);
+    sm.addFromVerticesFFB(2, 3, 0);
+    auto const &simplices = sm.getSimplices();
+
+    // sanity-check choices of simplices
+    for (size_t i = 0; i < simplices.size(); ++i) {
+      Dune::FieldMatrix<double, MY_DIM, MY_DIM> check;
+      for (size_t j = 0; j < MY_DIM; ++j)
+        check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
+      assert(check.determinant() > 0);
+      gridFactories[idx].insertElement(cell, simplices[i]);
+    }
+
+    grids[idx] = std::shared_ptr<Grid>(gridFactories[idx].createGrid());
+  }
+}
+
+template <class Grid>
+std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() {
+  return grids;
+}
+
+template <class Grid>
+template <class GridView>
+MyFaces<GridView> GridsConstructor<Grid>::constructFaces(
+    GridView const &gridView, CuboidGeometry const &cuboidGeometry) {
+  return MyFaces<GridView>(gridView, cuboidGeometry);
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<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 MyFaces<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 * cuboidGeometry.lengthScale > x[0] or
+      x[0] > minmax0.second + 1e-14 * cuboidGeometry.lengthScale)
+    return false;
+  if (minmax1.first - 1e-14 * cuboidGeometry.lengthScale > x[1] or
+      x[1] > minmax1.second + 1e-14 * cuboidGeometry.lengthScale)
+    return false;
+
+  return true;
+}
+
+template <class GridView>
+template <class Vector>
+bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
+                                  Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+template <class GridView>
+MyFaces<GridView>::MyFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry_)
+      :
+  #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
+        cuboidGeometry(cuboidGeometry_)
+  {
+    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
+  }
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+  return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+}
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter, double lengthScale) {
+  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 * lengthScale);
+      auto const admissibleDiameter =
+          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);
+
+      if (diameter(geometry) <= admissibleDiameter)
+        continue;
+
+      needRefine = true;
+      grid.mark(1, e);
+    }
+    if (!needRefine)
+      break;
+
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+  }
+}
+
+#include "mygrids_tmpl.cc"

src/multi-body-problem-data/mygrids.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_MYGRIDS_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_MYGRIDS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include "cuboidgeometry.hh"
+
+template <class GridView> struct MyFaces {
+  BoundaryPatch<GridView> lower;
+  BoundaryPatch<GridView> right;
+  BoundaryPatch<GridView> upper;
+  BoundaryPatch<GridView> left;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> front;
+  BoundaryPatch<GridView> back;
+#endif
+
+  MyFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry_);
+
+private:
+  CuboidGeometry const &cuboidGeometry;
+
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale;
+  };
+
+  bool isClose2(double a, double b) {
+    return std::abs(a - b) <
+           1e-14 * cuboidGeometry.lengthScale * cuboidGeometry.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);
+};
+
+class SimplexManager {
+public:
+  using SimplexList = std::vector<std::vector<unsigned int>>;
+
+#if MY_DIM == 3
+  SimplexManager(unsigned int shift);
+#endif
+
+  void addFromVerticesFBB(unsigned int U, unsigned int V, unsigned int W);
+  void addFromVerticesFFB(unsigned int U, unsigned int V, unsigned int W);
+
+  SimplexList const &getSimplices();
+
+private:
+  SimplexList simplices_;
+
+#if MY_DIM == 3
+  unsigned int const shift_;
+#endif
+};
+
+template <class Grid> class GridsConstructor {
+public:
+  GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+
+  std::vector<std::shared_ptr<Grid>>& getGrids();
+
+  template <class GridView>
+  MyFaces<GridView> constructFaces(GridView const &gridView, CuboidGeometry const &cuboidGeometry);
+
+private:
+  std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+  std::vector<Dune::GridFactory<Grid>> gridFactories;
+  std::vector<std::shared_ptr<Grid>> grids;
+};
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale);
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter, double lengthScale);
+
+#endif
+
+

src/multi-body-problem-data/mygrids_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+#include "cuboidgeometry.hh"
+
+template class GridsConstructor<Grid>;
+
+template struct MyFaces<LeafGridView>;
+template struct MyFaces<LevelGridView>;
+
+template MyFaces<LeafGridView> GridsConstructor<Grid>::constructFaces(
+    LeafGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+
+template MyFaces<LevelGridView> GridsConstructor<Grid>::constructFaces(
+    LevelGridView const &gridView, CuboidGeometry const &CuboidGeometry_);
+
+template void refine<Grid, LocalVector>(
+    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+    double smallestDiameter, double lengthScale);

src/multi-body-problem-data/patchfunction.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+class PatchFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  using Polyhedron = ConvexPolyhedron<Dune::FieldVector<double, MY_DIM>>;
+
+  double const v1_;
+  double const v2_;
+  Polyhedron const &segment_;
+
+public:
+  PatchFunction(double v1, double v2, Polyhedron const &segment)
+      : v1_(v1), v2_(v2), segment_(segment) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+                                                                      : v1_;
+  }
+};
+
+#endif

src/multi-body-problem-data/segmented-function.hh

+#ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include "cuboidgeometry.hh"
+
+class SegmentedFunction
+    : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
+                                   Dune::FieldVector<double, 1>> {
+private:
+  bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
+                 Dune::FieldVector<double, MY_DIM> const &y,
+                 Dune::FieldVector<double, MY_DIM> const &z) const {
+    return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
+  };
+  bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
+    return liesBelow(MyGeometry::K, MyGeometry::M, z);
+  };
+
+  double const _v1;
+  double const _v2;
+
+public:
+  SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
+
+  void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
+                Dune::FieldVector<double, 1> &y) const {
+    y = insideRegion2(x) ? _v2 : _v1;
+  }
+};
+#endif

src/multi-body-problem-data/weakpatch.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_WEAKPATCH_HH
+
+#include "cuboidgeometry.hh"
+
+template <class LocalVector>
+ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset, CuboidGeometry const &cuboidGeometry) {
+  ConvexPolyhedron<LocalVector> weakPatch;
+#if MY_DIM == 3
+  weakPatch.vertices.resize(4);
+  weakPatch.vertices[0] = weakPatch.vertices[2] = cuboidGeometry.X;
+  weakPatch.vertices[1] = weakPatch.vertices[3] = cuboidGeometry.Y;
+  for (size_t k = 0; k < 2; ++k) {
+    weakPatch.vertices[k][2] = -cuboidGeometry.depth / 2.0;
+    weakPatch.vertices[k + 2][2] = cuboidGeometry.depth / 2.0;
+  }
+  switch (parset.get<Config::PatchType>("patchType")) {
+    case Config::Rectangular:
+      break;
+    case Config::Trapezoidal:
+      weakPatch.vertices[1][0] += 0.05 * cuboidGeometry.lengthScale;
+      weakPatch.vertices[3][0] -= 0.05 * cuboidGeometry.lengthScale;
+      break;
+    default:
+      assert(false);
+  }
+#else
+  weakPatch.vertices.resize(2);
+  weakPatch.vertices[0] = cuboidGeometry.X;
+  weakPatch.vertices[1] = cuboidGeometry.Y;
+#endif
+  return weakPatch;
+};
+#endif

src/multi-body-problem.cfg

+# -*- mode:conf -*-
+gravity         = 9.81  # [m/s^2]
+
+[io]
+data.write      = true
+printProgress   = false
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = true
+vtk.write       = false
+
+[problem]
+finalTime       = 1000  # [s]
+bodyCount       = 2
+
+[body]
+bulkModulus     = 0.5e5 # [Pa]
+poissonRatio    = 0.3   # [1]
+[body.elastic]
+density         = 900   # [kg/m^3]
+shearViscosity  = 1e3   # [Pas]
+bulkViscosity   = 1e3   # [Pas]
+[body.viscoelastic]
+density         = 1000  # [kg/m^3]
+shearViscosity  = 1e4   # [Pas]
+bulkViscosity   = 1e4   # [Pas]
+
+[boundary.friction]
+C               = 10    # [Pa]
+mu0             = 0.7   # [ ]
+V0              = 5e-5  # [m/s]
+L               = 2.25e-5 # [m]
+initialAlpha    = 0     # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Regularised
+[boundary.friction.weakening]
+a               = 0.002 # [ ]
+b               = 0.017 # [ ]
+[boundary.friction.strengthening]
+a               = 0.020 # [ ]
+b               = 0.005 # [ ]
+
+[timeSteps]
+scheme = newmark
+
+[u0.solver]
+maximumIterations = 100000
+verbosity         = quiet
+
+[a0.solver]
+maximumIterations = 100000
+verbosity         = quiet
+
+[v.solver]
+maximumIterations = 100000
+verbosity         = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda            = 0.5
+
+[solver.tnnmg.linear]
+maxiumumIterations = 100000
+pre                = 3
+cycle              = 1  # 1 = V, 2 = W, etc.
+post               = 3
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0

src/one-body-problem.cc

@@ -31,9 +31,14 @@
 #include <dune/fufem/formatstring.hh>
 
 #include <dune/solvers/norms/energynorm.hh>
+
+
+/*
 #include <dune/solvers/solvers/loopsolver.hh>
 #include <dune/solvers/solvers/solver.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
+*/
+
 
 #include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/myblockproblem.hh>
@@ -62,13 +67,18 @@
 #include "time-stepping/updaters.hh"
 #include "vtk.hh"
 
+// for getcwd
+#include <unistd.h>
+
+#define USE_OLD_TNNMG
+
 size_t const dims = MY_DIM;
 
 Dune::ParameterTree getParameters(int argc, char *argv[]) {
   Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
+  Dune::ParameterTreeParser::readINITree("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem.cfg", parset);
   Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("one-body-problem-%dD.cfg", dims), parset);
+      Dune::Fufem::formatString("/home/mi/podlesny/software/dune/dune-tectonic/src/one-body-problem-%dD.cfg", dims), parset);
   Dune::ParameterTreeParser::readOptions(argc, argv, parset);
   return parset;
 }
@@ -79,6 +89,14 @@ void handleSignal(int signum) { terminationRequested = true; }
 int main(int argc, char *argv[]) {
   try {
     Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
     auto const parset = getParameters(argc, argv);
 
     MyGeometry::render();
@@ -140,6 +158,7 @@ int main(int argc, char *argv[]) {
 #endif
     }
 
+
     // Set up functions for time-dependent boundary conditions
     using Function = Dune::VirtualFunction<double, double>;
     Function const &velocityDirichletFunction = VelocityDirichletCondition();
@@ -184,8 +203,10 @@ int main(int argc, char *argv[]) {
     auto const writeData = parset.get<bool>("io.data.write");
     bool const handleRestarts = writeRestarts or firstRestart > 0;
 
+    /*
     auto dataFile =
         writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
+
     auto restartFile = handleRestarts
                            ? std::make_unique<HDF5::File>(
                                  "restarts.h5",
@@ -214,7 +235,7 @@ int main(int argc, char *argv[]) {
     Vector vertexCoordinates(leafVertexCount);
     {
       Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
+          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView, Dune::mcmgVertexLayout());
       for (auto &&v : vertices(leafView))
         vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
     }
@@ -229,6 +250,7 @@ int main(int argc, char *argv[]) {
     MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
         myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
 
+
     IterationRegister iterationCount;
     auto const report = [&](bool initial = false) {
       if (writeData) {
@@ -261,6 +283,7 @@ int main(int argc, char *argv[]) {
     };
     report(true);
 
+
     // Set up TNNMG solver
     using NonlinearFactory = SolverFactory<
         dims,
@@ -321,6 +344,8 @@ int main(int argc, char *argv[]) {
         break;
       }
     }
+
+    */
   } catch (Dune::Exception &e) {
     Dune::derr << "Dune reported error: " << e << std::endl;
   } catch (std::exception &e) {

src/program_state.hh

@@ -3,10 +3,14 @@
 
 #include <dune/common/parametertree.hh>
 
+#include <dune/matrix-vector/axpy.hh>
+
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
 #include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
 #include <dune/tectonic/body.hh>
 
 #include "assemblers.hh"
@@ -18,32 +22,45 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
   using Vector = VectorTEMPLATE;
   using ScalarVector = ScalarVectorTEMPLATE;
 
-  ProgramState(int leafVertexCount)
-      : u(leafVertexCount),
-        v(leafVertexCount),
-        a(leafVertexCount),
-        alpha(leafVertexCount),
-        weightedNormalStress(leafVertexCount) {}
+  ProgramState(const std::vector<int>& leafVertexCounts)
+    : bodyCount(leafVertexCounts.size()) {
+    u.resize(bodyCount);
+    v.resize(bodyCount);
+    a.resize(bodyCount);
+    alpha.resize(bodyCount);
+    weightedNormalStress.resize(bodyCount);
+
+    for (size_t i=0; i<bodyCount; i++) {
+        size_t leafVertexCount = leafVertexCounts[i];
+        u[i].resize(leafVertexCount);
+        v[i].resize(leafVertexCount);
+        a[i].resize(leafVertexCount);
+        alpha[i].resize(leafVertexCount);
+        weightedNormalStress[i].resize(leafVertexCount);
+    }
+  }
 
   // Set up initial conditions
   template <class Matrix, class GridView>
   void setupInitialConditions(
-      Dune::ParameterTree const &parset,
-      std::function<void(double, Vector &)> externalForces,
-      Matrices<Matrix> const matrices,
-      MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
-      Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
-      Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      Body<Vector::block_type::dimension> const &body) {
+      const Dune::ParameterTree& parset,
+      const Dune::Contact::NBodyAssembler<typename GridView::Grid, Vector>& nBodyAssembler,
+      std::vector<std::function<void(double, Vector &)>> externalForces,
+      const Matrices<Matrix>& matrices,
+      const std::vector<std::shared_ptr<MyAssembler<GridView, Vector::block_type::dimension>>>& assemblers,
+      const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& dirichletNodes,
+      const std::vector<Dune::BitSetVector<Vector::block_type::dimension>>& noNodes,
+      const std::vector<BoundaryPatch<GridView>>& frictionalBoundaries,
+      const Body<Vector::block_type::dimension>& body) {
 
     using LocalVector = typename Vector::block_type;
     using LocalMatrix = typename Matrix::block_type;
     auto constexpr dims = LocalVector::dimension;
 
+    /*
     // Solving a linear problem with a multigrid solver
     auto const solveLinearProblem = [&](
-        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
+        Dune::BitSetVector<dims> const &_dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrix,
         Vector const &_rhs, Vector &_x,
         Dune::ParameterTree const &_localParset) {
 
@@ -54,7 +71,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
 
       LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
-                            myAssembler.gridView.grid(), _dirichletNodes);
+                            assemblers.gridView.grid(), _dirichletNodes);
 
       typename LinearFactory::ConvexProblem convexProblem(
           1.0, _matrix, zeroNonlinearity, _rhs, _x);
@@ -62,7 +79,11 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
       auto multigridStep = factory.getStep();
       multigridStep->setProblem(_x, problem);
+      //multigridStep->setProblem(_x);
       EnergyNorm<Matrix, Vector> const norm(_matrix);
+
+
+
       LoopSolver<Vector> solver(
           multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
           _localParset.get<double>("tolerance"), &norm,
@@ -71,17 +92,28 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
       solver.preprocess();
       solver.solve();
+
+      Vector totalX = multigridStep->getSol();
+
+      // cleanup
+      delete(multigridStep);
+
+      nBodyAssembler.postprocess(totalX, x);
     };
+    */
 
     timeStep = 0;
     relativeTime = 0.0;
     relativeTau = 1e-6;
 
-    Vector ell0(u.size());
-    externalForces(relativeTime, ell0);
+    std::vector<Vector> ell0(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+      // Initial velocity
+      v[i] = 0.0;
 
-    // Initial velocity
-    v = 0.0;
+      ell0[i].resize(u[i].size());
+      externalForces[i](relativeTime, ell0[i]);
+    }
 
     // Initial displacement: Start from a situation of minimal stress,
     // which is automatically attained in the case [v = 0 = a].
@@ -91,29 +123,35 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
     // Initial acceleration: Computed in agreement with Ma = ell0 - Au
     // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
-    Vector accelerationRHS = ell0;
-    Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
-    solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
-                       parset.sub("a0.solver"));
+    std::vector<Vector> accelerationRHS = ell0;
+    for (size_t i=0; i<bodyCount; i++) {
+      // Initial state
+      alpha[i] = parset.get<double>("boundary.friction.initialAlpha");
+
+      // Initial normal stress
+      assemblers[i]->assembleWeightedNormalStress(
+        frictionalBoundaries[i], weightedNormalStress[i], body.getYoungModulus(),
+        body.getPoissonRatio(), u[i]);
 
-    // Initial state
-    alpha = parset.get<double>("boundary.friction.initialAlpha");
+      Dune::MatrixVector::subtractProduct(accelerationRHS[i], *matrices.elasticity[i], u[i]);
+    }
 
-    // Initial normal stress
-    myAssembler.assembleWeightedNormalStress(
-        frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
-        body.getPoissonRatio(), u);
+    solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
+                       parset.sub("a0.solver"));
   }
 
 public:
-  Vector u;
-  Vector v;
-  Vector a;
-  ScalarVector alpha;
-  ScalarVector weightedNormalStress;
+  std::vector<Vector> u;
+  std::vector<Vector> v;
+  std::vector<Vector> a;
+  std::vector<ScalarVector> alpha;
+  std::vector<ScalarVector> weightedNormalStress;
   double relativeTime;
   double relativeTau;
   size_t timeStep;
+
+private:
+  const size_t bodyCount;
 };
 
 #endif

src/spatial-solving/fixedpointiterator.cc

@@ -8,6 +8,19 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include <dune/contact/common/dualbasisadapter.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/functions/gridfunctions/gridfunction.hh>
+
+#include <dune/geometry/quadraturerules.hh>
+#include <dune/geometry/type.hh>
+#include <dune/geometry/referenceelements.hh>
+
+#include <dune/fufem/functions/basisgridfunction.hh>
+
 #include "../enums.hh"
 #include "../enumparser.hh"
 
@@ -46,37 +59,44 @@ FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
-  ScalarVector alpha;
+  std::vector<ScalarVector> alpha;
   updaters.state_->extractAlpha(alpha);
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
     // solve a velocity problem
-    globalFriction_->updateAlpha(alpha);
+   // globalFriction_->updateAlpha(alpha);
     ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                 velocityRHS, velocityIterate);
     BlockProblem velocityProblem(parset_, convexProblem);
     step_->setProblem(velocityIterate, velocityProblem);
+    //step_->setProblem(velocityIterate);
     velocityProblemSolver.preprocess();
     velocityProblemSolver.solve();
 
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
-    Vector v_m;
+    std::vector<Vector> v_m;
     updaters.rate_->extractOldVelocity(v_m);
-    v_m *= 1.0 - lambda_;
-    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
+
+    for (size_t i=0; i<v_m.size(); i++) {
+      v_m[i] *= 1.0 - lambda_;
+      //Arithmetic::addProduct(v_m[i], lambda_, velocityIterate[i]);
+    }
+
+    // compute relative velocities on contact boundaries
+    relativeVelocities(v_m);
 
     // solve a state problem
     updaters.state_->solve(v_m);
     ScalarVector newAlpha;
-    updaters.state_->extractAlpha(newAlpha);
+   /* updaters.state_->extractAlpha(newAlpha);
 
     if (lambda_ < 1e-12 or
         errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
       fixedPointIteration++;
       break;
     }
-    alpha = newAlpha;
+    alpha = newAlpha;*/
   }
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
@@ -97,4 +117,181 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
+template <class Factory, class Updaters, class ErrorNorm>
+void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
+  // needs assemblers to obtain basis
+  std::vector<std::shared_ptr<MyAssembler>> assemblers(bodyCount);
+
+  using field_type = typename Factory::Matrix::field_type;
+
+
+  // adaptation of DualMortarCoupling::setup()
+
+  const size_t dim = DeformedGrid::dimension;
+  typedef typename DeformedGrid::LeafGridView GridView;
+
+  //cache of local bases
+  typedef Dune::PQkLocalFiniteElementCache<typename DeformedGrid::ctype, field_type, dim,1> FiniteElementCache1;
+  FiniteElementCache1 cache1;
+
+  // cache for the dual functions on the boundary
+  using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
+  std::unique_ptr<DualCache> dualCache;
+  dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
+
+  // define FE grid functions
+  std::vector<BasisGridFunction<VertexBasis, Vector>* > gridFunctions(v_m.size());
+  for (size_t i=0; i<gridFunctions.size(); i++) {
+    gridFunctions[i] = new BasisGridFunction<MyAssembler::VertexBasis, Vector>(assemblers[i]->vertexBasis, v_m[i]);
+  }
+
+  /*
+  for (size_t i=0; i<nBodyAssembler_.nCouplings(); i++) {
+    const auto& coupling = nBodyAssembler_.getCoupling(i);
+    auto glue = coupling.backend();
+
+    const std::array<int, 2> gridIdx = coupling.gridIdx_;
+    const int nonmortarGridIdx = ;
+    const int mortarGridIdx = ;
+
+    // loop over all intersections
+    for (const auto& rIs : intersections(glue)) {
+        const auto& inside = rIs.inside();
+
+        if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+            continue;
+
+        const auto& outside = rIs.outside();
+
+        // types of the elements supporting the boundary segments in question
+        Dune::GeometryType nonmortarEType = inside.type();
+        Dune::GeometryType mortarEType    = outside.type();
+
+        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
+        const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);
+
+        int noOfMortarVec = targetRefElement.size(dim);
+
+        Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
+        Dune::GeometryType mFaceType  = targetRefElement.type(rIs.indexInOutside(),1);
+
+        // Select a quadrature rule
+        // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
+        // are quadrilaterals, then the quad order has to be risen to 3 (4).
+        int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
+        const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);
+
+        const auto& mortarFiniteElement = cache1.get(mortarEType);
+        dualCache->bind(inside, rIs.indexInInside());
+
+        std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;
+
+        const auto& rGeom = rIs.geometry();
+        const auto& rGeomOutside = rIs.geometryOutside();
+        const auto& rGeomInInside = rIs.geometryInInside();
+        const auto& rGeomInOutside = rIs.geometryInOutside();
+
+        int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
+        std::vector<int> nonmortarFaceNodes;
+        for (int i=0; i<nNonmortarFaceNodes; i++) {
+          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
+          nonmortarFaceNodes.push_back(faceIdxi);
+        }
+
+        for (const auto& quadPt : quadRule) {
+
+            // compute integration element of overlap
+            ctype integrationElement = rGeom.integrationElement(quadPt.position());
+
+            // quadrature point positions on the reference element
+            Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
+            Dune::FieldVector<ctype,dim> mortarQuadPos    = rGeomInOutside.global(quadPt.position());
+
+            // The current quadrature point in world coordinates
+            Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
+            Dune::FieldVector<field_type,dim> mortarQpWorld    = rGeomOutside.global(quadPt.position());;
+
+            // the gap direction (normal * gapValue)
+            Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld  - nonmortarQpWorld;
+
+            //evaluate all shapefunctions at the quadrature point
+            //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
+            mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
+            dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);
+
+            // loop over all Lagrange multiplier shape functions
+            for (int j=0; j<nNonmortarFaceNodes; j++) {
+
+                int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
+                int rowIdx = globalToLocal[globalDomainIdx];
+
+                weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
+                    * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);
+
+                // loop over all mortar shape functions
+                for (int k=0; k<noOfMortarVec; k++) {
+
+                    int colIdx  = indexSet1.subIndex(outside, k, dim);
+                    if (!mortarBoundary_.containsVertex(colIdx))
+                        continue;
+
+                    // Integrate over the product of two shape functions
+                    field_type mortarEntry =  integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];
+
+                    Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);
+
+                }
+
+            }
+
+        }
+
+
+
+
+
+
+  }
+
+
+      // Create mapping from the global set of block dofs to the ones on the contact boundary
+      std::vector<int> globalToLocal;
+      nonmortarBoundary_.makeGlobalToLocal(globalToLocal);
+
+      // loop over all intersections
+      for (const auto& rIs : intersections(glue)) {
+
+          if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
+              continue;
+
+          const auto& inside = rIs.inside();
+          const auto& outside = rIs.outside();
+
+          const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+          const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+          int nDomainVertices = domainRefElement.size(dim);
+          int nTargetVertices = targetRefElement.size(dim);
+
+          for (int j=0; j<nDomainVertices; j++) {
+
+              int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];
+
+              // if the vertex is not contained in the restricted contact boundary then dismiss it
+              if (localDomainIdx == -1)
+                  continue;
+
+              for (int k=0; k<nTargetVertices; k++) {
+                  int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
+                  if (!mortarBoundary_.containsVertex(globalTargetIdx))
+                      continue;
+
+                  mortarIndices.add(localDomainIdx, globalTargetIdx);
+              }
+          }
+      }
+    */
+
+}
+
 #include "fixedpointiterator_tmpl.cc"

src/spatial-solving/fixedpointiterator.hh

@@ -8,6 +8,8 @@
 #include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/solvers/solver.hh>
 
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
 struct FixedPointIterationCounter {
   size_t iterations = 0;
   size_t multigridIterations = 0;
@@ -27,10 +29,15 @@ class FixedPointIterator {
   using BlockProblem = typename Factory::BlockProblem;
   using Nonlinearity = typename ConvexProblem::NonlinearityType;
 
+  using DeformedGrid = typename Factory::DeformedGrid;
+
+private:
+  void relativeVelocities(std::vector<Vector>& v_m) const;
+
 public:
-  FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
+  FixedPointIterator(Factory &factory, const Dune::ParameterTree& parset,
                      std::shared_ptr<Nonlinearity> globalFriction,
-                     ErrorNorm const &errorNorm_);
+                     const ErrorNorm& errorNorm);
 
   FixedPointIterationCounter run(Updaters updaters,
                                  Matrix const &velocityMatrix,
@@ -38,6 +45,7 @@ class FixedPointIterator {
                                  Vector &velocityIterate);
 
 private:
+
   std::shared_ptr<typename Factory::Step> step_;
   Dune::ParameterTree const &parset_;
   std::shared_ptr<Nonlinearity> globalFriction_;

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -10,6 +10,8 @@
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
 
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/tectonic/myblockproblem.hh>
 
@@ -22,7 +24,7 @@ using Function = Dune::VirtualFunction<double, double>;
 using Factory = SolverFactory<
     MY_DIM,
     MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
+    DeformedGrid>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
 using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
 using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;

src/spatial-solving/solverfactory.cc

@@ -6,54 +6,73 @@
 #undef HAVE_IPOPT
 #endif
 
+#include <dune/common/bitsetvector.hh>
+
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/solvers/solvers/solver.hh>
 
 #include "solverfactory.hh"
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, Grid const &grid,
+template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
+SolverFactory<DeformedGrid, Matrix, Vector>::SolverFactory(
+    Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
     Dune::BitSetVector<dim> const &ignoreNodes)
-    : baseEnergyNorm(linearBaseSolverStep),
-      linearBaseSolver(&linearBaseSolverStep,
+    : nBodyAssembler_(nBodyAssembler),
+      baseEnergyNorm_(baseSolverStep_),
+      baseSolver_(&baseSolverStep_,
                        parset.get<size_t>("linear.maxiumumIterations"),
-                       parset.get<double>("linear.tolerance"), &baseEnergyNorm,
+                       parset.get<double>("linear.tolerance"), &baseEnergyNorm_,
                        Solver::QUIET),
-      transferOperators(grid.maxLevel()),
-      multigridStep(
-          std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
-  // linear iteration step
-  linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
-                                parset.get<int>("linear.pre"),
-                                parset.get<int>("linear.post"));
-  linearIterationStep.basesolver_ = &linearBaseSolver;
-  linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
-  // transfer operators
-  for (auto &&x : transferOperators)
-    x = new CompressedMultigridTransfer<Vector>;
-  TransferOperatorAssembler<Grid>(grid)
-      .assembleOperatorPointerHierarchy(transferOperators);
-  linearIterationStep.setTransferOperators(transferOperators);
+      transferOperators_(nBodyAssembler.getGrids().at(0)->maxLevel()) {
 
   // tnnmg iteration step
-  multigridStep->setSmoothingSteps(parset.get<int>("main.pre"),
-                                   parset.get<int>("main.multi"),
-                                   parset.get<int>("main.post"));
-  multigridStep->ignoreNodes_ = &ignoreNodes;
+  multigridStep_->setMGType(parset.get<int>("main.multi"), parset.get<int>("main.pre"), parset.get<int>("main.post"));
+  multigridStep_->setIgnore(ignoreNodes);
+  multigridStep_->setBaseSolver(baseSolver_);
+  multigridStep_->setSmoother(&presmoother_, &postsmoother_);
+  multigridStep_->setHasObstacle(nBodyAssembler_.totalHasObstacle_);
+  multigridStep_->setObstacles(nBodyAssembler_.totalObstacles_);
+
+  // create the transfer operators
+  const int nCouplings = nBodyAssembler_.nCouplings();
+  const auto grids = nBodyAssembler_.getGrids();
+  for (size_t i=0; i<transferOperators_.size(); i++) {
+    std::vector<const Dune::BitSetVector<1>*> coarseHasObstacle(nCouplings);
+    std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nCouplings);
+
+    std::vector<const Matrix*> mortarTransfer(nCouplings);
+    std::vector<std::array<int,2> > gridIdx(nCouplings);
+
+    for (int j=0; j<nCouplings; j++) {
+      coarseHasObstacle[j]  = nBodyAssembler_.nonmortarBoundary_[j][i].getVertices();
+      fineHasObstacle[j]    = nBodyAssembler_.nonmortarBoundary_[j][i+1].getVertices();
+
+      mortarTransfer[j] = &nBodyAssembler_.contactCoupling_[j].mortarLagrangeMatrix(i);
+      gridIdx[j]        = nBodyAssembler_.coupling_[j].gridIdx_;
+    }
+
+    transferOperators_[i] = new Dune::Contact::ContactMGTransfer<Vector>;
+    transferOperators_[i]->setup(grids, i, i+1,
+                                    nBodyAssembler_.localCoordSystems_[i],
+                                    nBodyAssembler_.localCoordSystems_[i+1],
+                                    coarseHasObstacle, fineHasObstacle,
+                                    mortarTransfer,
+                                    gridIdx);
+  }
+
+  multigridStep_->setTransferOperators(transferOperators_);
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
-  for (auto &&x : transferOperators)
+template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
+SolverFactory<DeformedGrid, Matrix, Vector>::~SolverFactory() {
+  for (auto &&x : transferOperators_)
     delete x;
 }
 
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
+template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
+auto SolverFactory<DeformedGrid, Matrix, Vector>::getStep()
     -> std::shared_ptr<Step> {
-  return multigridStep;
+  return multigridStep_;
 }
 
 #include "solverfactory_tmpl.cc"

src/spatial-solving/solverfactory.hh

@@ -5,29 +5,39 @@
 #include <dune/common/parametertree.hh>
 
 #include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+//#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
+//#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
+//#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
+#include <dune/contact/assemblers/nbodyassembler.hh>
+
+#include <dune/contact/estimators/hierarchiccontactestimator.hh>
+#include <dune/contact/solvers/nsnewtonmgstep.hh>
+#include <dune/contact/solvers/contacttransfer.hh>
+#include <dune/contact/solvers/contactobsrestrict.hh>
+#include <dune/contact/solvers/contacttransferoperatorassembler.hh>
+#include <dune/contact/solvers/nsnewtoncontacttransfer.hh>
+
+#define USE_OLD_TNNMG //needed for old bisection.hh in tnnmg
+
+template <class DeformedGridTEMPLATE, class MatrixType, class VectorType>
 class SolverFactory {
+protected:
+  const size_t dim = DeformedGrid::dimension;
+
 public:
-  using BlockProblem = BlockProblemTEMPLATE;
-  using ConvexProblem = typename BlockProblem::ConvexProblemType;
-  using Vector = typename BlockProblem::VectorType;
-  using Matrix = typename BlockProblem::MatrixType;
+  using Vector = VectorType;
+  using Matrix = MatrixType;
 
-private:
-  using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
+  using DeformedGrid = DeformedGridTEMPLATE;
 
 public:
-  using Step =
-      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
+  using Step = Dune::Contact::NonSmoothNewtonMGStep<Matrix, Vector>;
 
-  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
+  SolverFactory(Dune::ParameterTree const &parset, const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler,
                 Dune::BitSetVector<dim> const &ignoreNodes);
 
   ~SolverFactory();
@@ -35,14 +45,16 @@ class SolverFactory {
   std::shared_ptr<Step> getStep();
 
 private:
-  TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
-  EnergyNorm<Matrix, Vector> baseEnergyNorm;
-  LoopSolver<Vector> linearBaseSolver;
-  TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
-  TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
-  MultigridStep<Matrix, Vector> linearIterationStep;
-  std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
-  NonlinearSmoother nonlinearSmoother;
-  std::shared_ptr<Step> multigridStep;
+  const Dune::Contact::NBodyAssembler<DeformedGrid, Vector>& nBodyAssembler_;
+
+  ProjectedBlockGSStep<Matrix, Vector> baseSolverStep_;
+  EnergyNorm<Matrix, Vector> baseEnergyNorm_;
+  LoopSolver<Vector> baseSolver_;
+
+  ProjectedBlockGSStep<Matrix, Vector> presmoother_;
+  ProjectedBlockGSStep<Matrix, Vector> postsmoother_;
+  std::shared_ptr<Step> multigridStep_;
+
+  std::vector<Dune::Contact::ContactMGTransfer<Vector>* > transferOperators_;
 };
 #endif

src/time-stepping/rate.cc

@@ -9,20 +9,22 @@
 template <class Vector, class Matrix, class Function, int dimension>
 std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
 initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial) {
+                const std::vector<Function>& velocityDirichletFunctions,
+                const std::vector<Dune::BitSetVector<dimension>>& velocityDirichletNodes,
+                const Matrices<Matrix>& matrices,
+                const std::vector<Vector>& u_initial,
+                const std::vector<Vector>& v_initial,
+                const std::vector<Vector>& a_initial) {
   switch (scheme) {
     case Config::Newmark:
       return std::make_shared<Newmark<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     case Config::BackwardEuler:
       return std::make_shared<
           BackwardEuler<Vector, Matrix, Function, dimension>>(
           matrices, u_initial, v_initial, a_initial, velocityDirichletNodes,
-          velocityDirichletFunction);
+          velocityDirichletFunctions);
     default:
       assert(false);
   }