cuboidgridconstructor.hh

#ifndef SRC_CUBOIDGRIDCONSTRUCTOR_HH
#define SRC_CUBOIDGRIDCONSTRUCTOR_HH

#include <dune/common/fmatrix.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/grid/utility/structuredgridfactory.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/geometry/convexpolyhedron.hh>

#include "cuboidgeometry.hh"
#include "gridconstructor.hh"
#include "simplexmanager.hh"

// -------------------------------------------------
// MyFaces
// -------------------------------------------------

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_) :
    #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
  }

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) {
      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 Vector>
  bool 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 Vector>
  bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x) {
      return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
  }
};


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

// -------------------------------------------------
// CuboidGridConstructor
// -------------------------------------------------

template <class Grid>
class CuboidGridConstructor : public GridConstructor<Grid> {
public:
  CuboidGridConstructor(const CuboidGeometry<typename Grid::ctype>& cuboidGeometry) :
    cuboidGeometry_(cuboidGeometry) {

    createGrid();
  }

protected:
  void createGrid() {
    Dune::GridFactory<Grid> gridFactory;

    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)
                  gridFactory.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);
                  gridFactory.insertElement(cell, simplices[i]);
                }

                this->grid_ = std::shared_ptr<Grid>(gridFactory.createGrid());
                this->grid_->setRefinementType(Grid::RefinementType::COPY);
              }

    protected:
        const CuboidGeometry<typename Grid::ctype>& cuboidGeometry_;
};


// -------------------------------------------------
// DefaultCuboidGridConstructor
// -------------------------------------------------
template <class Grid>
class DefaultCuboidGridConstructor : public GridConstructor<Grid> {
  static const int dim = Grid::dimension;

public:
  DefaultCuboidGridConstructor(const CuboidGeometry<typename Grid::ctype>& cuboidGeometry, std::array<unsigned int, dim> elements) :
    cuboidGeometry_(cuboidGeometry),
    elements_(elements) {

    createGrid();
  }

protected:
  void createGrid() {
       this->grid_ = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
           cuboidGeometry_.lowerLeft(), cuboidGeometry_.upperRight(), elements_);
    }

    protected:
        const CuboidGeometry<typename Grid::ctype>& cuboidGeometry_;
        std::array<unsigned int, dim> elements_;
};

#endif