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