cubefaces.cc

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

#include "cubefaces.hh"

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

  return true;
}

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

template <class GridView>
CubeFaces<GridView>::CubeFaces(const GridView& gridView, const Cube& cube, double lengthScale)
      :
  #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
        cube_(cube),
        lengthScale_(lengthScale)
  {
    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
      return xyBetween(cuboidGeometry.A, cuboidGeometry.B, in.geometry().center());
    });

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

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

    left.insertFacesByProperty([&](typename GridView::Intersection const &in) {
      return xyBetween(cuboidGeometry.A, cuboidGeometry.D, 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
  }

#include "cubefaces_tmpl.cc"