dune/tectonic/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"

dune/tectonic/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
+
+

dune/tectonic/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>;

dune/tectonic/problem-data/grid/cubegridconstructor.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBEGRIDCONSTRUCTOR_HH
+
+#include "gridconstructor.hh"
+
+#include <dune/common/fmatrix.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+
+
+template <class GridType>
+class CubeGridConstructor : public GridConstructor {    
+public:
+    using Cube = Cube<GridType::dimensionworld>;
+    
+    CubeGridConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries_);
+
+    template <class GridView>
+    void constructFaces(const GridView& gridView, CuboidGeometry const &cuboidGeometry, CubeFaces<GridView>& cubeFaces);
+
+private:
+  std::vector<std::shared_ptr<CuboidGeometry>> const &cuboidGeometries;
+};
+
+
+#endif
+
+

dune/tectonic/problem-data/grid/cuboidgeometry.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <fstream>
+
+#ifdef HAVE_CAIROMM
+#include <cairomm/context.h>
+#include <cairomm/fontface.h>
+#include <cairomm/surface.h>
+#endif
+
+#include "cuboidgeometry.hh"
+
+#if MY_DIM == 3
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+                                      const double length = 1.00, const double height = 0.27, const double depth = 0.60) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    depth_(depth*lengthScale()),
+    lowerLeft_(origin),
+    lowerRight_({origin[0]+length_, origin[1], 0}),
+    upperRight_({origin[0]+length_, origin[1]+height_, 0}),
+    upperLeft_({origin[0], origin[1]+height_, 0})
+    {}
+#else
+
+template <class ctype>
+CuboidGeometry<ctype>::CuboidGeometry(const GlobalCoords& origin,
+                                      const double length, const double height) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    lowerLeft_(origin),
+    lowerRight_({origin[0]+length_, origin[1]}),
+    upperRight_({origin[0]+length_, origin[1]+height_}),
+    upperLeft_({origin[0], origin[1]+height_})
+    {}
+#endif
+
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningRegion(const WeakeningRegion& weakeningRegion) {
+    weakeningRegions_.emplace_back(weakeningRegion);
+}
+
+template <class ctype>
+void CuboidGeometry<ctype>::addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1) {
+    WeakeningRegion weakPatch;
+
+#if MY_DIM == 3 // TODO: Does not work yet
+    if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(4);
+        weakPatch.vertices[0] = weakPatch.vertices[2] = vertex0;
+        weakPatch.vertices[1] = weakPatch.vertices[3] = vertex1;
+
+        for (size_t k = 0; k < 2; ++k) {
+            weakPatch.vertices[k][2] = -depth_ / 2.0;
+            weakPatch.vertices[k + 2][2] = depth_ / 2.0;
+        }
+        switch (parset.get<Config::PatchType>("patchType")) {
+            case Config::Rectangular:
+                break;
+            case Config::Trapezoidal:
+                weakPatch.vertices[1][0] += 0.05 * lengthScale();
+                weakPatch.vertices[3][0] -= 0.05 * lengthScale();
+                break;
+            default:
+                assert(false);
+        }
+        addWeakeningRegion(weakPatch);
+    }
+#else
+    if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(2);
+        weakPatch.vertices[0] = vertex0;
+        weakPatch.vertices[1] = vertex1;
+        addWeakeningRegion(weakPatch);
+    }
+#endif
+}
+
+template <class ctype>
+void CuboidGeometry<ctype>::write() const {
+    std::fstream writer("geometry", std::fstream::out);
+    writer << "lowerLeft = " << lowerLeft_ << std::endl;
+    writer << "lowerRight = " << lowerRight_ << std::endl;
+    writer << "upperRight = " << upperRight_ << std::endl;
+    writer << "upperLeft = " << upperLeft_ << std::endl;
+}
+
+/*
+template <class ctype>
+void CuboidGeometry<ctype>::render() const {
+#ifdef HAVE_CAIROMM
+  std::string const filename = "geometry.png";
+  double const width = 600;
+  double const height = 400;
+  double const widthScale = 400;
+  double const heightScale = 400;
+
+  auto surface =
+      Cairo::ImageSurface::create(Cairo::FORMAT_ARGB32, width, height);
+  auto cr = Cairo::Context::create(surface);
+
+  auto const setRGBColor = [&](int colour) {
+    cr->set_source_rgb(((colour & 0xFF0000) >> 16) / 255.0,
+                       ((colour & 0x00FF00) >> 8) / 255.0,
+                       ((colour & 0x0000FF) >> 0) / 255.0);
+  };
+  auto const moveTo = [&](LocalVector2D const &v) { cr->move_to(v[0], -v[1]); };
+  auto const lineTo = [&](LocalVector2D const &v) { cr->line_to(v[0], -v[1]); };
+
+  cr->scale(widthScale, heightScale);
+  cr->translate(0.1, 0.1);
+  cr->set_line_width(0.0025);
+
+  // triangle
+  {
+    moveTo(reference::A);
+    lineTo(reference::B);
+    lineTo(reference::C);
+    cr->close_path();
+    cr->stroke();
+  }
+
+  // dashed lines
+  {
+    cr->save();
+    std::vector<double> dashPattern = { 0.005 };
+    cr->set_dash(dashPattern, 0);
+    moveTo(reference::Z);
+    lineTo(reference::Y);
+    moveTo(reference::U);
+    lineTo(reference::X);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // fill viscoelastic region
+  {
+    cr->save();
+    setRGBColor(0x0097E0);
+    moveTo(reference::B);
+    lineTo(reference::K);
+    lineTo(reference::M);
+    cr->fill();
+    cr->restore();
+  }
+
+  // mark weakening region
+  {
+    cr->save();
+    setRGBColor(0x7AD3FF);
+    cr->set_line_width(0.005);
+    moveTo(reference::X);
+    lineTo(reference::Y);
+    cr->stroke();
+    cr->restore();
+  }
+
+  // mark points
+  {
+    auto const drawCircle = [&](LocalVector2D const &v) {
+      cr->arc(v[0], -v[1], 0.0075, -M_PI, M_PI); // x,y,radius,angle1,angle2
+      cr->fill();
+    };
+
+    cr->save();
+    setRGBColor(0x002F47);
+    drawCircle(reference::A);
+    drawCircle(reference::B);
+    drawCircle(reference::C);
+    drawCircle(reference::Y);
+    drawCircle(reference::X);
+    drawCircle(reference::Z);
+    drawCircle(reference::U);
+    drawCircle(reference::K);
+    drawCircle(reference::M);
+    drawCircle(reference::G);
+    drawCircle(reference::H);
+    drawCircle(reference::J);
+    drawCircle(reference::I);
+    cr->restore();
+  }
+
+  // labels
+  {
+    auto const label = [&](LocalVector2D const &v, std::string l) {
+      moveTo(v);
+      cr->rel_move_to(0.005, -0.02);
+      cr->show_text(l);
+    };
+    auto font = Cairo::ToyFontFace::create(
+        "monospace", Cairo::FONT_SLANT_NORMAL, Cairo::FONT_WEIGHT_NORMAL);
+
+    cr->save();
+    cr->set_font_face(font);
+    cr->set_font_size(0.03);
+
+    label(A, "A");
+    label(B, "B");
+    label(C, "C");
+    label(D, "D");
+    label(X, "X");
+    label(Y, "Y");
+
+    cr->restore();
+  }
+
+  surface->write_to_png(filename);
+#endif
+}
+*/
+
+#include "cuboidgeometry_tmpl.cc"

dune/tectonic/problem-data/grid/cuboidgeometry.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+template <class ctype = double>
+class CuboidGeometry {
+public:
+    typedef Dune::FieldVector<ctype, MY_DIM> GlobalCoords;
+    using WeakeningRegion = ConvexPolyhedron<GlobalCoords>;
+
+    static constexpr double lengthScale() {
+        return 1.0;
+    } // scaling factor
+
+private:
+    const ctype length_;
+    const ctype height_;
+#if MY_DIM == 3
+    const ctype depth_;
+#endif
+
+    // corners of cube
+    const GlobalCoords lowerLeft_;
+    const GlobalCoords lowerRight_;
+    const GlobalCoords upperRight_;
+    const GlobalCoords upperLeft_;
+
+    // weakening regions
+    std::vector<WeakeningRegion> weakeningRegions_;
+
+public:
+#if MY_DIM == 3
+    CuboidGeometry(const GlobalCoords& origin,
+                   const double length = 1.00, const double height = 0.27, const double depth = 0.60);
+
+    const auto& depth() const {
+        return depth_;
+    }
+#else
+    CuboidGeometry(const GlobalCoords& origin,
+                       const double length = 1.00, const double height = 0.27);
+#endif
+
+    void addWeakeningRegion(const WeakeningRegion& weakeningRegion);
+    void addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1);
+
+    const auto& lowerLeft() const {
+        return lowerLeft_;
+    }
+
+    const auto& lowerRight() const {
+        return lowerRight_;
+    }
+
+    const auto& upperRight() const {
+        return upperRight_;
+    }
+
+    const auto& upperLeft() const {
+        return upperLeft_;
+    }
+
+    const auto& weakeningRegions() const {
+        return weakeningRegions_;
+    }
+
+    void write() const;
+
+   // void render() const;
+};
+#endif

dune/tectonic/problem-data/grid/cuboidgeometry_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "cuboidgeometry.hh"
+
+template class CuboidGeometry<typename Grid::ctype>;

dune/tectonic/problem-data/grid/gridconstructor.hh

+#ifndef SRC_GRIDCONSTRUCTOR_HH
+#define SRC_GRIDCONSTRUCTOR_HH
+
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "../../utils/diameter.hh"
+
+template <class GridType>
+class GridConstructor {
+    private:
+        double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
+            return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
+        }
+
+    public:
+        GridConstructor() {}
+
+        std::shared_ptr<GridType> grid() {
+            return grid_;
+        }
+
+        template <class LocalVector>
+        void refine(const ConvexPolyhedron<LocalVector>& 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();
+            }
+        }
+
+        virtual void createGrid() = 0;
+
+    private:
+        Dune::GridFactory<GridType> gridFactory_;
+        std::shared_ptr<GridType> grid_;
+};
+
+
+
+
+#endif
+
+

dune/tectonic/problem-data/grid/mygrids.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "mygrids.hh"
+#include "../midpoint.hh"
+#include "../../utils/diameter.hh"
+
+#include "simplexmanager.hh"
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& 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.lowerLeft();
+    const auto& B = cuboidGeometry.lowerRight();
+    const auto& C = cuboidGeometry.upperRight();
+    const auto& D = cuboidGeometry.upperLeft();
+
+    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());
+    grids[idx]->setRefinementType(Grid::RefinementType::COPY);
+  }
+}
+
+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, const CuboidGeometry<typename Grid::ctype>& 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, const CuboidGeometry<typename GridView::ctype>& 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.lowerLeft(), cuboidGeometry.lowerRight(), in.geometry().center());
+    });
+
+    right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.lowerRight(), cuboidGeometry.upperRight(), in.geometry().center());
+    });
+
+    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.upperLeft(), cuboidGeometry.upperRight(), in.geometry().center());
+    });
+
+    left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(cuboidGeometry.lowerLeft(), cuboidGeometry.upperLeft(), 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/one-body-problem-data/mygrid.hh → dune/tectonic/problem-data/grid/mygrids.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_MYGRID_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>
@@ -7,29 +7,32 @@
 #include <dune/fufem/boundarypatch.hh>
 #include <dune/fufem/geometry/convexpolyhedron.hh>
 
-#include "mygeometry.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);
+  MyFaces(GridView const &gridView, const CuboidGeometry<typename GridView::ctype>& cuboidGeometry_);
 
 private:
+  const CuboidGeometry<typename GridView::ctype>& cuboidGeometry;
+
   bool isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14 * MyGeometry::lengthScale;
-  };
+    return std::abs(a - b) < 1e-14 * cuboidGeometry.lengthScale();
+  }
 
   bool isClose2(double a, double b) {
     return std::abs(a - b) <
-           1e-14 * MyGeometry::lengthScale * MyGeometry::lengthScale;
-  };
+           1e-14 * cuboidGeometry.lengthScale() * cuboidGeometry.lengthScale();
+  }
 
   template <class Vector>
   bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x);
@@ -41,44 +44,27 @@ template <class GridView> struct MyFaces {
   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 GridConstructor {
+template <class Grid> class GridsConstructor {
 public:
-  GridConstructor();
+  GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_);
 
-  std::shared_ptr<Grid> getGrid();
+  std::vector<std::shared_ptr<Grid>>& getGrids();
 
   template <class GridView>
-  MyFaces<GridView> constructFaces(GridView const &gridView);
+  MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry);
 
 private:
-  Dune::GridFactory<Grid> gridFactory;
+  const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries;
+  std::vector<Dune::GridFactory<Grid>> gridFactories;
+  std::vector<std::shared_ptr<Grid>> grids;
 };
 
-double computeAdmissibleDiameter(double distance, double smallestDiameter);
+double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale);
 
 template <class Grid, class LocalVector>
 void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
-            double smallestDiameter);
+            double smallestDiameter, double lengthScale);
 
 #endif
+
+

dune/tectonic/problem-data/grid/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<DeformedGrid::LeafGridView>;
+template struct MyFaces<DeformedGrid::LevelGridView>;
+
+template MyFaces<DeformedGrid::LeafGridView> GridsConstructor<Grid>::constructFaces(
+    DeformedGrid::LeafGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
+
+template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
+    DeformedGrid::LevelGridView const &gridView, CuboidGeometry<typename Grid::ctype> const &CuboidGeometry_);
+
+template void refine<Grid, LocalVector>(
+    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+    double smallestDiameter, double lengthScale);

dune/tectonic/problem-data/grid/simplexmanager.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include "simplexmanager.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_;
+}

dune/tectonic/problem-data/grid/simplexmanager.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_SIMPLEXMANAGER_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_SIMPLEXMANAGER_HH
+
+#include <vector>
+
+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);
+
+  auto getSimplices() -> SimplexList const &;
+
+private:
+  SimplexList simplices_;
+
+#if MY_DIM == 3
+  unsigned int const shift_;
+#endif
+};
+
+#endif
+
+

src/one-body-problem-data/mygeometry.cc → dune/tectonic/problem-data/grid/trianglegeometry.cc

@@ -10,27 +10,84 @@
 #include <cairomm/surface.h>
 #endif
 
-#include "mygeometry.hh"
-
-void MyGeometry::write() {
-  std::fstream writer("geometry", std::fstream::out);
-  writer << "A = " << A << std::endl;
-  writer << "B = " << B << std::endl;
-  writer << "C = " << C << std::endl;
-  writer << "Y = " << Y << std::endl;
-  writer << "X = " << X << std::endl;
-  writer << "Z = " << Z << std::endl;
-  writer << "U = " << U << std::endl;
-  writer << "K = " << K << std::endl;
-  writer << "M = " << M << std::endl;
-  writer << "G = " << G << std::endl;
-  writer << "H = " << H << std::endl;
-  writer << "J = " << J << std::endl;
-  writer << "I = " << I << std::endl;
-  writer << "zenith = " << zenith << std::endl;
+#include "trianglegeometry.hh"
+
+#if MY_DIM == 3
+template <class ctype>
+TriangleGeometry<ctype>::TriangleGeometry(const GlobalCoords& origin,
+                                      const double length, const double height, const double depth) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    depth_(depth*lengthScale()),
+    A_(origin),
+    B_({origin[0]+length_, origin[1], 0}),
+    C_({origin[0]+length_, origin[1]+height_, 0})
+    {}
+#else
+
+template <class ctype>
+TriangleGeometry<ctype>::TriangleGeometry(const GlobalCoords& origin,
+                                      const double length, const double height) :
+    length_(length*lengthScale()),
+    height_(height*lengthScale()),
+    A_(origin),
+    B_({origin[0]+length_, origin[1]}),
+    C_({origin[0]+length_, origin[1]+height_})
+    {}
+#endif
+
+template <class ctype>
+void TriangleGeometry<ctype>::addWeakeningRegion(const WeakeningRegion& weakeningRegion) {
+    weakeningRegions_.emplace_back(weakeningRegion);
+}
+
+template <class ctype>
+void TriangleGeometry<ctype>::addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1) {
+    WeakeningRegion weakPatch;
+
+#if MY_DIM == 3 // TODO: Does not work yet
+    /*if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(4);
+        weakPatch.vertices[0] = weakPatch.vertices[2] = vertex0;
+        weakPatch.vertices[1] = weakPatch.vertices[3] = vertex1;
+
+        for (size_t k = 0; k < 2; ++k) {
+            weakPatch.vertices[k][2] = -depth_ / 2.0;
+            weakPatch.vertices[k + 2][2] = depth_ / 2.0;
+        }
+        switch (parset.get<Config::PatchType>("patchType")) {
+            case Config::Rectangular:
+                break;
+            case Config::Trapezoidal:
+                weakPatch.vertices[1][0] += 0.05 * lengthScale();
+                weakPatch.vertices[3][0] -= 0.05 * lengthScale();
+                break;
+            default:
+                assert(false);
+        }
+        addWeakeningRegion(weakPatch);
+    } */
+#else
+    if (vertex0 != vertex1) {
+        weakPatch.vertices.resize(2);
+        weakPatch.vertices[0] = vertex0;
+        weakPatch.vertices[1] = vertex1;
+        addWeakeningRegion(weakPatch);
+    }
+#endif
 }
 
-void MyGeometry::render() {
+template <class ctype>
+void TriangleGeometry<ctype>::write() const {
+    std::fstream writer("geometry", std::fstream::out);
+    writer << "A = " << A_ << std::endl;
+    writer << "B = " << B_ << std::endl;
+    writer << "C = " << C_ << std::endl;
+}
+
+/*
+template <class ctype>
+void CuboidGeometry<ctype>::render() const {
 #ifdef HAVE_CAIROMM
   std::string const filename = "geometry.png";
   double const width = 600;
@@ -137,22 +194,19 @@ void MyGeometry::render() {
     cr->set_font_face(font);
     cr->set_font_size(0.03);
 
-    label(reference::A, "A");
-    label(reference::B, "B");
-    label(reference::C, "C");
-    label(reference::K, "K");
-    label(reference::M, "M");
-    label(reference::U, "U");
-    label(reference::X, "X");
-    label(reference::Y, "Y");
-    label(reference::Z, "Z");
-    label(reference::G, "G");
-    label(reference::H, "H");
-    label(reference::J, "J");
-    label(reference::I, "I");
+    label(A, "A");
+    label(B, "B");
+    label(C, "C");
+    label(D, "D");
+    label(X, "X");
+    label(Y, "Y");
+
     cr->restore();
   }
 
   surface->write_to_png(filename);
 #endif
 }
+*/
+
+#include "trianglegeometry_tmpl.cc"

dune/tectonic/problem-data/grid/trianglegeometry.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGEOMETRY_HH
+
+#include <vector>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+//            C
+//          / |
+//         /  |
+//        /   |height
+//       /    |
+//      A --- B
+//       length
+
+template <class ctype = double>
+class TriangleGeometry {
+public:
+    typedef Dune::FieldVector<ctype, MY_DIM> GlobalCoords;
+    using WeakeningRegion = ConvexPolyhedron<GlobalCoords>;
+
+    static constexpr double lengthScale() {
+        return 1.0;
+    } // scaling factor
+
+private:
+    const ctype length_;
+    const ctype height_;
+#if MY_DIM == 3
+    const ctype depth_;
+#endif
+
+    // corners of triangle
+    const GlobalCoords A_;
+    const GlobalCoords B_;
+    const GlobalCoords C_;
+
+    // weakening regions
+    std::vector<WeakeningRegion> weakeningRegions_;
+
+public:
+#if MY_DIM == 3
+    TriangleGeometry(const GlobalCoords& origin,
+                   const double length = 0.5, const double height = 0.5, const double depth = 0.12);
+
+    const auto& depth() const {
+        return depth_;
+    }
+#else
+    TriangleGeometry(const GlobalCoords& origin,
+                       const double length = 0.5, const double height = 0.5);
+#endif
+
+    void addWeakeningRegion(const WeakeningRegion& weakeningRegion);
+    void addWeakeningPatch(const Dune::ParameterTree& parset, const GlobalCoords& vertex0, const GlobalCoords& vertex1);
+
+    const auto& A() const {
+        return A_;
+    }
+
+    const auto& B() const {
+        return B_;
+    }
+
+    const auto& C() const {
+        return C_;
+    }
+
+    const auto& weakeningRegions() const {
+        return weakeningRegions_;
+    }
+
+    const auto& length() const {
+        return length_;
+    }
+
+    const auto& height() const {
+        return height_;
+    }
+
+    void write() const;
+
+   // void render() const;
+};
+#endif

dune/tectonic/problem-data/grid/trianglegeometry_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../../explicitgrid.hh"
+#include "trianglegeometry.hh"
+
+template class TriangleGeometry<typename Grid::ctype>;

dune/tectonic/problem-data/grid/trianglegrids.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
+#include "trianglegrids.hh"
+#include "../midpoint.hh"
+#include "../../utils/diameter.hh"
+
+#include "simplexmanager.hh"
+
+// Fix: 3D case (still Elias code)
+template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries) :
+  triangleGeometries_(triangleGeometries)
+{
+  size_t const gridCount = triangleGeometries.size();
+  grids.resize(gridCount);
+  gridFactories.resize(gridCount);
+
+  for (size_t idx=0; idx<grids.size(); idx++) {
+    const auto& triangleGeometry = *triangleGeometries[idx];
+
+    const auto& A = triangleGeometry.A();
+    const auto& B = triangleGeometry.B();
+    const auto& C = triangleGeometry.C();
+
+    const size_t vc = 3;
+
+#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];
+        assert(k == j * vc);
+      }
+    }
+
+#if MY_DIM == 3
+    for (size_t k = 0; k < vc; ++k) {
+      vertices[k][2] = 0;
+      vertices[k + vc][2] = triangleGeometry.depth();
+    }
+#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);
+    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());
+    grids[idx]->setRefinementType(Grid::RefinementType::COPY);
+  }
+}
+
+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, const TriangleGeometry<typename Grid::ctype>& triangleGeometry) {
+  return MyFaces<GridView>(gridView, triangleGeometry);
+}
+
+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 * triangleGeometry_.lengthScale() > x[0] or
+      x[0] > minmax0.second + 1e-14 * triangleGeometry_.lengthScale())
+    return false;
+  if (minmax1.first - 1e-14 * triangleGeometry_.lengthScale() > x[1] or
+      x[1] > minmax1.second + 1e-14 * triangleGeometry_.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, const TriangleGeometry<typename GridView::ctype>& triangleGeometry)
+      :
+  #if MY_DIM == 3
+        a(gridView),
+        b(gridView),
+        c(gridView),
+        top(gridView),
+        bottom(gridView),
+  #else
+        a(gridView),
+        b(gridView),
+        c(gridView),
+  #endif
+        triangleGeometry_(triangleGeometry)
+  {
+    a.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(triangleGeometry_.B(), triangleGeometry_.C(), in.geometry().center());
+    });
+
+    b.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(triangleGeometry_.A(), triangleGeometry_.C(), in.geometry().center());
+    });
+
+    c.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return xyBetween(triangleGeometry_.A(), triangleGeometry_.B(), in.geometry().center());
+    });
+
+  #if MY_DIM == 3
+    top.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(triangleGeometry_.depth(), in.geometry().center()[2]);
+    });
+
+    bottom.insertFacesByProperty([&](typename GridView::Intersection const &in) {
+      return isClose(0.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 "trianglegrids_tmpl.cc"

dune/tectonic/problem-data/grid/trianglegrids.hh

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGRIDS_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_TRIANGLEGRIDS_HH
+
+#include <dune/common/fmatrix.hh>
+#include <dune/grid/common/gridfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+
+#include "trianglegeometry.hh"
+
+//            C
+//          / |
+//         /  |
+//      b /   | a
+//       /    |
+//      A --- B
+//         c
+
+template <class GridView> struct MyFaces {
+  BoundaryPatch<GridView> a;
+  BoundaryPatch<GridView> b;
+  BoundaryPatch<GridView> c;
+
+#if MY_DIM == 3
+  BoundaryPatch<GridView> top;
+  BoundaryPatch<GridView> bottom;
+#endif
+
+  MyFaces(GridView const &gridView, const TriangleGeometry<typename GridView::ctype>& triangleGeometry);
+
+private:
+  const TriangleGeometry<typename GridView::ctype>& triangleGeometry_;
+
+  bool isClose(double a, double b) {
+    return std::abs(a - b) < 1e-14 * triangleGeometry_.lengthScale();
+  }
+
+  bool isClose2(double a, double b) {
+    return std::abs(a - b) <
+           1e-14 * triangleGeometry_.lengthScale() * triangleGeometry_.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);
+};
+
+template <class Grid> class GridsConstructor {
+public:
+  GridsConstructor(const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries);
+
+  std::vector<std::shared_ptr<Grid>>& getGrids();
+
+  template <class GridView>
+  MyFaces<GridView> constructFaces(const GridView& gridView, const TriangleGeometry<typename Grid::ctype>& triangleGeometry);
+
+private:
+  const std::vector<std::shared_ptr<TriangleGeometry<typename Grid::ctype>>>& triangleGeometries_;
+  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
+
+

dune/tectonic/problem-data/grid/trianglegrids_tmpl.cc

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

src/one-body-problem-data/midpoint.hh → dune/tectonic/problem-data/midpoint.hh

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