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

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

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

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

// back-to-front, front-to-back, front-to-back
void SimplexManager::addFromVerticesFBB(unsigned int U, unsigned int V,
                                        unsigned int W) {
#if MY_DIM == 3
  unsigned int const U2 = U + shift_;
  unsigned int const V2 = V + shift_;
  unsigned int const W2 = W + shift_;

  simplices_.push_back({ U, V, W, U2 });
  simplices_.push_back({ V, V2, W2, U2 });
  simplices_.push_back({ W, W2, U2, V });
#else
  simplices_.push_back({ U, V, W });
#endif
}

// back-to-front, back-to-front, front-to-back
void SimplexManager::addFromVerticesFFB(unsigned int U, unsigned int V,
                                        unsigned int W) {
#if MY_DIM == 3
  unsigned int const U2 = U + shift_;
  unsigned int const V2 = V + shift_;
  unsigned int const W2 = W + shift_;

  simplices_.push_back({ U, V, W, U2 });
  simplices_.push_back({ V, V2, W, U2 });
  simplices_.push_back({ V2, W, U2, W2 });
#else
  simplices_.push_back({ U, V, W });
#endif
}

auto SimplexManager::getSimplices() -> SimplexList const & {
  return simplices_;
}

// Fix: 3D case (still Elias code)
template <class Grid> GridsConstructor<Grid>::GridsConstructor(std::vector<std::shared_ptr<CuboidGeometry>> const &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.A;
    const auto& B = cuboidGeometry.B;
    const auto& C = cuboidGeometry.C;
    const auto& D = cuboidGeometry.D;

    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());
  }
}

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, CuboidGeometry const &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, CuboidGeometry const &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.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
  }

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"