Connectivity2D.hpp

#ifndef CONNECTIVITY_2D_HPP
#define CONNECTIVITY_2D_HPP

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

#include <vector>
#include <map>
#include <algorithm>

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

private:
  const size_t  m_number_of_cells;
  size_t m_number_of_faces;
  size_t m_number_of_nodes;

  const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
  const Kokkos::View<const unsigned int**> m_cell_nodes;
  Kokkos::View<double*> m_inv_cell_nb_nodes;

  Kokkos::View<unsigned short*> m_cell_nb_faces;
  Kokkos::View<unsigned int**> m_cell_faces;

  Kokkos::View<unsigned short*> m_node_nb_cells;
  Kokkos::View<unsigned int**> m_node_cells;
  Kokkos::View<unsigned short**> m_node_cell_local_node;
  
  Kokkos::View<unsigned short*> m_face_nb_cells;
  Kokkos::View<unsigned int**> m_face_cells;
  Kokkos::View<unsigned short**> m_face_cell_local_face;

  size_t  m_max_nb_node_per_cell;

public:
  const size_t& numberOfNodes() const
  {
    return m_number_of_nodes;
  }

  const size_t& numberOfFaces() const
  {
    return m_number_of_faces;
  }

  const size_t& numberOfCells() const
  {
    return m_number_of_cells;
  }

  const size_t& maxNbNodePerCell() const
  {
    return m_max_nb_node_per_cell;
  }
  
  const Kokkos::View<const unsigned int**> cellNodes() const
  {
    return m_cell_nodes;
  }

  const Kokkos::View<const unsigned int**> cellFaces() const
  {
    return m_cell_faces;
  }

  const Kokkos::View<const unsigned short*> nodeNbCells() const
  {
    return m_node_nb_cells;
  }

  const Kokkos::View<const unsigned short*> cellNbNodes() const
  {
    return m_cell_nb_nodes;
  }

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

  const Kokkos::View<const unsigned short*> cellNbFaces() const
  {
    return m_cell_nb_faces;
  }

  const Kokkos::View<const unsigned short*> faceNbCells() const
  {
    return m_face_nb_cells;
  }

  const Kokkos::View<const unsigned int**> nodeCells() const
  {
    return m_node_cells;
  }

  const Kokkos::View<const unsigned int**> faceCells() const
  {
    return m_face_cells;
  }

  const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
  {
    return m_node_cell_local_node;
  }

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

  Connectivity2D(const Connectivity2D&) = delete;
  
  Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
		 const Kokkos::View<const unsigned int**> cell_nodes)
    : m_number_of_cells (cell_nodes.extent(0)),
      m_cell_nb_nodes(cell_nb_nodes),
      m_cell_nodes (cell_nodes)
  {
    {
      Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
      Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){
	  inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
	});
      m_inv_cell_nb_nodes = inv_cell_nb_nodes;
    }
    assert(m_number_of_cells>0);

    // Computes inefficiently node->cells connectivity [Version 0]
    
    std::multimap<unsigned int, unsigned int> node_cells_map;
    using namespace Kokkos::Experimental;
    Kokkos::parallel_reduce(m_number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
	const size_t n = m_cell_nb_nodes[j];
	if (n > nb_max) nb_max = n; 
      }, Max<size_t>(m_max_nb_node_per_cell));
    
    for (unsigned int j=0; j<m_number_of_cells; ++j) {
      for (unsigned int r=0; r<m_cell_nb_nodes[j]; ++r) {
	node_cells_map.insert(std::make_pair(cell_nodes(j,r),j));
      }
    }

    std::vector<unsigned int> node_ids;
    node_ids.reserve(node_cells_map.size());
    for (const auto& node_cell: node_cells_map) {
      node_ids.push_back(node_cell.first);
    }
    auto last_unique = std::unique(node_ids.begin(), node_ids.end());
    node_ids.resize(std::distance(node_ids.begin(), last_unique));

    m_number_of_nodes = node_ids.size();

    std::cout << "node_ids.size()=" << node_ids.size() << '\n';
    if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) {
      std::cerr << "sparse node numerotation NIY\n";
      for (int i=0; i<node_ids.size(); ++i) {
	std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
      }
      std::exit(0);
    }

    std::vector<std::vector<unsigned int>> node_cells_vector(node_ids.size());
    for (const auto& node_cell: node_cells_map) {
      node_cells_vector[node_cell.first].push_back(node_cell.second);
    }


    Kokkos::View<unsigned short*> node_nb_cells("node_nb_cells", node_ids.size());
    int max_node_cells = 0;
    for (int i=0; i<node_cells_vector.size(); ++i) {
      const auto& cells_vector = node_cells_vector[i];
      const size_t nb_cells = cells_vector.size();
      node_nb_cells[i] = nb_cells;
      if (nb_cells > max_node_cells) {
	max_node_cells = nb_cells;
      }
    }
    m_node_nb_cells = node_nb_cells;
    Kokkos::View<unsigned int**> node_cells("node_cells", node_ids.size(), max_node_cells);
    for (size_t i=0; i<node_cells_vector.size(); ++i) {
      const auto& cells_vector = node_cells_vector[i];
      for (size_t j=0; j<cells_vector.size(); ++j) {
	node_cells(i,j) = cells_vector[j];
      }
    }
    m_node_cells = node_cells;

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

  ~Connectivity2D()
  {
    ;
  }
};

#endif // CONNECTIVITY_2D_HPP