Connectivity2D.hpp

#ifndef CONNECTIVITY_2D_HPP
#define CONNECTIVITY_2D_HPP

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

#include <ConnectivityUtils.hpp>
#include <vector>
#include <map>
#include <algorithm>

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

class Connectivity2D
{
 public:
  static constexpr size_t dimension = 2;

  ConnectivityMatrix m_cell_to_node_matrix;

  ConnectivityMatrix m_face_to_cell_matrix;
  ConnectivityMatrix m_face_to_node_matrix;

  ConnectivityMatrix m_node_to_cell_matrix;
  ConnectivityMatrixShort m_node_to_cell_local_node_matrix;

  // Stores numbering of nodes of each cell.
  //
  // This is different from m_cell_to_node_matrix which return the global id of
  // a local node in a cell
  ConnectivityMatrix m_node_id_per_cell_matrix;

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

  Kokkos::View<double*> m_inv_cell_nb_nodes;

  Kokkos::View<unsigned short**> m_face_cell_local_face;
  Kokkos::View<unsigned short**> m_face_node_local_face;

  size_t  m_max_nb_node_per_cell;

  struct Face
  {
    const unsigned int m_node0_id;
    const unsigned int m_node1_id;

    KOKKOS_INLINE_FUNCTION
    bool operator<(const Face& 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
    Face& operator=(const Face&) = default;

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

    KOKKOS_INLINE_FUNCTION
    Face(const Face&) = default;

    KOKKOS_INLINE_FUNCTION
    Face(Face&&) = default;

    KOKKOS_INLINE_FUNCTION
    Face(unsigned int node0_id,
         unsigned int node1_id)
        : m_node0_id(node0_id),
          m_node1_id(node1_id)
    {
      ;
    }

    KOKKOS_INLINE_FUNCTION
    ~Face() = default;
  };

  void _computeFaceCellConnectivities()
  {
    // In 2D faces are simply define
    typedef std::pair<unsigned int, unsigned short> CellFaceId;
    std::map<Face, std::vector<CellFaceId>> face_cells_map;
    for (unsigned int j=0; j<this->numberOfCells(); ++j) {
      const auto& cell_nodes = m_cell_to_node_matrix.rowConst(j);
      for (unsigned short r=0; r<cell_nodes.length; ++r) {
        unsigned int node0_id = cell_nodes(r);
        unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
        if (node1_id<node0_id) {
          std::swap(node0_id, node1_id);
        }
        face_cells_map[Face(node0_id, node1_id)].push_back(std::make_pair(j, r));
      }
    }

    {
      std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
      int l=0;
      for (const auto& face_info : face_cells_map) {
        const Face& face = face_info.first;
        face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
        ++l;
      }
      m_face_to_node_matrix
          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_node_matrix", face_to_node_vector);
    }

    {
      std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
      int l=0;
      for (const auto& face_cells_vector : face_cells_map) {
        const auto& [face, cell_info_vector] = face_cells_vector;
        for (const auto& cell_info : cell_info_vector) {
          face_to_cell_vector[l].push_back(cell_info.second);
        }
        ++l;
      }
      m_face_to_cell_matrix
          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("face_to_cell_matrix", face_to_cell_vector);
    }

    Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
                                                        face_cells_map.size(), m_max_nb_node_per_cell);
    {
      int l=0;
      for (const auto& face_cells_vector : face_cells_map) {
        const auto& cells_vector = face_cells_vector.second;
        for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
          unsigned short cell_local_face = cells_vector[lj].second;
          face_cell_local_face(l,lj) = cell_local_face;
        }
        ++l;
      }
    }
    m_face_cell_local_face = face_cell_local_face;
  }

 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 size_t& maxNbNodePerCell() const
  {
    return m_max_nb_node_per_cell;
  }

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

  const Kokkos::View<const unsigned short**> faceCellLocalFace() const
  {
    return m_face_cell_local_face;
  }

  unsigned int getFaceNumber(const unsigned int node0_id,
                             const unsigned int node1_id) const
  {
#warning rewrite
    const unsigned int n0_id = std::min(node0_id, node1_id);
    const unsigned int n1_id = std::max(node0_id, node1_id);
    // worst code ever
    for (unsigned int l=0; l<this->numberOfFaces(); ++l) {
      const auto& face_nodes = m_face_to_node_matrix.rowConst(l);
      if ((face_nodes(0) == n0_id) and (face_nodes(1) == n1_id)) {
        return l;
      }
    }
    std::cerr << "Face not found!\n";
    std::exit(0);
    return 0;
  }

  Connectivity2D(const Connectivity2D&) = delete;

  Connectivity2D(std::vector<std::vector<unsigned int>> cell_by_node_vector)
  {
    Assert(cell_by_node_vector.size()>0);

    m_cell_to_node_matrix
        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("cell_to_node_matrix",
                                                            cell_by_node_vector);

    {
      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;
    }

    {
      std::vector<std::vector<unsigned int>> node_id_per_cell_vector(this->numberOfCells());
      unsigned int id=0;
      for (unsigned int j=0; j<this->numberOfCells(); ++j) {
        const auto& cell_to_node = m_cell_to_node_matrix.rowConst(j);
        auto& node_id_per_cell = node_id_per_cell_vector[j];
        node_id_per_cell.resize(cell_to_node.length);
        for (size_t r=0; r<cell_to_node.length; ++r) {
          node_id_per_cell[r] = id++;
        }
      }
      m_node_id_per_cell_matrix
          = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_id_per_cell_matrix",
                                                              cell_by_node_vector);
    }

    ConnectivityUtils utils;
    utils.computeNodeCellConnectivity(m_cell_to_node_matrix,
                                      m_max_nb_node_per_cell,
                                      m_node_to_cell_matrix,
                                      m_node_to_cell_local_node_matrix);

    this->_computeFaceCellConnectivities();
  }

  ~Connectivity2D()
  {
    ;
  }
};

#endif // CONNECTIVITY_2D_HPP