ConnectivityUtils.hpp

#ifndef CONNECTIVITY_UTILS_HPP
#define CONNECTIVITY_UTILS_HPP

#include <map>
#include <Kokkos_Core.hpp>
#include <Kokkos_StaticCrsGraph.hpp>

typedef Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace> ConnectivityMatrix;

class ConnectivityUtils
{
 public:
  void computeNodeCellConnectivity(const Kokkos::View<const unsigned int**> cell_nodes,
                                   const Kokkos::View<const unsigned short*> cell_nb_nodes,
                                   const size_t& number_of_cells,
                                   size_t& max_nb_node_per_cell,
                                   size_t& number_of_nodes,
                                   Kokkos::View<const unsigned short*>& node_nb_cells,
                                   Kokkos::View<const unsigned int**>& node_cells,
                                   Kokkos::View<const unsigned short**>& node_cell_local_node,
                                   ConnectivityMatrix& node_to_cell_matrix,
                                   ConnectivityMatrix& node_to_cell_local_node_matrix)
  {
    std::map<unsigned int, std::vector<unsigned int>> node_cells_map;
    using namespace Kokkos::Experimental;
    Kokkos::parallel_reduce(number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
        const size_t n = cell_nb_nodes[j];
        if (n > nb_max) nb_max = n;
      }, Kokkos::Max<size_t>(max_nb_node_per_cell));

    for (unsigned int j=0; j<number_of_cells; ++j) {
      for (unsigned int r=0; r<cell_nb_nodes[j]; ++r) {
        node_cells_map[cell_nodes(j,r)].push_back(j);
      }
    }

    {
      size_t i=0;
      for (const auto& i_cell_vector : node_cells_map) {
        const auto& [node_id, cell_vector] = i_cell_vector;
        if (node_id != i) {
          std::cerr << "sparse node numerotation NIY\n";
          std::exit(0);
        }
        ++i;
      }
    }

    number_of_nodes = node_cells_map.size();

    Kokkos::View<unsigned short*> built_node_nb_cells("node_nb_cells", node_cells_map.size());
    size_t max_node_cells = 0;
    for (const auto& i_cell_vector : node_cells_map) {
      const auto& [i, cells_vector] = i_cell_vector;
      const size_t nb_cells = cells_vector.size();
      built_node_nb_cells[i] = nb_cells;
      if (nb_cells > max_node_cells) {
        max_node_cells = nb_cells;
      }
    }
    node_nb_cells = built_node_nb_cells;

    Kokkos::View<unsigned int**> built_node_cells("node_cells", node_cells_map.size(), max_node_cells);
    for (const auto& i_cell_vector : node_cells_map) {
      const auto& [i, cells_vector] = i_cell_vector;
      for (size_t j=0; j<cells_vector.size(); ++j) {
        built_node_cells(i,j) = cells_vector[j];
      }
    }
    node_cells = built_node_cells;

    std::vector<std::vector<unsigned int>> node_to_cell_vector(node_cells_map.size());
    for (const auto& i_cell_vector : node_cells_map) {
      const auto& [r, cells_vector] = i_cell_vector;
      node_to_cell_vector[r] = cells_vector;
    }
    node_to_cell_matrix
        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_to_cell_matrix", node_to_cell_vector);

    Kokkos::View<unsigned short**> built_node_cell_local_node("node_cell_local_node",
                                                              node_cells_map.size(), max_node_cells);
    Kokkos::parallel_for(number_of_nodes, KOKKOS_LAMBDA(const unsigned int& r){
        for (unsigned short J=0; J<node_nb_cells[r]; ++J) {
          const unsigned int j = node_cells(r,J);
          for (unsigned int R=0; R<cell_nb_nodes[j]; ++R) {
            if (cell_nodes(j,R) == r) {
              built_node_cell_local_node(r,J)=R;
              break;
            }
          }
        }
      });
    node_cell_local_node = built_node_cell_local_node;

    std::vector<std::vector<unsigned int>> node_to_cell_local_node_vector(node_cells_map.size());
    for (unsigned int r=0; r<node_cells_map.size(); ++r) {
      node_to_cell_local_node_vector[r].resize(node_nb_cells[r]);
      for (unsigned short J=0; J<node_nb_cells[r]; ++J) {
        const unsigned int j = node_cells(r,J);
        for (unsigned int R=0; R<cell_nb_nodes[j]; ++R) {
          if (cell_nodes(j,R) == r) {
            node_to_cell_local_node_vector[r][J] = R;
            break;
          }
        }
      }
    }
    node_to_cell_local_node_matrix
        = Kokkos::create_staticcrsgraph<ConnectivityMatrix>("node_to_cell_local_node_matrix", node_to_cell_local_node_vector);

  }
};

#endif // CONNECTIVITY_UTILS_HPP