diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
index f74ba52787244cb2afa38d7c5e09cbe5397c1733..4c094045402ae97c866117a1eee28d55e05cc47d 100644
--- a/src/mesh/MeshFlatNodeBoundary.cpp
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -15,8 +15,9 @@ MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension
 
   bool is_bad = false;
 
-  parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
-    const Rd& x = xr[this->m_node_list[r]];
+  auto node_list = this->m_ref_node_list.list();
+  parallel_for(node_list.size(), [=, &is_bad](int r) {
+    const Rd& x = xr[node_list[r]];
     if (std::abs(dot(x - origin, normal)) > 1E-13 * length) {
       is_bad = true;
     }
@@ -39,8 +40,9 @@ MeshFlatNodeBoundary<1>::_getNormal(const Mesh<Connectivity<1>>& mesh)
   const size_t number_of_bc_nodes = [&]() {
     size_t number_of_bc_nodes = 0;
     auto node_is_owned        = mesh.connectivity().nodeIsOwned();
-    for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
-      number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
+    auto node_list            = m_ref_node_list.list();
+    for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+      number_of_bc_nodes += (node_is_owned[node_list[i_node]]);
     }
     return parallel::allReduceMax(number_of_bc_nodes);
   }();
@@ -149,14 +151,14 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
   const R normal = this->_getNormal(mesh);
 
   double max_height = 0;
-
-  if (m_node_list.size() > 0) {
+  auto node_list    = m_ref_node_list.list();
+  if (node_list.size() > 0) {
     const NodeValue<const R>& xr    = mesh.xr();
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
     const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
-    const NodeId r0      = m_node_list[0];
+    const NodeId r0      = node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
@@ -193,13 +195,14 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
 
   double max_height = 0;
 
-  if (m_node_list.size() > 0) {
+  auto node_list = m_ref_node_list.list();
+  if (node_list.size() > 0) {
     const NodeValue<const R2>& xr   = mesh.xr();
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
     const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
-    const NodeId r0      = m_node_list[0];
+    const NodeId r0      = node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
     for (size_t r = 0; r < j0_nodes.size(); ++r) {
@@ -234,14 +237,14 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
   const R3 normal = this->_getNormal(mesh);
 
   double max_height = 0;
-
-  if (m_node_list.size() > 0) {
+  auto node_list    = m_ref_node_list.list();
+  if (node_list.size() > 0) {
     const NodeValue<const R3>& xr   = mesh.xr();
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
     const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
-    const NodeId r0      = m_node_list[0];
+    const NodeId r0      = node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
 
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
index 71ff3e65a96e122a01f8a9a246837c07be7e6a14..87425c616f9accab96db214f13b79b18013bf290 100644
--- a/src/mesh/MeshLineNodeBoundary.cpp
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -17,9 +17,10 @@ MeshLineNodeBoundary<Dimension>::_checkBoundaryIsLine(const TinyVector<Dimension
 
   const Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
 
-  bool is_bad = false;
-  parallel_for(this->m_node_list.size(), [=, &is_bad](int i_node) {
-    const Rd& x    = xr[this->m_node_list[i_node]];
+  bool is_bad    = false;
+  auto node_list = this->m_ref_node_list.list();
+  parallel_for(node_list.size(), [=, &is_bad](int i_node) {
+    const Rd& x    = xr[node_list[i_node]];
     const Rd delta = P * (x - origin);
     if (dot(delta, delta) > 1E-13 * length) {
       is_bad = true;
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
index b992fc77b9cb947545e5d60f8687ee987b3b4675..649936f9b8e51a4aba76a8a8e266ca0f87b434b7 100644
--- a/src/mesh/MeshNodeBoundary.cpp
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -32,8 +32,9 @@ MeshNodeBoundary<2>::_getBounds(const Mesh<Connectivity<2>>& mesh) const
     }
   };
 
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R2& x = xr[m_node_list[r]];
+  auto node_list = m_ref_node_list.list();
+  for (size_t r = 0; r < node_list.size(); ++r) {
+    const R2& x = xr[node_list[r]];
     update_xmin(x, xmin);
     update_xmax(x, xmax);
   }
@@ -127,8 +128,9 @@ MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
 
   const NodeValue<const R3>& xr = mesh.xr();
 
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R3& x = xr[m_node_list[r]];
+  auto node_list = m_ref_node_list.list();
+  for (size_t r = 0; r < node_list.size(); ++r) {
+    const R3& x = xr[node_list[r]];
     update_xmin(x, xmin);
     update_ymin(x, ymin);
     update_zmin(x, zmin);
@@ -202,18 +204,17 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
     Array<NodeId> node_list(node_ids.size());
     parallel_for(
       node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
-    m_node_list = node_list;
+    m_ref_node_list = RefNodeList{ref_face_list.refId(), node_list, ref_face_list.isBoundary()};
   } else {
     Array<NodeId> node_list(face_list.size());
     parallel_for(
       face_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<FaceId::base_type>(face_list[r]); });
-    m_node_list = node_list;
+    m_ref_node_list = RefNodeList{ref_face_list.refId(), node_list, ref_face_list.isBoundary()};
   }
 
   // This is quite dirty but it allows a non negligible performance
   // improvement
-  const_cast<Connectivity<Dimension>&>(mesh.connectivity())
-    .addRefItemList(RefItemList<ItemType::node>(ref_face_list.refId(), m_node_list, ref_face_list.isBoundary()));
+  const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
 template <size_t Dimension>
@@ -249,23 +250,22 @@ MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>
     Array<NodeId> node_list(node_ids.size());
     parallel_for(
       node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
-    m_node_list = node_list;
+    m_ref_node_list = RefNodeList{ref_edge_list.refId(), node_list, ref_edge_list.isBoundary()};
   } else {
     Array<NodeId> node_list(edge_list.size());
     parallel_for(
       edge_list.size(), PUGS_LAMBDA(int r) { node_list[r] = static_cast<EdgeId::base_type>(edge_list[r]); });
-    m_node_list = node_list;
+    m_ref_node_list = RefNodeList{ref_edge_list.refId(), node_list, ref_edge_list.isBoundary()};
   }
 
   // This is quite dirty but it allows a non negligible performance
   // improvement
-  const_cast<Connectivity<Dimension>&>(mesh.connectivity())
-    .addRefItemList(RefItemList<ItemType::node>(ref_edge_list.refId(), m_node_list, ref_edge_list.isBoundary()));
+  const_cast<Connectivity<Dimension>&>(mesh.connectivity()).addRefItemList(m_ref_node_list);
 }
 
 template <size_t Dimension>
 MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list)
-  : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
+  : m_ref_node_list(ref_node_list), m_boundary_name(ref_node_list.refId().tagName())
 {
   if (not ref_node_list.isBoundary()) {
     std::ostringstream ost;
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 0422118ffc14abbc85030d4d56ca2f2da29eadb9..719ec08129b287f4925f20078e9ec268f4d83075 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -16,7 +16,7 @@ template <size_t Dimension>
 class [[nodiscard]] MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
-  Array<const NodeId> m_node_list;
+  RefNodeList m_ref_node_list;
   std::string m_boundary_name;
 
   std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const Mesh<Connectivity<Dimension>>& mesh)
@@ -30,9 +30,16 @@ class [[nodiscard]] MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
   MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
+  PUGS_INLINE
+  const RefNodeList& refNodeList() const
+  {
+    return m_ref_node_list;
+  }
+
+  PUGS_INLINE
   const Array<const NodeId>& nodeList() const
   {
-    return m_node_list;
+    return m_ref_node_list.list();
   }
 
  protected: