Connectivity.hpp

#ifndef CONNECTIVITY_HPP
#define CONNECTIVITY_HPP

#include <PastisAssert.hpp>
#include <TinyVector.hpp>

#include <Kokkos_Core.hpp>

#include <ConnectivityMatrix.hpp>
#include <ConnectivityComputer.hpp>
#include <SubItemValuePerItem.hpp>

#include <vector>
#include <unordered_map>
#include <algorithm>

#include <RefId.hpp>
#include <TypeOfItem.hpp>
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>

#include <tuple>
#include <algorithm>

#include <IConnectivity.hpp>

template <size_t Dimension>
class Connectivity;

template <size_t Dimension>
class ConnectivityFace;

template<>
class ConnectivityFace<1>
{
 public:
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const;
  };
};

template<>
class ConnectivityFace<2>
{
 public:
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const {
      size_t hash = 0;
      hash ^= std::hash<unsigned int>()(f.m_node0_id);
      hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
      return hash;
    }
  };

  unsigned int m_node0_id;
  unsigned int m_node1_id;

  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
  {
    os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
    return os;
  }

  KOKKOS_INLINE_FUNCTION
  bool operator==(const ConnectivityFace& f) const
  {
    return ((m_node0_id == f.m_node0_id) and
            (m_node1_id == f.m_node1_id));
  }

  KOKKOS_INLINE_FUNCTION
  bool operator<(const ConnectivityFace& f) const
  {
    return ((m_node0_id<f.m_node0_id) or
            ((m_node0_id == f.m_node0_id) and
             (m_node1_id<f.m_node1_id)));
  }

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace& operator=(const ConnectivityFace&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace& operator=(ConnectivityFace&&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
  {
    Assert(given_node_id_list.size()==2);
#warning rework this dirty constructor
    const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
    m_node0_id = min;
    m_node1_id = max;
  }

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(const ConnectivityFace&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(ConnectivityFace&&) = default;

  KOKKOS_INLINE_FUNCTION
  ~ConnectivityFace() = default;
};

template <>
class ConnectivityFace<3>
{
 private:
  friend class Connectivity<3>;
  friend struct Hash;
  struct Hash
  {
    size_t operator()(const ConnectivityFace& f) const {
      size_t hash = 0;
      for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
        hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
      }
      return hash;
    }
  };

  bool m_reversed;
  std::vector<unsigned int> m_node_id_list;

  friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
  {
    for (auto id : f.m_node_id_list) {
      std::cout << id << ' ';
    }
    return os;
  }

  KOKKOS_INLINE_FUNCTION
  const bool& reversed() const
  {
    return m_reversed;
  }

  KOKKOS_INLINE_FUNCTION
  const std::vector<unsigned int>& nodeIdList() const
  {
    return m_node_id_list;
  }

  KOKKOS_INLINE_FUNCTION
  std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
  {
    const auto min_id = std::min_element(node_list.begin(), node_list.end());
    const int shift = std::distance(node_list.begin(), min_id);

    std::vector<unsigned int> rotated_node_list(node_list.size());
    if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
      for (size_t i=0; i<node_list.size(); ++i) {
        rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
        m_reversed = true;
      }
    } else {
      for (size_t i=0; i<node_list.size(); ++i) {
        rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
      }
    }

    return rotated_node_list;
  }

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
      : m_reversed(false),
        m_node_id_list(_sort(given_node_id_list))
  {
    ;
  }

 public:
  bool operator==(const ConnectivityFace& f) const
  {
    if (m_node_id_list.size() == f.nodeIdList().size()) {
      for (size_t j=0; j<m_node_id_list.size(); ++j) {
        if (m_node_id_list[j] != f.nodeIdList()[j]) {
          return false;
        }
      }
      return true;
    }
    return false;
  }

  KOKKOS_INLINE_FUNCTION
  bool operator<(const ConnectivityFace& f) const
  {
    const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
    for (size_t i=0; i<min_nb_nodes; ++i) {
      if (m_node_id_list[i] <  f.m_node_id_list[i]) return true;
      if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
    }
    return m_node_id_list.size() < f.m_node_id_list.size();
  }

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace& operator=(const ConnectivityFace&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace& operator=(ConnectivityFace&&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(const ConnectivityFace&) = default;

  KOKKOS_INLINE_FUNCTION
  ConnectivityFace(ConnectivityFace&&) = default;


  KOKKOS_INLINE_FUNCTION
  ConnectivityFace() = delete;

  KOKKOS_INLINE_FUNCTION
  ~ConnectivityFace() = default;
};

template <size_t Dimension>
class Connectivity final
    : public IConnectivity
{
public:
  static constexpr size_t dimension = Dimension;
  static constexpr bool has_edges = (dimension>2);
  static constexpr bool has_face = (dimension>1);

  ConnectivityMatrix m_cell_to_node_matrix;
  NodeValuePerCell<unsigned short> m_cell_to_node_local_cell;

  ConnectivityMatrix m_cell_to_face_matrix;
  FaceValuePerCell<bool> m_cell_face_is_reversed;

  ConnectivityMatrix m_face_to_cell_matrix;
  CellValuePerFace<unsigned short> m_face_to_cell_local_face;
  ConnectivityMatrix m_face_to_node_matrix;

  ConnectivityMatrix m_node_to_cell_matrix;
  CellValuePerNode<unsigned short> m_node_to_cell_local_node;

  template <TypeOfItem SubItemType,
            TypeOfItem ItemType>
  const ConnectivityMatrix& itemToItemMatrix() const = delete;

  const ConnectivityMatrix& itemToItemMatrix(const TypeOfItem& item_type_0,
                                             const TypeOfItem& item_type_1) const final;

private:
  ConnectivityComputer m_connectivity_computer;

  std::vector<RefFaceList> m_ref_face_list;
  std::vector<RefNodeList> m_ref_node_list;

  Kokkos::View<double*> m_inv_cell_nb_nodes;

  Kokkos::View<const unsigned short*> m_node_nb_faces;
  Kokkos::View<const unsigned int**> m_node_faces;

  using Face = ConnectivityFace<Dimension>;

  std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map;

  void _computeFaceCellConnectivities();

 public:
  void addRefFaceList(const RefFaceList& ref_face_list)
  {
    m_ref_face_list.push_back(ref_face_list);
  }

  size_t numberOfRefFaceList() const
  {
    return m_ref_face_list.size();
  }

  const RefFaceList& refFaceList(const size_t& i) const
  {
    return m_ref_face_list[i];
  }

  void addRefNodeList(const RefNodeList& ref_node_list)
  {
    m_ref_node_list.push_back(ref_node_list);
  }

  size_t numberOfRefNodeList() const
  {
    return m_ref_node_list.size();
  }

  const RefNodeList& refNodeList(const size_t& i) const
  {
    return m_ref_node_list[i];
  }

  KOKKOS_INLINE_FUNCTION
  size_t numberOfNodes() const
  {
    return m_node_to_cell_matrix.numRows();
  }

  KOKKOS_INLINE_FUNCTION
  size_t numberOfFaces() const
  {
    return m_face_to_cell_matrix.numRows();
  }

  KOKKOS_INLINE_FUNCTION
  size_t numberOfCells() const
  {
    return m_cell_to_node_matrix.numRows();
  }

  const Kokkos::View<const double*> invCellNbNodes() const
  {
    return m_inv_cell_nb_nodes;
  }

  unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const
  {
    const Face face(face_nodes);
    auto i_face = m_face_number_map.find(face);
    if (i_face == m_face_number_map.end()) {
      std::cerr << "Face " << face << " not found!\n";
      throw std::exception();
      std::exit(0);
    }
    return i_face->second;
  }

  Connectivity(const Connectivity&) = delete;

  Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
  {
    m_cell_to_node_matrix = cell_by_node_vector;

    Assert(this->numberOfCells()>0);

    {
      Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
      Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
          const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
          inv_cell_nb_nodes[j] = 1./cell_nodes.length;
        });
      m_inv_cell_nb_nodes = inv_cell_nb_nodes;
    }

    m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_node_matrix,
                                                             m_node_to_cell_matrix);

    m_node_to_cell_local_node = CellValuePerNode<unsigned short>(*this);

    m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_node_matrix,
                                                              m_node_to_cell_matrix,
                                                              m_node_to_cell_local_node);

    m_cell_to_node_local_cell = NodeValuePerCell<unsigned short>(*this);

    m_connectivity_computer.computeLocalChildItemNumberInItem(m_node_to_cell_matrix,
                                                              m_cell_to_node_matrix,
                                                              m_cell_to_node_local_cell);
    if constexpr (Dimension>1) {
      this->_computeFaceCellConnectivities();
    }
  }

  ~Connectivity()
  {
    ;
  }
};


using Connectivity3D = Connectivity<3>;

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<3>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::face>() const
{
  return m_cell_to_face_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<3>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::node>() const
{
  return m_cell_to_node_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<3>::itemToItemMatrix<TypeOfItem::face,
                                  TypeOfItem::cell>() const
{
  return m_face_to_cell_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<3>::itemToItemMatrix<TypeOfItem::node,
                                  TypeOfItem::cell>() const
{
  return m_node_to_cell_matrix;
}


using Connectivity2D = Connectivity<2>;

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<2>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::face>() const
{
  return m_cell_to_face_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<2>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::node>() const
{
  return m_cell_to_node_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<2>::itemToItemMatrix<TypeOfItem::face,
                                  TypeOfItem::cell>() const
{
  return m_face_to_cell_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<2>::itemToItemMatrix<TypeOfItem::node,
                                  TypeOfItem::cell>() const
{
  return m_node_to_cell_matrix;
}

using Connectivity1D = Connectivity<1>;

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<1>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::node>() const
{
  return m_cell_to_node_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<1>::itemToItemMatrix<TypeOfItem::cell,
                                  TypeOfItem::face>() const
{
#warning in 1d, faces and node are the same
  return m_cell_to_face_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<1>::itemToItemMatrix<TypeOfItem::face,
                                  TypeOfItem::cell>() const
{
#warning in 1d, faces and node are the same
  return m_face_to_cell_matrix;
}

template <>
template <>
inline const ConnectivityMatrix&
Connectivity<1>::itemToItemMatrix<TypeOfItem::node,
                                  TypeOfItem::cell>() const
{
#warning in 1d, faces and node are the same
  return m_node_to_cell_matrix;
}

template <size_t Dimension>
const ConnectivityMatrix&
Connectivity<Dimension>::
itemToItemMatrix(const TypeOfItem& item_type_0,
                 const TypeOfItem& item_type_1) const
{
  switch (item_type_0) {
    case TypeOfItem::cell: {
      switch (item_type_1) {
        case TypeOfItem::node: {
          return itemToItemMatrix<TypeOfItem::cell, TypeOfItem::node>();
        }
        case TypeOfItem::face: {
          return itemToItemMatrix<TypeOfItem::cell, TypeOfItem::face>();
        }
        default: {
          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
          std::exit(1);
        }
      }
    }
    case TypeOfItem::face: {
      switch (item_type_1) {
        case TypeOfItem::cell: {
          return itemToItemMatrix<TypeOfItem::face, TypeOfItem::cell>();
        }
        default: {
          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
          std::exit(1);
        }
      }
    }
    case TypeOfItem::node: {
      switch (item_type_1) {
        case TypeOfItem::cell: {
          return itemToItemMatrix<TypeOfItem::node, TypeOfItem::cell>();
        }
        default: {
          std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n";
          std::exit(1);
        }
      }
    }
    default: {
      std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_0) << "\n";
      std::exit(1);
    }
  }
}


#endif // CONNECTIVITY_HPP