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

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

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

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

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

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

src/multi-body-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
+
+

src/multi-body-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"
+
+/*
+const CuboidGeometry::LocalVector CuboidGeometry::lift(const double v0, const double v1) {
+  vc = 0;
+  vc[0] = vc2D[0];
+  vc[1] = vc2D[1];
+}
+
+const CuboidGeometry::LocalVector& CuboidGeometry::A() const { return A_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::B() const { return B_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::C() const { return C_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::D() const { return D_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::X() const { return X_; }
+const CuboidGeometry::LocalVector& CuboidGeometry::Y() const { return Y_; }
+double CuboidGeometry::depth() const { return depth_; }
+
+void CuboidGeometry::setupWeak(const LocalVector& weakOrigin, const double weakLength) {
+    lift(weakOrigin, X_);
+    const LocalVector2D Y({X_[0]+weakLength, X_[1]});
+    lift(Y, Y_);
+}
+*/
+
+#if MY_DIM == 3
+CuboidGeometry::CuboidGeometry(const LocalVector& origin,
+                               const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
+                               const double lowerWeakLength, const double upperWeakLength,
+                               const double length const double width, const double depth):
+    length_(length*lengthScale),
+    width_(width*lengthScale),
+    A(origin),
+    B({origin[0]+length_, origin[1], 0}),
+    C({origin[0]+length_, origin[1]+width_, 0}),
+    D({origin[0], origin[1]+width_, 0}),
+    V(lowerWeakOrigin),
+    W({V[0]+lowerWeakLength*lengthScale, V[1], 0}),
+    X(upperWeakOrigin),
+    Y({X[0]+upperWeakLength*lengthScale, X[1], 0}),
+    depth(depth_*lengthScale) {}
+#else
+CuboidGeometry::CuboidGeometry(const LocalVector& origin,
+                               const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
+                               const double lowerWeakLength, const double upperWeakLength,
+                               const double length, const double width):
+    length_(length*lengthScale),
+    width_(width*lengthScale),
+    A(origin),
+    B({origin[0]+length_, origin[1]}),
+    C({origin[0]+length_, origin[1]+width_}),
+    D({origin[0], origin[1]+width_}),
+    V(lowerWeakOrigin),
+    W({V[0]+lowerWeakLength*lengthScale, V[1]}),
+    X(upperWeakOrigin),
+    Y({X[0]+upperWeakLength*lengthScale, X[1]}) {}
+#endif
+
+
+void CuboidGeometry::write() const {
+  std::fstream writer("geometry", std::fstream::out);
+  writer << "A = " << A << std::endl;
+  writer << "B = " << B << std::endl;
+  writer << "C = " << C << std::endl;
+  writer << "D = " << D << std::endl;
+  writer << "V = " << V << std::endl;
+  writer << "W = " << W << std::endl;
+  writer << "X = " << X << std::endl;
+  writer << "Y = " << Y << std::endl;
+}
+
+/*
+void CuboidGeometry::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
+}
+*/

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

+#ifndef SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+#define SRC_MULTI_BODY_PROBLEM_DATA_CUBOIDGEOMETRY_HH
+
+#include <dune/common/fvector.hh>
+
+#include "../midpoint.hh"
+
+class CuboidGeometry {
+  public:
+    typedef Dune::FieldVector<double, MY_DIM> LocalVector;
+
+    constexpr static double const lengthScale = 1.0; // scaling factor
+
+  private:
+    const double length_;
+    const double width_;
+   // const double lowerWeakLength_;     // length of the lower weak zone
+   // const double upperWeakLength_;     // length of the lower weak zone
+
+  public:
+    // corners of cube
+    const LocalVector A;
+    const LocalVector B;
+    const LocalVector C;
+    const LocalVector D;
+
+    // corners of lower weakening patch
+    const LocalVector V;
+    const LocalVector W;
+
+    // corners of upper weakening patch
+    const LocalVector X;
+    const LocalVector Y;
+
+#if MY_DIM == 3
+    const double depth;
+
+    CuboidGeometry(const LocalVector& origin,
+                   const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
+                   const double lowerWeakLength = 0.20, const double upperWeakLength = 0.20,
+                   const double length = 1.00, const double width = 0.27, const double depth = 0.60);
+#else
+        CuboidGeometry(const LocalVector& origin,
+                       const LocalVector& lowerWeakOrigin, const LocalVector& upperWeakOrigin,
+                       const double lowerWeakLength = 0.20, const double upperWeakLength = 0.20,
+                       const double length = 1.00, const double width = 0.27);
+#endif
+
+    void write() const;
+
+   // void render() const;
+};
+#endif

src/multi-body-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 const &CuboidGeometry_);
+
+template MyFaces<DeformedGrid::LevelGridView> GridsConstructor<Grid>::constructFaces(
+    DeformedGrid::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/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
+
+