mygrid.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "mygrid.hh"
#include "midpoint.hh"

class SimplexManager {
public:
  using SimplexList = std::vector<std::vector<unsigned int>>;

#if DIM == 3
  SimplexManager(unsigned int shift) : shift_(shift) {}
#endif

  void addFromVertices(unsigned int U, unsigned int V, unsigned int W) {
#if DIM == 3
    unsigned int const U2 = U + shift_;
    unsigned int const V2 = V + shift_;
    unsigned int const W2 = W + shift_;

    // back-to-front, back-to-front, back-to-front
    simplices_.push_back({ U, V, W, U2 });
    simplices_.push_back({ V, V2, W2, U2 });
    simplices_.push_back({ W, W2, U2, V });

// back-to-front, front-to-back, back-to-front would be
/*
  simplices_.push_back({ U, V, W, U2 });
  simplices_.push_back({ V, V2, W, U2 });
  simplices_.push_back({ V2, W2, U2, W });
*/
#else
    simplices_.push_back({ U, V, W });
#endif
  }

  SimplexList const &getSimplices() { return simplices_; }

private:
  SimplexList simplices_;

#if DIM == 3
  unsigned int const shift_;
#endif
};

template <class Grid> class GridConstructor {
public:
  GridConstructor() {
    auto const &A = MyGeometry::A;
    auto const &B = MyGeometry::B;
    auto const &C = MyGeometry::C;

    auto const AB = midPoint(A, B);
    auto const AC = midPoint(A, C);
    auto const BC = midPoint(B, C);

    auto const AAB = midPoint(A, AB);
    auto const AAC = midPoint(A, AC);
    auto const ABB = midPoint(AB, B);
    auto const ABC = midPoint(AB, C); // the interiour point
    auto const ACC = midPoint(AC, C);

    unsigned int const vc = 11;

#if DIM == 3
    Dune::FieldMatrix<double, 2 * vc, DIM> vertices;
#else
    Dune::FieldMatrix<double, vc, DIM> vertices;
#endif
    for (size_t i = 0; i < 2; ++i) {
      size_t k = 0;
      vertices[k++][i] = A[i];
      vertices[k++][i] = AAB[i];
      vertices[k++][i] = AB[i];
      vertices[k++][i] = ABB[i];
      vertices[k++][i] = B[i];
      vertices[k++][i] = AAC[i];
      vertices[k++][i] = AC[i];
      vertices[k++][i] = ABC[i];
      vertices[k++][i] = BC[i];
      vertices[k++][i] = ACC[i];
      vertices[k++][i] = C[i];
      assert(k == vc);
#if DIM == 3
      vertices[k++][i] = A[i];
      vertices[k++][i] = AAB[i];
      vertices[k++][i] = AB[i];
      vertices[k++][i] = ABB[i];
      vertices[k++][i] = B[i];
      vertices[k++][i] = AAC[i];
      vertices[k++][i] = AC[i];
      vertices[k++][i] = ABC[i];
      vertices[k++][i] = BC[i];
      vertices[k++][i] = ACC[i];
      vertices[k++][i] = C[i];
      assert(k == 2 * vc);
#endif
    }

#if DIM == 3
    for (size_t k = 0; k < vc; ++k) {
      vertices[k][2] = -MyGeometry::depth / 2.0;
      vertices[k + vc][2] = MyGeometry::depth / 2.0;
    }
#endif

    for (size_t i = 0; i < vertices.N(); ++i)
      gridFactory.insertVertex(vertices[i]);

    Dune::GeometryType cell;
#if DIM == 3
    cell.makeTetrahedron();
#else
    cell.makeTriangle();
#endif

#if DIM == 3
    SimplexManager sm(vc);
#else
    SimplexManager sm;
#endif
    sm.addFromVertices(1, 5, 0);
    sm.addFromVertices(1, 6, 5);
    sm.addFromVertices(2, 6, 1);
    sm.addFromVertices(7, 6, 2);
    sm.addFromVertices(7, 2, 3);
    sm.addFromVertices(7, 3, 8);
    sm.addFromVertices(7, 8, 10);
    sm.addFromVertices(7, 10, 9);
    sm.addFromVertices(7, 9, 6);
    sm.addFromVertices(8, 3, 4);
    auto const &simplices = sm.getSimplices();

    // sanity-check choices of simplices
    for (size_t i = 0; i < simplices.size(); ++i) {
      Dune::FieldMatrix<double, DIM, DIM> check;
      for (size_t j = 0; j < DIM; ++j)
        check[j] = vertices[simplices[i][j + 1]] - vertices[simplices[i][j]];
      assert(check.determinant() > 0);
      gridFactory.insertElement(cell, simplices[i]);
    }
  }

  std::shared_ptr<Grid> getGrid() {
    return std::shared_ptr<Grid>(gridFactory.createGrid());
  }

  template <class GridView>
  MyFaces<GridView> constructFaces(GridView const &gridView) {
    return MyFaces<GridView>(gridView);
  }

private:
  Dune::GridFactory<Grid> gridFactory;
};

template <class GridView>
template <class Vector>
bool MyFaces<GridView>::xyCollinear(Vector const &a, Vector const &b,
                                    Vector const &c) {
  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
}

template <class GridView>
template <class Vector>
bool MyFaces<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 > x[0] or x[0] > minmax0.second + 1e-14)
    return false;
  if (minmax1.first - 1e-14 > x[1] or x[1] > minmax1.second + 1e-14)
    return false;

  return true;
}

template <class GridView>
template <class Vector>
bool MyFaces<GridView>::xyBetween(Vector const &v1, Vector const &v2,
                                  Vector const &x) {
  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
}

template <class GridView>
MyFaces<GridView>::MyFaces(GridView const &gridView)
    :
#if DIM == 3
      lower(gridView),
      right(gridView),
      upper(gridView),
      front(gridView),
      back(gridView)
#else
      lower(gridView),
      right(gridView),
      upper(gridView)
#endif
{
  using Grid = typename GridView::Grid;
  assert(isClose(MyGeometry::A[1], 0));
  assert(isClose(MyGeometry::B[1], 0));
  lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
    return isClose(0, in.geometry().center()[1]);
  });
#if DIM == 3
  front.insertFacesByProperty([&](typename GridView::Intersection const &in) {
    return isClose(-MyGeometry::depth / 2.0, in.geometry().center()[2]);
  });

  back.insertFacesByProperty([&](typename GridView::Intersection const &in) {
    return isClose(MyGeometry::depth / 2.0, in.geometry().center()[2]);
  });
#endif

  upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
    return xyBetween(MyGeometry::A, MyGeometry::C, in.geometry().center());
  });

  right.insertFacesByProperty([&](typename GridView::Intersection const &in) {
    return xyBetween(MyGeometry::B, MyGeometry::C, in.geometry().center());
  });
}

#include "mygrid_tmpl.cc"