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, 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 * 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);
};

template <class Grid> class GridsConstructor {
public:
  GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_);

  std::vector<std::shared_ptr<Grid>>& getGrids();

  template <class GridView>
  MyFaces<GridView> constructFaces(const GridView& gridView, const CuboidGeometry<typename Grid::ctype>& cuboidGeometry);

private:
  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 lengthScale);

template <class Grid, class LocalVector>
void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
            double smallestDiameter, double lengthScale);

#endif