#ifndef MESH_NODE_BOUNDARY_HPP
#define MESH_NODE_BOUNDARY_HPP

#include <Kokkos_Vector.hpp>
#include <algebra/TinyVector.hpp>
#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/IConnectivity.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/RefItemList.hpp>
#include <utils/Array.hpp>
#include <utils/Exceptions.hpp>
#include <utils/Messenger.hpp>

template <size_t Dimension>
class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
{
 protected:
  Array<const NodeId> m_node_list;
  std::string m_boundary_name;

 public:
  MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
  MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;

  const Array<const NodeId>&
  nodeList() const
  {
    return m_node_list;
  }

  template <typename MeshType>
  MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
    : m_boundary_name(ref_face_list.refId().tagName())
  {
    static_assert(Dimension == MeshType::Dimension);
    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();

    const Array<const FaceId>& face_list = ref_face_list.list();
    parallel_for(
      face_list.size(), PUGS_LAMBDA(int l) {
        const auto& face_cells = face_to_cell_matrix[face_list[l]];
        if (face_cells.size() > 1) {
          std::ostringstream ost;
          ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
              << ": inner faces cannot be used to define mesh boundaries";
          throw NormalError(ost.str());
        }
      });

    Kokkos::vector<unsigned int> node_ids;
    // not enough but should reduce significantly the number of resizing
    node_ids.reserve(Dimension * face_list.size());
    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();

    for (size_t l = 0; l < face_list.size(); ++l) {
      const FaceId face_number = face_list[l];
      const auto& face_nodes   = face_to_node_matrix[face_number];

      for (size_t r = 0; r < face_nodes.size(); ++r) {
        node_ids.push_back(face_nodes[r]);
      }
    }
    std::sort(node_ids.begin(), node_ids.end());
    auto last = std::unique(node_ids.begin(), node_ids.end());
    node_ids.resize(std::distance(node_ids.begin(), last));

    Array<NodeId> node_list(node_ids.size());
    parallel_for(
      node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
    m_node_list = node_list;
  }

  template <typename MeshType>
  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())

  {
    static_assert(Dimension == MeshType::Dimension);
  }

  MeshNodeBoundary()          = default;
  virtual ~MeshNodeBoundary() = default;

 protected:
  MeshNodeBoundary(const MeshNodeBoundary&) = default;
  MeshNodeBoundary(MeshNodeBoundary&&)      = default;
};

template <size_t Dimension>
class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
{
 public:
  using Rd = TinyVector<Dimension, double>;

 private:
  const Rd m_outgoing_normal;

  template <typename MeshType>
  PUGS_INLINE Rd _getNormal(const MeshType& mesh);

  template <typename MeshType>
  PUGS_INLINE void
  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
                       const TinyVector<Dimension, double>& origin,
                       const double length,
                       const MeshType& mesh) const
  {
    static_assert(MeshType::Dimension == Dimension);

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

    parallel_for(
      this->m_node_list.size(), PUGS_LAMBDA(const size_t r) {
        const Rd& x = xr[this->m_node_list[r]];
        if (dot(x - origin, normal) > 1E-13 * length) {
          std::ostringstream ost;
          ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
              << ": boundary is not flat!";
          throw NormalError(ost.str());
        }
      });
  }

  template <typename MeshType>
  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);

 public:
  const Rd&
  outgoingNormal() const
  {
    return m_outgoing_normal;
  }

  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;

  template <typename MeshType>
  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
  {
    ;
  }

  template <typename MeshType>
  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
  {
    ;
  }

  MeshFlatNodeBoundary()                            = default;
  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
  virtual ~MeshFlatNodeBoundary()                   = default;
};

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<1, double>
MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 1);
  using R = TinyVector<1, double>;

  const size_t number_of_bc_nodes = [&]() {
    size_t number_of_bc_nodes = 0;
    auto node_is_owned        = mesh.connectivity().nodeIsOwned();
    for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
      number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
    }
    return parallel::allReduceMax(number_of_bc_nodes);
  }();

  if (number_of_bc_nodes != 1) {
    std::ostringstream ost;
    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
        << ": node boundaries in 1D require to have exactly 1 node";
    throw NormalError(ost.str());
  }

  return R{1};
}

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<2, double>
MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 2);
  using R2 = TinyVector<2, double>;

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

  R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
  R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());

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

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

  for (size_t r = 0; r < m_node_list.size(); ++r) {
    const R2& x = xr[m_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);
    }
  }

  if (xmin == xmax) {
    std::ostringstream ost;
    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
        << ": unable to compute normal";
    throw NormalError(ost.str());
  }

  R2 dx = xmax - xmin;
  dx *= 1. / l2Norm(dx);

  const R2 normal(-dx[1], dx[0]);

  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);

  return normal;
}

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<3, double>
MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 3);
  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[1]) 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[2] < zmin[2])) or
        ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
      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[2] > zmax[2])) or
        ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
      zmax = x;
    }
  };

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

  R3 xmax = -xmin;
  R3 ymax = xmax;
  R3 zmax = xmax;

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

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

  if (parallel::size() > 0) {
    Array<const R3> xmin_array = parallel::allGather(xmin);
    Array<const R3> xmax_array = parallel::allGather(xmax);
    Array<const R3> ymin_array = parallel::allGather(ymin);
    Array<const R3> ymax_array = parallel::allGather(ymax);
    Array<const R3> zmin_array = parallel::allGather(zmin);
    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);
    }
  }

  const R3 u = xmax - xmin;
  const R3 v = ymax - ymin;
  const R3 w = zmax - zmin;

  const R3 uv        = crossProduct(u, v);
  const double uv_l2 = dot(uv, uv);

  R3 normal        = uv;
  double normal_l2 = uv_l2;

  const R3 uw        = crossProduct(u, w);
  const double uw_l2 = dot(uw, uw);

  if (uw_l2 > uv_l2) {
    normal    = uw;
    normal_l2 = uw_l2;
  }

  const R3 vw        = crossProduct(v, w);
  const double vw_l2 = dot(vw, vw);

  if (vw_l2 > normal_l2) {
    normal    = vw;
    normal_l2 = vw_l2;
  }

  if (normal_l2 == 0) {
    std::ostringstream ost;
    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
        << ": unable to compute normal";
    throw NormalError(ost.str());
  }

  const double length = sqrt(normal_l2);

  normal *= 1. / length;

  this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);

  return normal;
}

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<1, double>
MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 1);
  using R = TinyVector<1, double>;

  const R normal = this->_getNormal(mesh);

  double max_height = 0;

  if (m_node_list.size() > 0) {
    const NodeValue<const R>& xr    = mesh.xr();
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();

    const NodeId r0      = m_node_list[0];
    const CellId j0      = node_to_cell_matrix[r0][0];
    const auto& j0_nodes = cell_to_node_matrix[j0];

    for (size_t r = 0; r < j0_nodes.size(); ++r) {
      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
      if (std::abs(height) > std::abs(max_height)) {
        max_height = height;
      }
    }
  }

  Array<double> max_height_array = parallel::allGather(max_height);
  for (size_t i = 0; i < max_height_array.size(); ++i) {
    const double height = max_height_array[i];
    if (std::abs(height) > std::abs(max_height)) {
      max_height = height;
    }
  }

  if (max_height > 0) {
    return -normal;
  } else {
    return normal;
  }
}

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<2, double>
MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 2);
  using R2 = TinyVector<2, double>;

  const R2 normal = this->_getNormal(mesh);

  double max_height = 0;

  if (m_node_list.size() > 0) {
    const NodeValue<const R2>& xr   = mesh.xr();
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();

    const NodeId r0      = m_node_list[0];
    const CellId j0      = node_to_cell_matrix[r0][0];
    const auto& j0_nodes = cell_to_node_matrix[j0];
    for (size_t r = 0; r < j0_nodes.size(); ++r) {
      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
      if (std::abs(height) > std::abs(max_height)) {
        max_height = height;
      }
    }
  }

  Array<double> max_height_array = parallel::allGather(max_height);
  for (size_t i = 0; i < max_height_array.size(); ++i) {
    const double height = max_height_array[i];
    if (std::abs(height) > std::abs(max_height)) {
      max_height = height;
    }
  }

  if (max_height > 0) {
    return -normal;
  } else {
    return normal;
  }
}

template <>
template <typename MeshType>
PUGS_INLINE TinyVector<3, double>
MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
{
  static_assert(MeshType::Dimension == 3);
  using R3 = TinyVector<3, double>;

  const R3 normal = this->_getNormal(mesh);

  double max_height = 0;

  if (m_node_list.size() > 0) {
    const NodeValue<const R3>& xr   = mesh.xr();
    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();

    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();

    const NodeId r0      = m_node_list[0];
    const CellId j0      = node_to_cell_matrix[r0][0];
    const auto& j0_nodes = cell_to_node_matrix[j0];

    for (size_t r = 0; r < j0_nodes.size(); ++r) {
      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
      if (std::abs(height) > std::abs(max_height)) {
        max_height = height;
      }
    }
  }

  Array<double> max_height_array = parallel::allGather(max_height);
  for (size_t i = 0; i < max_height_array.size(); ++i) {
    const double height = max_height_array[i];
    if (std::abs(height) > std::abs(max_height)) {
      max_height = height;
    }
  }

  if (max_height > 0) {
    return -normal;
  } else {
    return normal;
  }
}

#endif   // MESH_NODE_BOUNDARY_HPP