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

#include <dune/fufem/geometry/polyhedrondistance.hh>

#include "mygrids.hh"
#include "../midpoint.hh"
#include "../../utils/diameter.hh"

#include "simplexmanager.hh"

// Fix: 3D case (still Elias code)
template <class Grid> GridsConstructor<Grid>::GridsConstructor(const std::vector<std::shared_ptr<CuboidGeometry<typename Grid::ctype>>>& cuboidGeometries_) :
  cuboidGeometries(cuboidGeometries_)
{
  size_t const gridCount = cuboidGeometries.size();
  grids.resize(gridCount);
  gridFactories.resize(gridCount);

  for (size_t idx=0; idx<grids.size(); idx++) {
    const auto& cuboidGeometry = *cuboidGeometries[idx];

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

    grids[idx] = std::shared_ptr<Grid>(gridFactories[idx].createGrid());
    grids[idx]->setRefinementType(Grid::RefinementType::COPY);
  }
}

template <class Grid>
std::vector<std::shared_ptr<Grid>>& GridsConstructor<Grid>::getGrids() {
  return grids;
}

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

template <class GridView>
template <class Vector>
bool MyFaces<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 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 * 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 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, 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
  }

double computeAdmissibleDiameter(double distance, double smallestDiameter, double lengthScale) {
  return (distance / 0.0125 / lengthScale + 1.0) * smallestDiameter;
}

template <class Grid, class LocalVector>
void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
            double smallestDiameter, double lengthScale) {
  bool needRefine = true;
  while (true) {
    needRefine = false;
    for (auto &&e : elements(grid.leafGridView())) {
      auto const geometry = e.geometry();

      auto const weakeningRegionDistance =
          distance(weakPatch, geometry, 1e-6 * lengthScale);
      auto const admissibleDiameter =
          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter, lengthScale);

      if (diameter(geometry) <= admissibleDiameter)
        continue;

      needRefine = true;
      grid.mark(1, e);
    }
    if (!needRefine)
      break;

    grid.preAdapt();
    grid.adapt();
    grid.postAdapt();
  }
}

#include "mygrids_tmpl.cc"