diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 5e1dc1b4da0860836f1e6e9b99ac6e2e293e0155..d9ad48e86a8a7edc006a8d063fd51092dc88b118 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -227,37 +227,6 @@ class Connectivity
  public:
   static constexpr size_t dimension = Dimension;
 
- public:
-  KOKKOS_INLINE_FUNCTION
-  ConnectivityMatrix cellToNodeMatrix() const
-  {
-    return m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::node)];
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  ConnectivityMatrix cellToFaceMatrix() const
-  {
-    return m_item_to_item_matrix[itemId(TypeOfItem::cell)][itemId(TypeOfItem::face)];
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  ConnectivityMatrix faceToCellMatrix() const
-  {
-    return m_item_to_item_matrix[itemId(TypeOfItem::face)][itemId(TypeOfItem::cell)];
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  ConnectivityMatrix faceToNodeMatrix() const
-  {
-    return m_item_to_item_matrix[itemId(TypeOfItem::face)][itemId(TypeOfItem::node)];
-  }
-
-  KOKKOS_INLINE_FUNCTION
-  ConnectivityMatrix nodeToCellMatrix() const
-  {
-    return m_item_to_item_matrix[itemId(TypeOfItem::node)][itemId(TypeOfItem::cell)];
-  }
-
   NodeValuePerCell<unsigned short> m_cell_to_node_local_cell;
 
   FaceValuePerCell<bool> m_cell_face_is_reversed;
@@ -266,19 +235,12 @@ class Connectivity
 
   CellValuePerNode<unsigned short> m_node_to_cell_local_node;
 
-  template <TypeOfItem item_type_0,
-            TypeOfItem item_type_1>
+  template <TypeOfItem item_type_0, TypeOfItem item_type_1>
   const ConnectivityMatrix& itemToItemMatrix() const
   {
-    const auto& item0_to_item1_matrix
-        = m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
-    return item0_to_item1_matrix;
+    return m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
   };
 
-  KOKKOS_INLINE_FUNCTION
-  const ConnectivityMatrix& itemToItemMatrix(const TypeOfItem& item_type_0,
-                                             const TypeOfItem& item_type_1) const;
-
 private:
   ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
 
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
index 50b788a4e360a87f29df46b0d6335816ab189f90..1150348d34f45c516a10e7ce82ec2d26c6520713 100644
--- a/src/mesh/ConnectivityUtils.hpp
+++ b/src/mesh/ConnectivityUtils.hpp
@@ -7,7 +7,7 @@
 template <TypeOfItem item_type_0,
           TypeOfItem item_type_1,
           typename ConnectivityType>
-const ConnectivityMatrix& getConnectivityMatrix(const ConnectivityType& connectivity)
+inline const ConnectivityMatrix& getConnectivityMatrix(const ConnectivityType& connectivity)
 {
   return connectivity.template itemToItemMatrix<item_type_0, item_type_1>();
 }
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 4b10a428e23b8357332abdceac72ff0ed2611426..9e46108737416efc4d594bda227b52565329c904 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -35,7 +35,10 @@ class MeshData
     if constexpr (dimension == 1) {
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
-      const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+      const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(m_mesh.connectivity());
+
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
           m_xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
@@ -44,7 +47,10 @@ class MeshData
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
       const Kokkos::View<const double*>& inv_cell_nb_nodes
           = m_mesh.connectivity().invCellNbNodes();
-      const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
+      const auto& cell_to_node_matrix
+          = getConnectivityMatrix<TypeOfItem::cell,
+                                  TypeOfItem::node>(m_mesh.connectivity());
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           Rd X = zero;
@@ -61,7 +67,9 @@ class MeshData
   void _updateVolume()
   {
     const Kokkos::View<const Rd*> xr = m_mesh.xr();
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(m_mesh.connectivity());
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
         double sum_cjr_xr = 0;
@@ -81,7 +89,9 @@ class MeshData
     }
     else if constexpr (dimension == 2) {
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
-      const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+      const auto& cell_to_node_matrix
+          = getConnectivityMatrix<TypeOfItem::cell,
+                                  TypeOfItem::node>(m_mesh.connectivity());
 
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -107,7 +117,10 @@ class MeshData
       const Kokkos::View<const Rd*> xr = m_mesh.xr();
 
       NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
-      const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+      const auto& face_to_node_matrix
+          = getConnectivityMatrix<TypeOfItem::face,
+                                  TypeOfItem::node>(m_mesh.connectivity());
+
       Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
           const auto& face_nodes = face_to_node_matrix.rowConst(l);
           const size_t nb_nodes = face_nodes.length;
@@ -132,8 +145,14 @@ class MeshData
           m_Cjr[jr] = zero;
         });
 
-      const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
-      const auto& cell_to_face_matrix = m_mesh.connectivity().cellToFaceMatrix();
+      const auto& cell_to_node_matrix
+          = getConnectivityMatrix<TypeOfItem::cell,
+                                  TypeOfItem::node>(m_mesh.connectivity());
+
+      const auto& cell_to_face_matrix
+          = getConnectivityMatrix<TypeOfItem::cell,
+                                  TypeOfItem::face>(m_mesh.connectivity());
+
       Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
           const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 56d3d254e314f7e18002f1f3b1624d03acc7218d..a1e3e9562df266fff75d0f25a6ed2b0bb90fdd99 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -8,6 +8,8 @@
 #include <RefNodeList.hpp>
 #include <RefFaceList.hpp>
 
+#include <ConnectivityMatrix.hpp>
+#include <ConnectivityUtils.hpp>
 #include <iostream>
 
 template <size_t dimension>
@@ -29,7 +31,10 @@ class MeshNodeBoundary
                    const RefFaceList& ref_face_list)
   {
     static_assert(dimension == MeshType::dimension);
-    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+    const auto& face_to_cell_matrix
+        = getConnectivityMatrix<TypeOfItem::face,
+                                TypeOfItem::cell>(mesh.connectivity());
+
     const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
     Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
         const auto& face_cells = face_to_cell_matrix.rowConst(face_list[l]);
@@ -42,7 +47,9 @@ class MeshNodeBoundary
     Kokkos::vector<unsigned int> node_ids;
     // not enough but should reduce significantly the number of resizing
     node_ids.reserve(dimension*face_list.extent(0));
-    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+    const auto& face_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::face,
+                                TypeOfItem::node>(mesh.connectivity());
 
     for (size_t l=0; l<face_list.extent(0); ++l) {
       const size_t face_number = face_list[l];
@@ -309,8 +316,14 @@ _getOutgoingNormal(const MeshType& mesh)
   const R normal = this->_getNormal(mesh);
 
   const Kokkos::View<const R*>& xr = mesh.xr();
-  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-  const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+  const auto& cell_to_node_matrix
+      = getConnectivityMatrix<TypeOfItem::cell,
+                              TypeOfItem::node>(mesh.connectivity());
+
+  const auto& node_to_cell_matrix
+          = getConnectivityMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>(mesh.connectivity());
+
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
   const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
@@ -340,8 +353,13 @@ _getOutgoingNormal(const MeshType& mesh)
   const R2 normal = this->_getNormal(mesh);
 
   const Kokkos::View<const R2*>& xr = mesh.xr();
-  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-  const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+  const auto& cell_to_node_matrix
+      = getConnectivityMatrix<TypeOfItem::cell,
+                              TypeOfItem::node>(mesh.connectivity());
+
+  const auto& node_to_cell_matrix
+          = getConnectivityMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>(mesh.connectivity());
 
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
@@ -372,8 +390,13 @@ _getOutgoingNormal(const MeshType& mesh)
   const R3 normal = this->_getNormal(mesh);
 
   const Kokkos::View<const R3*>& xr = mesh.xr();
-  const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
-  const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+  const auto& cell_to_node_matrix
+      = getConnectivityMatrix<TypeOfItem::cell,
+                              TypeOfItem::node>(mesh.connectivity());
+
+  const auto& node_to_cell_matrix
+          = getConnectivityMatrix<TypeOfItem::node,
+                                  TypeOfItem::cell>(mesh.connectivity());
 
   const size_t r0 = m_node_list[0];
   const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index d4ea1f33a91abfa67b5aa32863a8c4eb449e75e5..623576e67ee329d0a353062a9d323b8a98046a4f 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -6,6 +6,7 @@
 #include <iomanip>
 #include <sstream>
 #include <TinyVector.hpp>
+#include <ConnectivityUtils.hpp>
 
 class VTKWriter
 {
@@ -69,7 +70,11 @@ class VTKWriter
     fout << "<Cells>\n";
 
     fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
-    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
+
+    const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(mesh.connectivity());
+
     for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
       const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
       for (unsigned short r=0; r<cell_nodes.length; ++r) {
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 24953e541be234f8df618768ded573900c482387..39a43eddd6ac2a646d1007cc7527488c276b34ae 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -98,7 +98,9 @@ class AcousticSolver
   KOKKOS_INLINE_FUNCTION
   const Kokkos::View<const Rdd*>
   computeAr(const NodeValuePerCell<Rdd>& Ajr) {
-    const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
+    const auto& node_to_cell_matrix
+        = getConnectivityMatrix<TypeOfItem::node,
+                                TypeOfItem::cell>(m_connectivity);
 
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
         Rdd sum = zero;
@@ -121,7 +123,9 @@ class AcousticSolver
             const NodeValuePerCell<Rd>& Cjr,
             const Kokkos::View<const Rd*>& uj,
             const Kokkos::View<const double*>& pj) {
-    const auto& node_to_cell_matrix = m_connectivity.nodeToCellMatrix();
+    const auto& node_to_cell_matrix
+        = getConnectivityMatrix<TypeOfItem::node,
+                                TypeOfItem::cell>(m_connectivity);
 
     Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
         Rd& br = m_br(r);
@@ -200,7 +204,10 @@ class AcousticSolver
              const Kokkos::View<const Rd*>& uj,
              const Kokkos::View<const double*>& pj)
   {
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(m_mesh.connectivity());
+
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
 
@@ -276,7 +283,9 @@ class AcousticSolver
                      const Kokkos::View<const double*>& cj) const
   {
     const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(m_mesh.connectivity());
 
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
         const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
@@ -318,7 +327,9 @@ class AcousticSolver
 
     const NodeValuePerCell<Rd>& Fjr = m_Fjr;
     const Kokkos::View<const Rd*> ur = m_ur;
-    const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix
+        = getConnectivityMatrix<TypeOfItem::cell,
+                                TypeOfItem::node>(m_mesh.connectivity());
 
     const Kokkos::View<const double*> inv_mj = unknowns.invMj();
     Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {