#include <mesh/MeshNodeBoundaryUtils.hpp>

#include <mesh/Connectivity.hpp>
#include <mesh/Mesh.hpp>
#include <utils/Messenger.hpp>

template <>
std::array<TinyVector<2>, 2>
getBounds(const Mesh<2>& mesh, const RefNodeList& ref_node_list)
{
  using R2 = TinyVector<2, double>;

  const NodeValue<const R2>& xr = mesh.xr();

  std::array<R2, 2> bounds;
  R2& xmin = bounds[0];
  R2& xmax = bounds[1];

  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
  xmax = R2{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest()};

  auto update_xmin = [](const R2& x, R2& x_min) {
    if ((x[0] < x_min[0]) or ((x[0] == x_min[0]) and (x[1] < x_min[1]))) {
      x_min = x;
    }
  };

  auto update_xmax = [](const R2& x, R2& x_max) {
    if ((x[0] > x_max[0]) or ((x[0] == x_max[0]) and (x[1] > x_max[1]))) {
      x_max = x;
    }
  };

  auto node_list = ref_node_list.list();
  for (size_t r = 0; r < node_list.size(); ++r) {
    const R2& x = xr[node_list[r]];
    update_xmin(x, xmin);
    update_xmax(x, xmax);
  }

  if (parallel::size() > 1) {
    Array<R2> xmin_array = parallel::allGather(xmin);
    Array<R2> xmax_array = parallel::allGather(xmax);
    for (size_t i = 0; i < xmin_array.size(); ++i) {
      update_xmin(xmin_array[i], xmin);
    }
    for (size_t i = 0; i < xmax_array.size(); ++i) {
      update_xmax(xmax_array[i], xmax);
    }
  }

  return bounds;
}

template <>
std::array<TinyVector<3>, 6>
getBounds(const Mesh<3>& mesh, const RefNodeList& ref_node_list)
{
  using R3 = TinyVector<3, double>;

  auto update_xmin = [](const R3& x, R3& xmin) {
    // XMIN: X.xmin X.ymax X.zmax
    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
      xmin = x;
    }
  };

  auto update_xmax = [](const R3& x, R3& xmax) {
    // XMAX: X.xmax X.ymin X.zmin
    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
      xmax = x;
    }
  };

  auto update_ymin = [](const R3& x, R3& ymin) {
    // YMIN: X.ymin X.zmax X.xmin
    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
      ymin = x;
    }
  };

  auto update_ymax = [](const R3& x, R3& ymax) {
    // YMAX: X.ymax X.zmin X.xmax
    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
        ((x[1] == ymax[1]) and (x[2] == ymax[2]) and (x[0] > ymax[0]))) {
      ymax = x;
    }
  };

  auto update_zmin = [](const R3& x, R3& zmin) {
    // ZMIN: X.zmin X.xmin X.ymin
    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[0] < zmin[0])) or
        ((x[2] == zmin[2]) and (x[0] == zmin[0]) and (x[1] < zmin[1]))) {
      zmin = x;
    }
  };

  auto update_zmax = [](const R3& x, R3& zmax) {
    // ZMAX: X.zmax X.xmax X.ymax
    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[0] > zmax[0])) or
        ((x[2] == zmax[2]) and (x[0] == zmax[0]) and (x[1] > zmax[1]))) {
      zmax = x;
    }
  };

  std::array<R3, 6> bounds;
  R3& xmin = bounds[0];
  R3& ymin = bounds[1];
  R3& zmin = bounds[2];
  R3& xmax = bounds[3];
  R3& ymax = bounds[4];
  R3& zmax = bounds[5];

  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
  ymin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
  zmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};

  xmax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
            std::numeric_limits<double>::lowest()};
  ymax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
            std::numeric_limits<double>::lowest()};
  zmax = R3{std::numeric_limits<double>::lowest(), std::numeric_limits<double>::lowest(),
            std::numeric_limits<double>::lowest()};

  const NodeValue<const R3>& xr = mesh.xr();

  auto node_list = ref_node_list.list();
  for (size_t r = 0; r < node_list.size(); ++r) {
    const R3& x = xr[node_list[r]];
    update_xmin(x, xmin);
    update_ymin(x, ymin);
    update_zmin(x, zmin);
    update_xmax(x, xmax);
    update_ymax(x, ymax);
    update_zmax(x, zmax);
  }

  if (parallel::size() > 1) {
    Array<const R3> xmin_array = parallel::allGather(xmin);
    Array<const R3> ymin_array = parallel::allGather(ymin);
    Array<const R3> zmin_array = parallel::allGather(zmin);
    Array<const R3> xmax_array = parallel::allGather(xmax);
    Array<const R3> ymax_array = parallel::allGather(ymax);
    Array<const R3> zmax_array = parallel::allGather(zmax);

    for (size_t i = 0; i < xmin_array.size(); ++i) {
      update_xmin(xmin_array[i], xmin);
    }
    for (size_t i = 0; i < ymin_array.size(); ++i) {
      update_ymin(ymin_array[i], ymin);
    }
    for (size_t i = 0; i < zmin_array.size(); ++i) {
      update_zmin(zmin_array[i], zmin);
    }
    for (size_t i = 0; i < xmax_array.size(); ++i) {
      update_xmax(xmax_array[i], xmax);
    }
    for (size_t i = 0; i < ymax_array.size(); ++i) {
      update_ymax(ymax_array[i], ymax);
    }
    for (size_t i = 0; i < zmax_array.size(); ++i) {
      update_zmax(zmax_array[i], zmax);
    }
  }

  return bounds;
}