From 0f44ba245ab5f263c88d4ec6fa24df3a9cdefe88 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Thu, 9 Aug 2018 16:55:28 +0200
Subject: [PATCH] Compute reverse connectivity on demand

- should compute some base connectivities on demand (eg. edge->cell,...)
- remove old connectivity matrices access interface
---
 src/mesh/Connectivity.cpp         | 23 ++------
 src/mesh/Connectivity.hpp         | 30 ++++++++---
 src/mesh/ConnectivityComputer.cpp | 87 ++++++++++++++++++++++---------
 src/mesh/ConnectivityComputer.hpp | 19 +++++--
 src/mesh/ConnectivityMatrix.hpp   | 10 +++-
 src/mesh/ConnectivityUtils.hpp    | 16 ------
 src/mesh/MeshData.hpp             | 28 +++++-----
 src/mesh/MeshNodeBoundary.hpp     | 33 ++++++------
 src/mesh/SubItemValuePerItem.hpp  |  1 -
 src/output/VTKWriter.hpp          |  6 +--
 src/scheme/AcousticSolver.hpp     | 20 +++----
 11 files changed, 154 insertions(+), 119 deletions(-)
 delete mode 100644 src/mesh/ConnectivityUtils.hpp

diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 5ce0c5aab..93354d6de 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -4,13 +4,12 @@
 template<>
 void Connectivity<3>::_computeFaceCellConnectivities()
 {
-  Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
-
   typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
 
   const auto& cell_to_node_matrix
-      = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
+      = this->getMatrix(ItemType::cell, ItemType::node);
 
+  Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
   std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
   for (unsigned int j=0; j<this->numberOfCells(); ++j) {
     const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -151,16 +150,6 @@ void Connectivity<3>::_computeFaceCellConnectivities()
     }
   }
 
-  const auto& cell_to_face_matrix
-      = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::face)];
-  auto& face_to_cell_matrix
-      = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::cell)];
-  m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_face_matrix,
-                                                           face_to_cell_matrix);
-
-  m_face_local_numbers_in_their_cells
-      = m_connectivity_computer.computeLocalItemNumberInChildItem<ItemType::face,
-                                                                  ItemType::cell>(*this);
 #warning check that the number of cell per faces is <=2
 }
 
@@ -168,7 +157,7 @@ template<>
 void Connectivity<2>::_computeFaceCellConnectivities()
 {
   const auto& cell_to_node_matrix
-      = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
+      = this->getMatrix(ItemType::cell, ItemType::node);
 
   // In 2D faces are simply define
   typedef std::pair<unsigned int, unsigned short> CellFaceId;
@@ -243,12 +232,6 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector)
     m_inv_cell_nb_nodes = inv_cell_nb_nodes;
   }
 
-  auto& node_to_cell_matrix
-      = m_item_to_item_matrix[itemId(ItemType::node)][itemId(ItemType::cell)];
-
-  m_connectivity_computer.computeInverseConnectivityMatrix(cell_to_node_matrix,
-                                                           node_to_cell_matrix);
-
   if constexpr (Dimension>1) {
     this->_computeFaceCellConnectivities();
   }
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 420dbb8d0..6647cecb9 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -275,16 +275,29 @@ class Connectivity final
   }
 
  public:
-  template <ItemType item_type_0, ItemType item_type_1>
-  const ConnectivityMatrix& itemToItemMatrix() const
+
+  KOKKOS_INLINE_FUNCTION
+  const bool& isConnectivityMatrixBuilt(const ItemType& item_type_0,
+                                        const ItemType& item_type_1) const
   {
-    return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
-  };
+    const ConnectivityMatrix& connectivity_matrix
+        = m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
+    return connectivity_matrix.isBuilt();
+  }
 
+  KOKKOS_INLINE_FUNCTION
   const ConnectivityMatrix& getMatrix(const ItemType& item_type_0,
                                       const ItemType& item_type_1) const
   {
-    return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
+    const ConnectivityMatrix& connectivity_matrix
+        = m_item_to_item_matrix[itemId(item_type_0)][itemId(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;
   }
 
   const auto& cellFaceIsReversed() const
@@ -352,7 +365,7 @@ class Connectivity final
   size_t numberOfNodes() const
   {
     const auto& node_to_cell_matrix
-        = m_item_to_item_matrix[itemId(ItemType::node)][itemId(ItemType::cell)];
+        = this->getMatrix(ItemType::node,ItemType::cell);
     return node_to_cell_matrix.numRows();
   }
 
@@ -360,7 +373,7 @@ class Connectivity final
   size_t numberOfFaces() const
   {
     const auto& face_to_node_matrix
-        = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::node)];
+        = this->getMatrix(ItemType::face,ItemType::cell);
     return face_to_node_matrix.numRows();
   }
 
@@ -368,12 +381,13 @@ class Connectivity final
   size_t numberOfCells() const
   {
     const auto& cell_to_node_matrix
-        = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
+        = this->getMatrix(ItemType::cell,ItemType::node);
     return cell_to_node_matrix.numRows();
   }
 
   const Kokkos::View<const double*> invCellNbNodes() const
   {
+#warning add calculation on demand when variables will be defined
     return m_inv_cell_nb_nodes;
   }
 
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index a128d1ed1..336c88170 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -1,26 +1,61 @@
 #include <Connectivity.hpp>
-#include <ConnectivityUtils.hpp>
 
 #include <ConnectivityComputer.hpp>
 #include <map>
 
-void ConnectivityComputer::
-computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_matrix,
-                                 ConnectivityMatrix& child_item_to_item_matrix) const
+
+template <typename ConnectivityType>
+KOKKOS_INLINE_FUNCTION
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const ConnectivityType& connectivity,
+                          const ItemType& item_type,
+                          const ItemType& child_item_type) const
 {
-  std::map<unsigned int, std::vector<unsigned int>> child_item_to_item_vector_map;
-  const size_t& number_of_items = item_to_child_item_matrix.numRows();
+  ConnectivityMatrix item_to_child_item_matrix;
+  if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
+    const ConnectivityMatrix& child_to_item_matrix
+        = connectivity.getMatrix(child_item_type, item_type);
+
+    item_to_child_item_matrix
+        = this->_computeInverse(child_to_item_matrix);
+  } else {
+    std::cerr << "unable to compute connectivity " << '\n';
+    std::exit(0);
+  }
+
+  return item_to_child_item_matrix;
+}
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity1D&, const ItemType&, const ItemType&) const;
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity2D&, const ItemType&, const ItemType&) const;
+
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity3D&, const ItemType&, const ItemType&) const;
+
+ConnectivityMatrix
+ConnectivityComputer::
+_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
+{
+  std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map;
+  const size_t& number_of_items = item_to_child_matrix.numRows();
 
   for (unsigned int j=0; j<number_of_items; ++j) {
-    const auto& item_to_child_items = item_to_child_item_matrix.rowConst(j);
-    for (unsigned int r=0; r<item_to_child_items.length; ++r) {
-      child_item_to_item_vector_map[item_to_child_items(r)].push_back(j);
+    const auto& item_to_child = item_to_child_matrix.rowConst(j);
+    for (unsigned int r=0; r<item_to_child.length; ++r) {
+      child_to_item_vector_map[item_to_child(r)].push_back(j);
     }
   }
 
   {
     size_t i=0;
-    for (const auto& [child_item_id, item_vector] : child_item_to_item_vector_map) {
+    for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
       if (child_item_id != i) {
         std::cerr << "sparse item numerotation NIY\n";
         std::exit(0);
@@ -29,11 +64,12 @@ computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_ma
     }
   }
 
-  std::vector<std::vector<unsigned int>> child_item_to_items_vector(child_item_to_item_vector_map.size());
-  for (const auto& [child_item_id, item_vector] : child_item_to_item_vector_map) {
-    child_item_to_items_vector[child_item_id] = item_vector;
+  std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size());
+  for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
+    child_to_items_vector[child_item_id] = item_vector;
   }
-  child_item_to_item_matrix = child_item_to_items_vector;
+
+  return ConnectivityMatrix(child_to_items_vector);
 }
 
 template <ItemType item_type,
@@ -43,9 +79,9 @@ SubItemValuePerItem<const unsigned short, child_item_type, item_type>
 ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
 {
   const ConnectivityMatrix& child_item_to_items_matrix
-      = getConnectivityMatrix<child_item_type, item_type>(connectivity);
+      = connectivity.getMatrix(child_item_type, item_type);
   const ConnectivityMatrix& item_to_child_items_matrix
-      = getConnectivityMatrix<item_type, child_item_type>(connectivity);
+      = connectivity.getMatrix(item_type, child_item_type);
 
   SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
   for (unsigned int r=0; r<item_to_child_items_matrix.numRows(); ++r) {
@@ -71,49 +107,48 @@ ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType&
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::node,
-                                  ItemType::cell>(const Connectivity1D& connectivity) const;
+                                  ItemType::cell>(const Connectivity1D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::face,
-                                  ItemType::cell>(const Connectivity1D& connectivity) const;
+                                  ItemType::cell>(const Connectivity1D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::cell,
-                                  ItemType::node>(const Connectivity1D& connectivity) const;
+                                  ItemType::node>(const Connectivity1D&) const;
 
 // 2D
 
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::node,
-                                  ItemType::cell>(const Connectivity2D& connectivity) const;
+                                  ItemType::cell>(const Connectivity2D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::face,
-                                  ItemType::cell>(const Connectivity2D& connectivity) const;
+                                  ItemType::cell>(const Connectivity2D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::cell,
-                                  ItemType::node>(const Connectivity2D& connectivity) const;
+                                  ItemType::node>(const Connectivity2D&) const;
 
 // 3D
 
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::node,
-                                  ItemType::cell>(const Connectivity3D& connectivity) const;
-
+                                  ItemType::cell>(const Connectivity3D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::face,
-                                  ItemType::cell>(const Connectivity3D& connectivity) const;
+                                  ItemType::cell>(const Connectivity3D&) const;
 
 template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
 ConnectivityComputer::
 computeLocalItemNumberInChildItem<ItemType::cell,
-                                  ItemType::node>(const Connectivity3D& connectivity) const;
+                                  ItemType::node>(const Connectivity3D&) const;
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
index 6a726e52d..8b4fdfc3a 100644
--- a/src/mesh/ConnectivityComputer.hpp
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -4,16 +4,29 @@
 #include <ConnectivityMatrix.hpp>
 #include <SubItemValuePerItem.hpp>
 
-struct ConnectivityComputer
+class ConnectivityComputer
 {
-  void computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_matrix,
-                                        ConnectivityMatrix& child_item_to_item_matrix) const;
+ private:
+  ConnectivityMatrix
+  _computeInverse(const ConnectivityMatrix& item_to_child_matrix) const;
+
+ public:
+  template <typename ConnectivityType>
+  ConnectivityMatrix
+  computeConnectivityMatrix(const ConnectivityType& connectivity,
+                            const ItemType& item_type,
+                            const ItemType& child_item_type) const;
+
 
   template <ItemType item_type,
             ItemType child_item_type,
             typename ConnectivityType>
   SubItemValuePerItem<const unsigned short, child_item_type, item_type>
   computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const;
+
+  ConnectivityComputer(const ConnectivityComputer&) = default;
+  ConnectivityComputer() = default;
+  ~ConnectivityComputer() = default;
 };
 
 #endif // CONNECTIVITY_COMPUTER_HPP
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 30a40e180..1a0ab4880 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -10,9 +10,16 @@ class ConnectivityMatrix
   typedef Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace> HostMatrix;
   HostMatrix m_host_matrix;
 
+  bool m_is_built{false};
+
  public:
   typedef HostMatrix::row_map_type HostRowType;
 
+  const bool& isBuilt() const
+  {
+    return m_is_built;
+  }
+
   const auto& entries() const
   {
     return m_host_matrix.entries;
@@ -49,7 +56,8 @@ class ConnectivityMatrix
 
   KOKKOS_INLINE_FUNCTION
   ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer)
-      : m_host_matrix(Kokkos::create_staticcrsgraph<HostMatrix>("connecticity_matrix", initializer))
+      : m_host_matrix{Kokkos::create_staticcrsgraph<HostMatrix>("connecticity_matrix", initializer)},
+        m_is_built{true}
   {
     ;
   }
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
deleted file mode 100644
index e073d5521..000000000
--- a/src/mesh/ConnectivityUtils.hpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#ifndef CONNECTIVITY_UTILS_HPP
-#define CONNECTIVITY_UTILS_HPP
-
-#include <ItemType.hpp>
-#include <ConnectivityMatrix.hpp>
-
-template <ItemType item_type_0,
-          ItemType item_type_1,
-          typename ConnectivityType>
-inline
-const ConnectivityMatrix& getConnectivityMatrix(const ConnectivityType& connectivity)
-{
-  return connectivity.template itemToItemMatrix<item_type_0, item_type_1>();
-}
-
-#endif // CONNECTIVITY_UTILS_HPP
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 66ba3626a..0fe98256b 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -36,8 +36,8 @@ class MeshData
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
       const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                            ItemType::node);
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -49,8 +49,8 @@ class MeshData
           = m_mesh.connectivity().invCellNbNodes();
 
       const auto& cell_to_node_matrix
-          = getConnectivityMatrix<ItemType::cell,
-                                  ItemType::node>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                            ItemType::node);
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           Rd X = zero;
@@ -68,8 +68,8 @@ class MeshData
   {
     const Kokkos::View<const Rd*> xr = m_mesh.xr();
     const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(m_mesh.connectivity());
+        = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                          ItemType::node);
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
         double sum_cjr_xr = 0;
@@ -90,8 +90,8 @@ class MeshData
     else if constexpr (dimension == 2) {
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
       const auto& cell_to_node_matrix
-          = getConnectivityMatrix<ItemType::cell,
-                                  ItemType::node>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                            ItemType::node);
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -118,8 +118,8 @@ class MeshData
 
       NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
       const auto& face_to_node_matrix
-          = getConnectivityMatrix<ItemType::face,
-                                  ItemType::node>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::face,
+                                            ItemType::node);
 
       Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
           const auto& face_nodes = face_to_node_matrix.rowConst(l);
@@ -146,12 +146,12 @@ class MeshData
         });
 
       const auto& cell_to_node_matrix
-          = getConnectivityMatrix<ItemType::cell,
-                                  ItemType::node>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                            ItemType::node);
 
       const auto& cell_to_face_matrix
-          = getConnectivityMatrix<ItemType::cell,
-                                  ItemType::face>(m_mesh.connectivity());
+          = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                            ItemType::face);
 
       const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index cf9df04f0..87c537439 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -9,7 +9,7 @@
 #include <RefFaceList.hpp>
 
 #include <ConnectivityMatrix.hpp>
-#include <ConnectivityUtils.hpp>
+#include <IConnectivity.hpp>
 #include <iostream>
 
 template <size_t dimension>
@@ -32,8 +32,7 @@ class MeshNodeBoundary
   {
     static_assert(dimension == MeshType::dimension);
     const auto& face_to_cell_matrix
-        = getConnectivityMatrix<ItemType::face,
-                                ItemType::cell>(mesh.connectivity());
+        = mesh.connectivity().getMatrix(ItemType::face, ItemType::cell);
 
     const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
     Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
@@ -48,8 +47,8 @@ class MeshNodeBoundary
     // not enough but should reduce significantly the number of resizing
     node_ids.reserve(dimension*face_list.extent(0));
     const auto& face_to_node_matrix
-        = getConnectivityMatrix<ItemType::face,
-                                ItemType::node>(mesh.connectivity());
+        = mesh.connectivity().getMatrix(ItemType::face,
+                                        ItemType::node);
 
     for (size_t l=0; l<face_list.extent(0); ++l) {
       const size_t face_number = face_list[l];
@@ -317,12 +316,12 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const Kokkos::View<const R*>& xr = mesh.xr();
   const auto& cell_to_node_matrix
-      = getConnectivityMatrix<ItemType::cell,
-                              ItemType::node>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::cell,
+                                      ItemType::node);
 
   const auto& node_to_cell_matrix
-          = getConnectivityMatrix<ItemType::node,
-                                  ItemType::cell>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::node,
+                                      ItemType::cell);
 
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
@@ -354,12 +353,12 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const Kokkos::View<const R2*>& xr = mesh.xr();
   const auto& cell_to_node_matrix
-      = getConnectivityMatrix<ItemType::cell,
-                              ItemType::node>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::cell,
+                                      ItemType::node);
 
   const auto& node_to_cell_matrix
-          = getConnectivityMatrix<ItemType::node,
-                                  ItemType::cell>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::node,
+                                      ItemType::cell);
 
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
@@ -391,12 +390,12 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const Kokkos::View<const R3*>& xr = mesh.xr();
   const auto& cell_to_node_matrix
-      = getConnectivityMatrix<ItemType::cell,
-                              ItemType::node>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::cell,
+                                      ItemType::node);
 
   const auto& node_to_cell_matrix
-          = getConnectivityMatrix<ItemType::node,
-                                  ItemType::cell>(mesh.connectivity());
+      = mesh.connectivity().getMatrix(ItemType::node,
+                                      ItemType::cell);
 
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index c6e98dcf6..3ab0ac487 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -7,7 +7,6 @@
 
 #include <ConnectivityMatrix.hpp>
 #include <IConnectivity.hpp>
-#include <ConnectivityUtils.hpp>
 
 template <typename DataType,
           ItemType sub_item_type,
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index a2454cb0f..a33fef0bf 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -6,7 +6,7 @@
 #include <iomanip>
 #include <sstream>
 #include <TinyVector.hpp>
-#include <ConnectivityUtils.hpp>
+#include <IConnectivity.hpp>
 
 class VTKWriter
 {
@@ -72,8 +72,8 @@ class VTKWriter
     fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
 
     const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(mesh.connectivity());
+        = mesh.connectivity().getMatrix(ItemType::cell,
+                                        ItemType::node);
 
     for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
       const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 4145c01ea..0cc6ec0df 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -99,8 +99,8 @@ class AcousticSolver
   const Kokkos::View<const Rdd*>
   computeAr(const NodeValuePerCell<Rdd>& Ajr) {
     const auto& node_to_cell_matrix
-        = getConnectivityMatrix<ItemType::node,
-                                ItemType::cell>(m_connectivity);
+        = m_connectivity.getMatrix(ItemType::node,
+                                   ItemType::cell);
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
@@ -128,8 +128,8 @@ class AcousticSolver
             const Kokkos::View<const Rd*>& uj,
             const Kokkos::View<const double*>& pj) {
     const auto& node_to_cell_matrix
-        = getConnectivityMatrix<ItemType::node,
-                                ItemType::cell>(m_connectivity);
+        = m_connectivity.getMatrix(ItemType::node,
+                                   ItemType::cell);
     const auto& node_local_numbers_in_their_cells
         = m_connectivity.nodeLocalNumbersInTheirCells();
 
@@ -212,8 +212,8 @@ class AcousticSolver
              const Kokkos::View<const double*>& pj)
   {
     const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(m_mesh.connectivity());
+        = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                          ItemType::node);
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -291,8 +291,8 @@ class AcousticSolver
   {
     const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
     const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(m_mesh.connectivity());
+        = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                          ItemType::node);
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -335,8 +335,8 @@ class AcousticSolver
     const NodeValuePerCell<Rd>& Fjr = m_Fjr;
     const Kokkos::View<const Rd*> ur = m_ur;
     const auto& cell_to_node_matrix
-        = getConnectivityMatrix<ItemType::cell,
-                                ItemType::node>(m_mesh.connectivity());
+        = m_mesh.connectivity().getMatrix(ItemType::cell,
+                                          ItemType::node);
 
     const Kokkos::View<const double*> inv_mj = unknowns.invMj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
-- 
GitLab