Connectivity.hpp

#ifndef CONNECTIVITY_HPP
#define CONNECTIVITY_HPP

#include <algebra/TinyVector.hpp>
#include <mesh/CellType.hpp>
#include <mesh/ConnectivityComputer.hpp>
#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/IConnectivity.hpp>
#include <mesh/ItemToItemMatrix.hpp>
#include <mesh/ItemType.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/RefId.hpp>
#include <mesh/RefItemList.hpp>
#include <mesh/SubItemValuePerItem.hpp>
#include <utils/CRSGraph.hpp>
#include <utils/Exceptions.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
#include <utils/PugsTraits.hpp>
#include <utils/PugsUtils.hpp>

#include <algorithm>
#include <iostream>
#include <set>
#include <vector>

class ConnectivityDescriptor;

template <size_t Dim>
class Connectivity final : public IConnectivity
{
 public:
  PUGS_INLINE
  static std::shared_ptr<Connectivity<Dim>> build(const ConnectivityDescriptor&);

 private:
  constexpr static auto& itemTId = ItemTypeId<Dim>::itemTId;

 public:
  static constexpr size_t Dimension = Dim;

  PUGS_INLINE
  size_t
  dimension() const final
  {
    return Dimension;
  }

 private:
  size_t m_number_of_cells;
  size_t m_number_of_faces;
  size_t m_number_of_edges;
  size_t m_number_of_nodes;

  ConnectivityMatrix m_item_to_item_matrix[Dimension + 1][Dimension + 1];
  WeakCellValue<const CellType> m_cell_type;

  WeakCellValue<const int> m_cell_global_index;

  WeakCellValue<const int> m_cell_number;
  WeakFaceValue<const int> m_face_number;
  WeakEdgeValue<const int> m_edge_number;
  WeakNodeValue<const int> m_node_number;

  WeakCellValue<const int> m_cell_owner;
  WeakCellValue<const bool> m_cell_is_owned;

  WeakFaceValue<const int> m_face_owner;
  WeakFaceValue<const bool> m_face_is_owned;
  WeakFaceValue<const bool> m_is_boundary_face;

  WeakEdgeValue<const int> m_edge_owner;
  WeakEdgeValue<const bool> m_edge_is_owned;
  WeakEdgeValue<const bool> m_is_boundary_edge;

  WeakNodeValue<const int> m_node_owner;
  WeakNodeValue<const bool> m_node_is_owned;
  WeakNodeValue<const bool> m_is_boundary_node;

  WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
  WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;

  WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
  WeakFaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces;
  WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;

  WeakCellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells;
  WeakEdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges;
  WeakNodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes;

  WeakCellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells;
  WeakFaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces;
  WeakNodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes;

  WeakCellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells;
  WeakFaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces;
  WeakEdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges;

  ConnectivityComputer m_connectivity_computer;

  std::vector<RefCellList> m_ref_cell_list_vector;
  std::vector<RefFaceList> m_ref_face_list_vector;
  std::vector<RefEdgeList> m_ref_edge_list_vector;
  std::vector<RefNodeList> m_ref_node_list_vector;

  WeakCellValue<const double> m_inv_cell_nb_nodes;

  void _computeCellFaceAndFaceNodeConnectivities();

  void _buildIsBoundaryFace() const;
  void _buildIsBoundaryEdge() const;
  void _buildIsBoundaryNode() const;

  template <typename SubItemValuePerItemType>
  PUGS_INLINE const SubItemValuePerItemType&
  _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
  {
    using ReversedItemOfItem = typename SubItemValuePerItemType::ItemOfItemType::Reversed;
    if (not sub_item_value_per_item.isBuilt()) {
      const_cast<SubItemValuePerItemType&>(sub_item_value_per_item) =
        m_connectivity_computer.computeLocalItemNumberInChildItem<ReversedItemOfItem>(*this);
    }
    return sub_item_value_per_item;
  }

  friend class ConnectivityComputer;

  PUGS_INLINE
  const ConnectivityMatrix&
  _getMatrix(const ItemType& item_type_0, const ItemType& item_type_1) const
  {
    const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
    if (not connectivity_matrix.isBuilt()) {
      const_cast<ConnectivityMatrix&>(connectivity_matrix) =
        m_connectivity_computer.computeConnectivityMatrix(*this, item_type_0, item_type_1);
    }
    return connectivity_matrix;
  }

 public:
  PUGS_INLINE
  CellValue<const CellType>
  cellType() const
  {
    return m_cell_type;
  }

  PUGS_INLINE
  CellValue<const int>
  cellNumber() const
  {
    return m_cell_number;
  }

  PUGS_INLINE
  FaceValue<const int>
  faceNumber() const
  {
    return m_face_number;
  }

  PUGS_INLINE
  EdgeValue<const int>
  edgeNumber() const
  {
    return m_edge_number;
  }

  PUGS_INLINE
  NodeValue<const int>
  nodeNumber() const
  {
    return m_node_number;
  }

  template <ItemType item_type>
  PUGS_INLINE auto
  number() const
  {
    if constexpr (item_type == ItemType::cell) {
      return m_cell_number;
    } else if constexpr (item_type == ItemType::face) {
      return m_face_number;
    } else if constexpr (item_type == ItemType::edge) {
      return m_edge_number;
    } else if constexpr (item_type == ItemType::node) {
      return m_node_number;
    } else {
      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
    }
  }

  PUGS_INLINE
  CellValue<const int>
  cellOwner() const
  {
    return m_cell_owner;
  }

  PUGS_INLINE
  FaceValue<const int>
  faceOwner() const
  {
    return m_face_owner;
  }

  PUGS_INLINE
  EdgeValue<const int>
  edgeOwner() const
  {
    return m_edge_owner;
  }

  PUGS_INLINE
  NodeValue<const int>
  nodeOwner() const
  {
    return m_node_owner;
  }

  template <ItemType item_type>
  PUGS_INLINE const auto&
  owner() const
  {
    if constexpr (item_type == ItemType::cell) {
      return m_cell_owner;
    } else if constexpr (item_type == ItemType::face) {
      return m_face_owner;
    } else if constexpr (item_type == ItemType::edge) {
      return m_edge_owner;
    } else if constexpr (item_type == ItemType::node) {
      return m_node_owner;
    } else {
      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
    }
  }

  PUGS_INLINE
  CellValue<const bool>
  cellIsOwned() const
  {
    return m_cell_is_owned;
  }

  PUGS_INLINE
  FaceValue<const bool>
  faceIsOwned() const
  {
    return m_face_is_owned;
  }

  PUGS_INLINE
  EdgeValue<const bool>
  edgeIsOwned() const
  {
    return m_edge_is_owned;
  }

  PUGS_INLINE
  NodeValue<const bool>
  nodeIsOwned() const
  {
    return m_node_is_owned;
  }

  template <ItemType item_type>
  PUGS_INLINE auto
  isOwned() const
  {
    if constexpr (item_type == ItemType::cell) {
      return m_cell_is_owned;
    } else if constexpr (item_type == ItemType::face) {
      return m_face_is_owned;
    } else if constexpr (item_type == ItemType::edge) {
      return m_edge_is_owned;
    } else if constexpr (item_type == ItemType::node) {
      return m_node_is_owned;
    } else {
      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
    }
  }

  PUGS_INLINE
  FaceValue<const bool>
  isBoundaryFace() const
  {
    if (not m_is_boundary_face.isBuilt()) {
      this->_buildIsBoundaryFace();
    }
    return m_is_boundary_face;
  }

  PUGS_INLINE
  EdgeValue<const bool>
  isBoundaryEdge() const
  {
    if (not m_is_boundary_edge.isBuilt()) {
      this->_buildIsBoundaryEdge();
    }
    return m_is_boundary_edge;
  }

  PUGS_INLINE
  NodeValue<const bool>
  isBoundaryNode() const
  {
    if (not m_is_boundary_node.isBuilt()) {
      this->_buildIsBoundaryNode();
    }
    return m_is_boundary_node;
  }

  template <ItemType item_type>
  PUGS_INLINE auto
  isBoundary() const
  {
    if constexpr (item_type == ItemType::face) {
      return isBoundaryFace();
    } else if constexpr (item_type == ItemType::edge) {
      return isBoundaryEdge();
    } else if constexpr (item_type == ItemType::node) {
      return isBoundaryNode();
    } else {
      static_assert(item_type != ItemType::cell, "cell boundary makes no sense");
      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
    }
  }

  PUGS_INLINE
  const bool&
  isConnectivityMatrixBuilt(ItemType item_type_0, ItemType item_type_1) const
  {
    const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
    return connectivity_matrix.isBuilt();
  }

  template <ItemType source_item_type, ItemType target_item_type>
  PUGS_INLINE ItemToItemMatrix<source_item_type, target_item_type>
  getItemToItemMatrix() const
  {
    return ItemToItemMatrix<source_item_type, target_item_type>(_getMatrix(source_item_type, target_item_type));
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::cell, ItemType::face>
  cellToFaceMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::cell, ItemType::edge>
  cellToEdgeMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::cell, ItemType::node>
  cellToNodeMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::face, ItemType::cell>
  faceToCellMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::face, ItemType::edge>
  faceToEdgeMatrix() const
  {
    static_assert(Dimension > 2, "face to edge matrix makes sense in dimension > 2");
    return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::face, ItemType::node>
  faceToNodeMatrix() const
  {
    static_assert(Dimension > 1, "face to node matrix makes sense in dimension > 1");
    return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::edge, ItemType::cell>
  edgeToCellMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::edge, ItemType::face>
  edgeToFaceMatrix() const
  {
    static_assert(Dimension > 2, "edge to face matrix makes sense in dimension > 2");
    return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::edge, ItemType::node>
  edgeToNodeMatrix() const
  {
    static_assert(Dimension > 1, "edge to node matrix makes sense in dimension > 1");
    return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::node, ItemType::cell>
  nodeToCellMatrix() const
  {
    return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::node, ItemType::face>
  nodeToFaceMatrix() const
  {
    static_assert(Dimension > 1, "node to face matrix makes sense in dimension > 1");
    return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
  }

  PUGS_INLINE
  ItemToItemMatrix<ItemType::node, ItemType::edge>
  nodeToEdgeMatrix() const
  {
    static_assert(Dimension > 1, "node to edge matrix makes sense in dimension > 1");
    return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
  }

  PUGS_INLINE
  FaceValuePerCell<const bool>
  cellFaceIsReversed() const
  {
    static_assert(Dimension > 1, "reversed faces makes no sense in dimension 1");
    return m_cell_face_is_reversed;
  }

  PUGS_INLINE
  EdgeValuePerFace<const bool>
  faceEdgeIsReversed() const
  {
    static_assert(Dimension > 2, "reversed edges makes no sense in dimension 1 or 2");
    return m_face_edge_is_reversed;
  }

  PUGS_INLINE
  const auto&
  cellLocalNumbersInTheirNodes() const
  {
    return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
  }

  PUGS_INLINE
  const auto&
  cellLocalNumbersInTheirEdges() const
  {
    if constexpr (Dimension > 2) {
      return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges);
    } else {
      return cellLocalNumbersInTheirFaces();
    }
  }

  PUGS_INLINE
  const auto&
  cellLocalNumbersInTheirFaces() const
  {
    if constexpr (Dimension > 1) {
      return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces);
    } else {
      return cellLocalNumbersInTheirNodes();
    }
  }

  PUGS_INLINE
  const auto&
  faceLocalNumbersInTheirCells() const
  {
    if constexpr (Dimension > 1) {
      return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells);
    } else {
      return nodeLocalNumbersInTheirCells();
    }
  }

  PUGS_INLINE
  const auto&
  faceLocalNumbersInTheirEdges() const
  {
    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
    return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges);
  }

  PUGS_INLINE
  const auto&
  faceLocalNumbersInTheirNodes() const
  {
    static_assert(Dimension > 1, "this function has no meaning in 1d");
    return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes);
  }

  PUGS_INLINE
  const auto&
  edgeLocalNumbersInTheirCells() const
  {
    if constexpr (Dimension > 2) {
      return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells);
    } else {
      return faceLocalNumbersInTheirCells();
    }
  }

  PUGS_INLINE
  const auto&
  edgeLocalNumbersInTheirFaces() const
  {
    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
    return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces);
  }

  PUGS_INLINE
  const auto&
  edgeLocalNumbersInTheirNodes() const
  {
    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
    return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes);
  }

  PUGS_INLINE
  const auto&
  nodeLocalNumbersInTheirCells() const
  {
    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
  }

  PUGS_INLINE
  const auto&
  nodeLocalNumbersInTheirEdges() const
  {
    static_assert(Dimension > 2, "this function has no meaning in 1d or 2d");
    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges);
  }

  PUGS_INLINE
  const auto&
  nodeLocalNumbersInTheirFaces() const
  {
    static_assert(Dimension > 1, "this function has no meaning in 1d");
    return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
  }

  template <ItemType item_type>
  size_t
  numberOfRefItemList() const
  {
    if constexpr (item_type == ItemType::cell) {
      return m_ref_cell_list_vector.size();
    } else if constexpr (item_type == ItemType::face) {
      return m_ref_face_list_vector.size();
    } else if constexpr (item_type == ItemType::edge) {
      return m_ref_edge_list_vector.size();
    } else if constexpr (item_type == ItemType::node) {
      return m_ref_node_list_vector.size();
    } else {
      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
    }
  }

  template <ItemType item_type>
  const RefItemList<item_type>&
  refItemList(size_t i) const
  {
    if constexpr (item_type == ItemType::cell) {
      return m_ref_cell_list_vector[i];
    } else if constexpr (item_type == ItemType::face) {
      return m_ref_face_list_vector[i];
    } else if constexpr (item_type == ItemType::edge) {
      return m_ref_edge_list_vector[i];
    } else if constexpr (item_type == ItemType::node) {
      return m_ref_node_list_vector[i];
    } else {
      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
    }
  }

  template <ItemType item_type>
  void
  addRefItemList(const RefItemList<item_type>& ref_item_list)
  {
    if constexpr (item_type == ItemType::cell) {
      m_ref_cell_list_vector.push_back(ref_item_list);
    } else if constexpr (item_type == ItemType::face) {
      m_ref_face_list_vector.push_back(ref_item_list);
    } else if constexpr (item_type == ItemType::edge) {
      m_ref_edge_list_vector.push_back(ref_item_list);
    } else if constexpr (item_type == ItemType::node) {
      m_ref_node_list_vector.push_back(ref_item_list);
    } else {
      static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
    }
  }

  PUGS_INLINE
  CRSGraph
  cellToCellGraph() const
  {
    std::vector<std::set<int>> cell_cells(this->numberOfCells());
    const auto& face_to_cell_matrix = this->faceToCellMatrix();

    for (FaceId l = 0; l < this->numberOfFaces(); ++l) {
      const auto& face_to_cell = face_to_cell_matrix[l];
      if (face_to_cell.size() > 1) {
        const CellId cell_0 = face_to_cell[0];
        const CellId cell_1 = face_to_cell[1];

        cell_cells[cell_0].insert(cell_1);
        cell_cells[cell_1].insert(cell_0);
      }
    }

    Array<int> entries(this->numberOfCells() + 1);
    entries[0] = 0;
    for (size_t j = 0; j < this->numberOfCells(); ++j) {
      entries[j + 1] = entries[j] + cell_cells[j].size();
    }
    Array<int> neighbors(entries[this->numberOfCells()]);
    {
      size_t k = 0;
      for (size_t j = 0; j < this->numberOfCells(); ++j) {
        for (CellId cell_id : cell_cells[j]) {
          neighbors[k] = m_cell_global_index[cell_id];
          ++k;
        }
      }
    }
    return CRSGraph(entries, neighbors);
  }

  PUGS_INLINE
  size_t
  numberOfNodes() const final
  {
    return m_number_of_nodes;
  }

  PUGS_INLINE
  size_t
  numberOfEdges() const final
  {
    return m_number_of_edges;
  }

  PUGS_INLINE
  size_t
  numberOfFaces() const final
  {
    return m_number_of_faces;
  }

  PUGS_INLINE
  size_t
  numberOfCells() const final
  {
    return m_number_of_cells;
  }

  template <ItemType item_type>
  PUGS_INLINE size_t
  numberOf() const
  {
    if constexpr (item_type == ItemType::cell) {
      return this->numberOfCells();
    } else if constexpr (item_type == ItemType::face) {
      return this->numberOfFaces();
    } else if constexpr (item_type == ItemType::edge) {
      return this->numberOfEdges();
    } else if constexpr (item_type == ItemType::node) {
      return this->numberOfNodes();
    } else {
      static_assert(is_false_item_type_v<item_type>, "unknown ItemType");
    }
  }

  CellValue<const double>
  invCellNbNodes() const
  {
    if (not m_inv_cell_nb_nodes.isBuilt()) {
      const auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];

      WeakCellValue<double> inv_cell_nb_nodes(*this);
      parallel_for(
        this->numberOfCells(), PUGS_LAMBDA(CellId j) {
          const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
          inv_cell_nb_nodes[j]   = 1. / cell_nodes.length;
        });
      const_cast<WeakCellValue<const double>&>(m_inv_cell_nb_nodes) = inv_cell_nb_nodes;
    }

    return m_inv_cell_nb_nodes;
  }

  Connectivity(const Connectivity&) = delete;

 private:
  Connectivity();
  void _buildFrom(const ConnectivityDescriptor& descriptor);

 public:
  ~Connectivity() = default;
};

template <size_t Dimension>
PUGS_INLINE std::shared_ptr<Connectivity<Dimension>>
Connectivity<Dimension>::build(const ConnectivityDescriptor& descriptor)
{
  // Cannot use std::make_shared in this context since the default constructor is private
  std::shared_ptr<Connectivity<Dimension>> connectivity_ptr(new Connectivity<Dimension>);
  connectivity_ptr->_buildFrom(descriptor);

  return connectivity_ptr;
}

using Connectivity3D = Connectivity<3>;
using Connectivity2D = Connectivity<2>;
using Connectivity1D = Connectivity<1>;

#endif   // CONNECTIVITY_HPP