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

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/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