diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 634ffe04fe88817042f028319ea1d22e8104d98c..81812defdaf2d71101a7c35fb36b466b3a20b33a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -7,6 +7,7 @@
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Connectivity.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
@@ -23,7 +24,6 @@
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/NamedBoundaryDescriptor.hpp>
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 6c0088831c86023e31914b953c9ea137e5411d21..97b03eda3487e6506c2d18f77a0ac5fd13d9b780 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -17,6 +17,10 @@ add_library(
   LogicalConnectivityBuilder.cpp
   MeshBuilderBase.cpp
   MeshDataManager.cpp
+  MeshFaceBoundary.cpp
+  MeshFlatNodeBoundary.cpp
+  MeshLineNodeBoundary.cpp
+  MeshNodeBoundary.cpp
   SynchronizerManager.cpp)
 
 # Additional dependencies
diff --git a/src/scheme/IBoundaryDescriptor.hpp b/src/mesh/IBoundaryDescriptor.hpp
similarity index 100%
rename from src/scheme/IBoundaryDescriptor.hpp
rename to src/mesh/IBoundaryDescriptor.hpp
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index c15bb1a66ddacc34fe392a08d133e6f683b0e757..b312e4a8e3296fceb7ca0b747bdf5841200b8274 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -199,85 +199,83 @@ class MeshData : public IMeshData
     const CellValue<const Rd> cell_iso_barycenter = this->cellIsoBarycenter();
     if constexpr (Dimension == 1) {
       m_cell_centroid = cell_iso_barycenter;
-    } else {
-      if constexpr (Dimension == 2) {
-        const CellValue<const double> Vj = this->Vj();
-        const NodeValue<const Rd>& xr    = m_mesh.xr();
-        const auto& cell_to_node_matrix  = m_mesh.connectivity().cellToNodeMatrix();
+    } else if constexpr (Dimension == 2) {
+      const CellValue<const double> Vj = this->Vj();
+      const NodeValue<const Rd>& xr    = m_mesh.xr();
+      const auto& cell_to_node_matrix  = m_mesh.connectivity().cellToNodeMatrix();
 
-        CellValue<Rd> cell_centroid{m_mesh.connectivity()};
-        parallel_for(
-          m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-            Rd sum = zero;
+      CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+          Rd sum = zero;
 
-            const auto& cell_nodes = cell_to_node_matrix[j];
+          const auto& cell_nodes = cell_to_node_matrix[j];
 
-            for (size_t R = 0; R < cell_nodes.size(); ++R) {
-              const Rd& xr0 = xr[cell_nodes[R]];
-              const Rd& xr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]];
+          for (size_t R = 0; R < cell_nodes.size(); ++R) {
+            const Rd& xr0 = xr[cell_nodes[R]];
+            const Rd& xr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]];
 
-              Rd xjxr0 = xr[cell_nodes[R]] - cell_iso_barycenter[j];
-              Rd xjxr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]] - cell_iso_barycenter[j];
+            Rd xjxr0 = xr[cell_nodes[R]] - cell_iso_barycenter[j];
+            Rd xjxr1 = xr[cell_nodes[(R + 1) % cell_nodes.size()]] - cell_iso_barycenter[j];
 
-              const double Vjl = 0.5 * (xjxr0[0] * xjxr1[1] - xjxr0[1] * xjxr1[0]);
+            const double Vjl = 0.5 * (xjxr0[0] * xjxr1[1] - xjxr0[1] * xjxr1[0]);
 
-              sum += Vjl * (xr0 + xr1 + cell_iso_barycenter[j]);
-            }
+            sum += Vjl * (xr0 + xr1 + cell_iso_barycenter[j]);
+          }
 
-            sum *= 1 / (3 * Vj[j]);
+          sum *= 1 / (3 * Vj[j]);
 
-            cell_centroid[j] = sum;
-          });
-        m_cell_centroid = cell_centroid;
-      } else {
-        const auto& face_center           = this->xl();
-        const CellValue<const double> Vj  = this->Vj();
-        const NodeValue<const Rd>& xr     = m_mesh.xr();
-        const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
-        const auto& face_to_node_matrix   = m_mesh.connectivity().faceToNodeMatrix();
-        const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
+          cell_centroid[j] = sum;
+        });
+      m_cell_centroid = cell_centroid;
+    } else {
+      const auto& face_center           = this->xl();
+      const CellValue<const double> Vj  = this->Vj();
+      const NodeValue<const Rd>& xr     = m_mesh.xr();
+      const auto& cell_to_face_matrix   = m_mesh.connectivity().cellToFaceMatrix();
+      const auto& face_to_node_matrix   = m_mesh.connectivity().faceToNodeMatrix();
+      const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
 
-        CellValue<Rd> cell_centroid{m_mesh.connectivity()};
-        parallel_for(
-          m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
-            const Rd xj = m_cell_iso_barycenter[j];
+      CellValue<Rd> cell_centroid{m_mesh.connectivity()};
+      parallel_for(
+        m_mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
+          const Rd xj = m_cell_iso_barycenter[j];
 
-            const auto& cell_faces = cell_to_face_matrix[j];
+          const auto& cell_faces = cell_to_face_matrix[j];
 
-            Rd sum = zero;
-            for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-              const FaceId face_id = cell_faces[i_face];
+          Rd sum = zero;
+          for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
+            const FaceId face_id = cell_faces[i_face];
 
-              const Rd xl = face_center[face_id];
+            const Rd xl = face_center[face_id];
 
-              const Rd xjxl = xl - xj;
+            const Rd xjxl = xl - xj;
 
-              const auto& face_nodes = face_to_node_matrix[face_id];
+            const auto& face_nodes = face_to_node_matrix[face_id];
 
-              const Rd xl_plus_xj = xl + xj;
+            const Rd xl_plus_xj = xl + xj;
 
-              double sign = (cell_face_is_reversed(j, i_face)) ? -1 : 1;
+            double sign = (cell_face_is_reversed(j, i_face)) ? -1 : 1;
 
-              for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
-                const Rd& xr0 = xr[face_nodes[i_face_node]];
-                const Rd& xr1 = xr[face_nodes[(i_face_node + 1) % face_nodes.size()]];
+            for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+              const Rd& xr0 = xr[face_nodes[i_face_node]];
+              const Rd& xr1 = xr[face_nodes[(i_face_node + 1) % face_nodes.size()]];
 
-                const Rd xjxr0 = xr0 - xj;
-                const Rd xjxr1 = xr1 - xj;
+              const Rd xjxr0 = xr0 - xj;
+              const Rd xjxr1 = xr1 - xj;
 
-                const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl);
+              const double Vjlr = dot(crossProduct(xjxr0, xjxr1), xjxl);
 
-                sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
-              }
+              sum += (sign * Vjlr) * (xl_plus_xj + xr0 + xr1);
             }
+          }
 
-            sum *= 1 / (24 * Vj[j]);
+          sum *= 1 / (24 * Vj[j]);
 
-            cell_centroid[j] = sum;
-          });
+          cell_centroid[j] = sum;
+        });
 
-        m_cell_centroid = cell_centroid;
-      }
+      m_cell_centroid = cell_centroid;
     }
   }
 
diff --git a/src/mesh/MeshFaceBoundary.cpp b/src/mesh/MeshFaceBoundary.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3690c5c0c614a844060fa4c266e75fdcea42e377
--- /dev/null
+++ b/src/mesh/MeshFaceBoundary.cpp
@@ -0,0 +1,36 @@
+#include <mesh/MeshFaceBoundary.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+MeshFaceBoundary<Dimension>::MeshFaceBoundary(const Mesh<Connectivity<Dimension>>&, const RefFaceList& ref_face_list)
+  : m_face_list(ref_face_list.list()), m_boundary_name(ref_face_list.refId().tagName())
+{}
+
+template MeshFaceBoundary<2>::MeshFaceBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshFaceBoundary<3>::MeshFaceBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+
+template <size_t Dimension>
+MeshFaceBoundary<Dimension>
+getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+  for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+       ++i_ref_face_list) {
+    const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+    const RefId& ref          = ref_face_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshFaceBoundary<Dimension>{mesh, ref_face_list};
+    }
+  }
+
+  std::ostringstream ost;
+  ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+  throw NormalError(ost.str());
+}
+
+template MeshFaceBoundary<1> getMeshFaceBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshFaceBoundary<2> getMeshFaceBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshFaceBoundary<3> getMeshFaceBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFaceBoundary.hpp b/src/mesh/MeshFaceBoundary.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..afd409d939b038009944a98ce57bfd38b89fda4d
--- /dev/null
+++ b/src/mesh/MeshFaceBoundary.hpp
@@ -0,0 +1,55 @@
+#ifndef MESH_FACE_BOUNDARY_HPP
+#define MESH_FACE_BOUNDARY_HPP
+
+#include <algebra/TinyVector.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <mesh/ItemValue.hpp>
+#include <mesh/RefItemList.hpp>
+#include <utils/Array.hpp>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+template <size_t Dimension>
+class MeshFaceBoundary   // clazy:exclude=copyable-polymorphic
+{
+ protected:
+  Array<const FaceId> m_face_list;
+  std::string m_boundary_name;
+
+  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(
+    const Mesh<Connectivity<Dimension>>& mesh) const;
+
+ public:
+  template <size_t MeshDimension>
+  friend MeshFaceBoundary<MeshDimension> getMeshFaceBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
+
+  MeshFaceBoundary& operator=(const MeshFaceBoundary&) = default;
+  MeshFaceBoundary& operator=(MeshFaceBoundary&&) = default;
+
+  const Array<const FaceId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+ protected:
+  MeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+
+ public:
+  MeshFaceBoundary(const MeshFaceBoundary&) = default;
+  MeshFaceBoundary(MeshFaceBoundary&&)      = default;
+
+  MeshFaceBoundary()          = default;
+  virtual ~MeshFaceBoundary() = default;
+};
+
+template <size_t Dimension>
+MeshFaceBoundary<Dimension> getMeshFaceBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                                const IBoundaryDescriptor& boundary_descriptor);
+
+#endif   // MESH_FACE_BOUNDARY_HPP
diff --git a/src/mesh/MeshFlatNodeBoundary.cpp b/src/mesh/MeshFlatNodeBoundary.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..059d6b004e1f3eaff9c9940dec3f2af7b8fbf57c
--- /dev/null
+++ b/src/mesh/MeshFlatNodeBoundary.cpp
@@ -0,0 +1,300 @@
+#include <mesh/MeshFlatNodeBoundary.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+void
+MeshFlatNodeBoundary<Dimension>::_checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+                                                      const TinyVector<Dimension, double>& origin,
+                                                      const double length,
+                                                      const Mesh<Connectivity<Dimension>>& mesh) const
+{
+  const NodeValue<const Rd>& xr = mesh.xr();
+
+  bool is_bad = false;
+
+  parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+    const Rd& x = xr[this->m_node_list[r]];
+    if (dot(x - origin, normal) > 1E-13 * length) {
+      is_bad = true;
+    }
+  });
+
+  if (parallel::allReduceOr(is_bad)) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": boundary is not flat!";
+    throw NormalError(ost.str());
+  }
+}
+
+template <>
+TinyVector<1, double>
+MeshFlatNodeBoundary<1>::_getNormal(const Mesh<Connectivity<1>>& mesh)
+{
+  using R = TinyVector<1, double>;
+
+  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]]);
+    }
+    return parallel::allReduceMax(number_of_bc_nodes);
+  }();
+
+  if (number_of_bc_nodes != 1) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
+        << ": node boundaries in 1D require to have exactly 1 node";
+    throw NormalError(ost.str());
+  }
+
+  return R{1};
+}
+
+template <>
+TinyVector<2, double>
+MeshFlatNodeBoundary<2>::_getNormal(const Mesh<Connectivity<2>>& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  if (xmin == xmax) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
+  }
+
+  R2 dx = xmax - xmin;
+  dx *= 1. / l2Norm(dx);
+
+  const R2 normal(-dx[1], dx[0]);
+
+  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
+
+  return normal;
+}
+
+template <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getNormal(const Mesh<Connectivity<3>>& mesh)
+{
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
+
+  const R3 u = xmax - xmin;
+  const R3 v = ymax - ymin;
+  const R3 w = zmax - zmin;
+
+  const R3 uv        = crossProduct(u, v);
+  const double uv_l2 = dot(uv, uv);
+
+  R3 normal        = uv;
+  double normal_l2 = uv_l2;
+
+  const R3 uw        = crossProduct(u, w);
+  const double uw_l2 = dot(uw, uw);
+
+  if (uw_l2 > uv_l2) {
+    normal    = uw;
+    normal_l2 = uw_l2;
+  }
+
+  const R3 vw        = crossProduct(v, w);
+  const double vw_l2 = dot(vw, vw);
+
+  if (vw_l2 > normal_l2) {
+    normal    = vw;
+    normal_l2 = vw_l2;
+  }
+
+  if (normal_l2 == 0) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
+  }
+
+  const double length = sqrt(normal_l2);
+
+  normal *= 1. / length;
+
+  this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
+
+  return normal;
+}
+
+template <>
+TinyVector<1, double>
+MeshFlatNodeBoundary<1>::_getOutgoingNormal(const Mesh<Connectivity<1>>& mesh)
+{
+  using R = TinyVector<1, double>;
+
+  const R normal = this->_getNormal(mesh);
+
+  double max_height = 0;
+
+  if (m_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 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) {
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
+      if (std::abs(height) > std::abs(max_height)) {
+        max_height = height;
+      }
+    }
+  }
+
+  Array<double> max_height_array = parallel::allGather(max_height);
+  for (size_t i = 0; i < max_height_array.size(); ++i) {
+    const double height = max_height_array[i];
+    if (std::abs(height) > std::abs(max_height)) {
+      max_height = height;
+    }
+  }
+
+  if (max_height > 0) {
+    return -normal;
+  } else {
+    return normal;
+  }
+}
+
+template <>
+TinyVector<2, double>
+MeshFlatNodeBoundary<2>::_getOutgoingNormal(const Mesh<Connectivity<2>>& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  const R2 normal = this->_getNormal(mesh);
+
+  double max_height = 0;
+
+  if (m_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 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) {
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
+      if (std::abs(height) > std::abs(max_height)) {
+        max_height = height;
+      }
+    }
+  }
+
+  Array<double> max_height_array = parallel::allGather(max_height);
+  for (size_t i = 0; i < max_height_array.size(); ++i) {
+    const double height = max_height_array[i];
+    if (std::abs(height) > std::abs(max_height)) {
+      max_height = height;
+    }
+  }
+
+  if (max_height > 0) {
+    return -normal;
+  } else {
+    return normal;
+  }
+}
+
+template <>
+TinyVector<3, double>
+MeshFlatNodeBoundary<3>::_getOutgoingNormal(const Mesh<Connectivity<3>>& mesh)
+{
+  using R3 = TinyVector<3, double>;
+
+  const R3 normal = this->_getNormal(mesh);
+
+  double max_height = 0;
+
+  if (m_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 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) {
+      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
+      if (std::abs(height) > std::abs(max_height)) {
+        max_height = height;
+      }
+    }
+  }
+
+  Array<double> max_height_array = parallel::allGather(max_height);
+  for (size_t i = 0; i < max_height_array.size(); ++i) {
+    const double height = max_height_array[i];
+    if (std::abs(height) > std::abs(max_height)) {
+      max_height = height;
+    }
+  }
+
+  if (max_height > 0) {
+    return -normal;
+  } else {
+    return normal;
+  }
+}
+
+template <size_t Dimension>
+MeshFlatNodeBoundary<Dimension>
+getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+  for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+       ++i_ref_node_list) {
+    const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+    const RefId& ref          = ref_node_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshFlatNodeBoundary<Dimension>{mesh, ref_node_list};
+    }
+  }
+  for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+       ++i_ref_face_list) {
+    const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+    const RefId& ref          = ref_face_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshFlatNodeBoundary<Dimension>{mesh, ref_face_list};
+    }
+  }
+
+  std::ostringstream ost;
+  ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+  throw NormalError(ost.str());
+}
+
+template MeshFlatNodeBoundary<1> getMeshFlatNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<2> getMeshFlatNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshFlatNodeBoundary<3> getMeshFlatNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshFlatNodeBoundary.hpp b/src/mesh/MeshFlatNodeBoundary.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..334cf03aac844670c084b8385ebbb00d1cc9507d
--- /dev/null
+++ b/src/mesh/MeshFlatNodeBoundary.hpp
@@ -0,0 +1,60 @@
+#ifndef MESH_FLAT_NODE_BOUNDARY_HPP
+#define MESH_FLAT_NODE_BOUNDARY_HPP
+
+#include <mesh/MeshNodeBoundary.hpp>
+
+template <size_t Dimension>
+class MeshFlatNodeBoundary final : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Rd m_outgoing_normal;
+
+  Rd _getNormal(const Mesh<Connectivity<Dimension>>& mesh);
+
+  void _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+                            const TinyVector<Dimension, double>& origin,
+                            const double length,
+                            const Mesh<Connectivity<Dimension>>& mesh) const;
+
+  Rd _getOutgoingNormal(const Mesh<Connectivity<Dimension>>& mesh);
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_outgoing_normal;
+  }
+
+  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
+
+  template <size_t MeshDimension>
+  friend MeshFlatNodeBoundary<MeshDimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+                                                                     const IBoundaryDescriptor& boundary_descriptor);
+
+ private:
+  template <typename MeshType>
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+  {}
+
+  template <typename MeshType>
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
+  {}
+
+ public:
+  MeshFlatNodeBoundary()                            = default;
+  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
+  ~MeshFlatNodeBoundary()                           = default;
+};
+
+template <size_t Dimension>
+MeshFlatNodeBoundary<Dimension> getMeshFlatNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                                        const IBoundaryDescriptor& boundary_descriptor);
+
+#endif   // MESH_FLAT_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshLineNodeBoundary.cpp b/src/mesh/MeshLineNodeBoundary.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..cae862f33dd9a60d5b0e895aba4907cf5a5b102a
--- /dev/null
+++ b/src/mesh/MeshLineNodeBoundary.cpp
@@ -0,0 +1,147 @@
+#include <mesh/MeshLineNodeBoundary.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+void
+MeshLineNodeBoundary<Dimension>::_checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+                                                      const TinyVector<Dimension, double>& origin,
+                                                      const double length,
+                                                      const Mesh<Connectivity<Dimension>>& mesh) const
+{
+  using Rdxd = TinyMatrix<Dimension>;
+
+  const NodeValue<const Rd>& xr = mesh.xr();
+
+  Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
+
+  bool is_bad = false;
+  parallel_for(this->m_node_list.size(), [=, &is_bad](int node_id) {
+    const Rd& x    = xr[this->m_node_list[node_id]];
+    const Rd delta = x - origin;
+    if (dot(P * delta, direction) > 1E-13 * length) {
+      is_bad = true;
+    }
+  });
+
+  if (parallel::allReduceOr(is_bad)) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": boundary is not a line!";
+    throw NormalError(ost.str());
+  }
+}
+
+template <>
+template <>
+TinyVector<1>
+MeshLineNodeBoundary<1>::_getDirection(const Mesh<Connectivity<1>>&)
+{
+  throw UnexpectedError("MeshLineNodeBoundary makes no sense in dimension 1");
+}
+
+template <>
+template <>
+TinyVector<2>
+MeshLineNodeBoundary<2>::_getDirection(const Mesh<Connectivity<2>>& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  if (xmin == xmax) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute direction";
+    throw NormalError(ost.str());
+  }
+
+  R2 direction        = xmax - xmin;
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
+template <>
+template <>
+TinyVector<3>
+MeshLineNodeBoundary<3>::_getDirection(const Mesh<Connectivity<3>>& mesh)
+{
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
+
+  const R3 u = xmax - xmin;
+  const R3 v = ymax - ymin;
+  const R3 w = zmax - zmin;
+
+  R3 direction = u;
+  if (dot(v, v) > dot(direction, direction)) {
+    direction = v;
+  }
+
+  if (dot(w, w) > dot(direction, direction)) {
+    direction = w;
+  }
+
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
+template <size_t Dimension>
+MeshLineNodeBoundary<Dimension>
+getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+  for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+       ++i_ref_node_list) {
+    const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+    const RefId& ref          = ref_node_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshLineNodeBoundary<Dimension>{mesh, ref_node_list};
+    }
+  }
+  for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
+       ++i_ref_edge_list) {
+    const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
+    const RefId& ref          = ref_edge_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshLineNodeBoundary<Dimension>{mesh, ref_edge_list};
+    }
+  }
+  for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+       ++i_ref_face_list) {
+    const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+    const RefId& ref          = ref_face_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshLineNodeBoundary<Dimension>{mesh, ref_face_list};
+    }
+  }
+
+  std::ostringstream ost;
+  ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+  throw NormalError(ost.str());
+}
+
+template MeshLineNodeBoundary<1> getMeshLineNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<2> getMeshLineNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshLineNodeBoundary<3> getMeshLineNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshLineNodeBoundary.hpp b/src/mesh/MeshLineNodeBoundary.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0ae35d058d3b749453463bd347edb3f69c2a109e
--- /dev/null
+++ b/src/mesh/MeshLineNodeBoundary.hpp
@@ -0,0 +1,65 @@
+#ifndef MESH_LINE_NODE_BOUNDARY_HPP
+#define MESH_LINE_NODE_BOUNDARY_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+
+template <size_t Dimension>
+class MeshLineNodeBoundary final : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Rd m_direction;
+
+  template <size_t MeshDimension>
+  TinyVector<MeshDimension> _getDirection(const Mesh<Connectivity<MeshDimension>>&);
+
+  void _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+                            const TinyVector<Dimension, double>& origin,
+                            const double length,
+                            const Mesh<Connectivity<Dimension>>& mesh) const;
+
+ public:
+  template <size_t MeshDimension>
+  friend MeshLineNodeBoundary<MeshDimension> getMeshLineNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+                                                                     const IBoundaryDescriptor& boundary_descriptor);
+
+  PUGS_INLINE
+  const Rd&
+  direction() const
+  {
+    return m_direction;
+  }
+
+  MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
+
+ private:
+  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
+  {}
+
+  MeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
+  {}
+
+ public:
+  MeshLineNodeBoundary()                            = default;
+  MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary(MeshLineNodeBoundary&&)      = default;
+  ~MeshLineNodeBoundary()                           = default;
+};
+
+template <size_t Dimension>
+MeshLineNodeBoundary<Dimension> getMeshLineNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                                        const IBoundaryDescriptor& boundary_descriptor);
+
+#endif   // MESH_LINE_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshNodeBoundary.cpp b/src/mesh/MeshNodeBoundary.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..6a54afb554c1d1454bd01934d2dd08f4c023d62f
--- /dev/null
+++ b/src/mesh/MeshNodeBoundary.cpp
@@ -0,0 +1,328 @@
+#include <mesh/MeshNodeBoundary.hpp>
+
+#include <Kokkos_Vector.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/Mesh.hpp>
+#include <utils/Messenger.hpp>
+
+template <>
+std::array<TinyVector<2>, 2>
+MeshNodeBoundary<2>::_getBounds(const Mesh<Connectivity<2>>& mesh) const
+{
+  using R2 = TinyVector<2, double>;
+
+  const NodeValue<const R2>& xr = mesh.xr();
+
+  std::array<R2, 2> bounds;
+  R2& xmin = bounds[0];
+  R2& xmax = bounds[1];
+
+  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
+
+  auto update_xmin = [](const R2& x, R2& xmin) {
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R2& x, R2& xmax) {
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
+      xmax = x;
+    }
+  };
+
+  for (size_t r = 0; r < m_node_list.size(); ++r) {
+    const R2& x = xr[m_node_list[r]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+  }
+
+  if (parallel::size() > 1) {
+    Array<R2> xmin_array = parallel::allGather(xmin);
+    Array<R2> xmax_array = parallel::allGather(xmax);
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+  }
+
+  return bounds;
+}
+
+template <>
+std::array<TinyVector<3>, 6>
+MeshNodeBoundary<3>::_getBounds(const Mesh<Connectivity<3>>& mesh) const
+{
+  using R3 = TinyVector<3, double>;
+
+  auto update_xmin = [](const R3& x, R3& xmin) {
+    // XMIN: X.xmin X.ymax X.zmax
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
+        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R3& x, R3& xmax) {
+    // XMAX: X.xmax X.ymin X.zmin
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
+        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
+      xmax = x;
+    }
+  };
+
+  auto update_ymin = [](const R3& x, R3& ymin) {
+    // YMIN: X.ymin X.zmax X.xmin
+    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
+        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
+      ymin = x;
+    }
+  };
+
+  auto update_ymax = [](const R3& x, R3& ymax) {
+    // YMAX: X.ymax X.zmin X.xmax
+    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
+        ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
+      ymax = x;
+    }
+  };
+
+  auto update_zmin = [](const R3& x, R3& zmin) {
+    // ZMIN: X.zmin X.xmin X.ymin
+    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
+        ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
+      zmin = x;
+    }
+  };
+
+  auto update_zmax = [](const R3& x, R3& zmax) {
+    // ZMAX: X.zmax X.xmax X.ymax
+    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
+        ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
+      zmax = x;
+    }
+  };
+
+  std::array<R3, 6> bounds;
+  R3& xmin = bounds[0];
+  R3& ymin = bounds[1];
+  R3& zmin = bounds[2];
+  R3& xmax = bounds[3];
+  R3& ymax = bounds[4];
+  R3& zmax = bounds[5];
+
+  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymin = xmin;
+  zmin = xmin;
+
+  xmax = -xmin;
+  ymax = xmax;
+  zmax = xmax;
+
+  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]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+    update_ymin(x, ymin);
+    update_ymax(x, ymax);
+    update_zmin(x, zmin);
+    update_zmax(x, zmax);
+  }
+
+  if (parallel::size() > 0) {
+    Array<const R3> xmin_array = parallel::allGather(xmin);
+    Array<const R3> xmax_array = parallel::allGather(xmax);
+    Array<const R3> ymin_array = parallel::allGather(ymin);
+    Array<const R3> ymax_array = parallel::allGather(ymax);
+    Array<const R3> zmin_array = parallel::allGather(zmin);
+    Array<const R3> zmax_array = parallel::allGather(zmax);
+
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < ymin_array.size(); ++i) {
+      update_ymin(ymin_array[i], ymin);
+    }
+    for (size_t i = 0; i < zmin_array.size(); ++i) {
+      update_zmin(zmin_array[i], zmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+    for (size_t i = 0; i < ymax_array.size(); ++i) {
+      update_ymax(ymax_array[i], ymax);
+    }
+    for (size_t i = 0; i < zmax_array.size(); ++i) {
+      update_zmax(zmax_array[i], zmax);
+    }
+  }
+
+  return bounds;
+}
+
+template <size_t Dimension>
+MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                              const RefFaceList& ref_face_list)
+  : m_boundary_name(ref_face_list.refId().tagName())
+{
+  const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+  const Array<const FaceId>& face_list = ref_face_list.list();
+
+  bool is_bad = false;
+  parallel_for(face_list.size(), [=, &is_bad](int l) {
+    const auto& face_cells = face_to_cell_matrix[face_list[l]];
+    if (face_cells.size() > 1) {
+      is_bad = true;
+    }
+  });
+
+  if (parallel::allReduceOr(is_bad)) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": inner faces cannot be used to define mesh boundaries";
+    throw NormalError(ost.str());
+  }
+
+  Kokkos::vector<unsigned int> node_ids;
+  // not enough but should reduce significantly the number of resizing
+  node_ids.reserve(Dimension * face_list.size());
+  const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
+
+  for (size_t l = 0; l < face_list.size(); ++l) {
+    const FaceId face_number = face_list[l];
+    const auto& face_nodes   = face_to_node_matrix[face_number];
+
+    for (size_t r = 0; r < face_nodes.size(); ++r) {
+      node_ids.push_back(face_nodes[r]);
+    }
+  }
+  std::sort(node_ids.begin(), node_ids.end());
+  auto last = std::unique(node_ids.begin(), node_ids.end());
+  node_ids.resize(std::distance(node_ids.begin(), last));
+
+  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;
+}
+
+template <size_t Dimension>
+MeshNodeBoundary<Dimension>::MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                              const RefEdgeList& ref_edge_list)
+  : m_boundary_name(ref_edge_list.refId().tagName())
+{
+  const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
+  const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+  auto edge_is_owned = mesh.connectivity().edgeIsOwned();
+
+  const Array<const EdgeId>& edge_list = ref_edge_list.list();
+
+  bool is_bad = false;
+  parallel_for(edge_list.size(), [=, &is_bad](int l) {
+    const EdgeId edge_id = edge_list[l];
+    if (edge_is_owned[edge_id]) {
+      const auto& edge_faces             = edge_to_face_matrix[edge_id];
+      bool is_connected_to_boundary_face = false;
+      for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) {
+        const FaceId edge_face_id = edge_faces[i_edge_face];
+        if (face_to_cell_matrix[edge_face_id].size() == 1) {
+          is_connected_to_boundary_face = true;
+        }
+      }
+      if (not is_connected_to_boundary_face) {
+        is_bad = true;
+      }
+    }
+  });
+
+  if (parallel::allReduceOr(is_bad)) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": inner edges cannot be used to define mesh boundaries";
+    throw NormalError(ost.str());
+  }
+
+  Kokkos::vector<unsigned int> node_ids;
+  node_ids.reserve(2 * edge_list.size());
+  const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
+
+  for (size_t l = 0; l < edge_list.size(); ++l) {
+    const EdgeId edge_number = edge_list[l];
+    const auto& edge_nodes   = edge_to_node_matrix[edge_number];
+
+    for (size_t r = 0; r < edge_nodes.size(); ++r) {
+      node_ids.push_back(edge_nodes[r]);
+    }
+  }
+  std::sort(node_ids.begin(), node_ids.end());
+  auto last = std::unique(node_ids.begin(), node_ids.end());
+  node_ids.resize(std::distance(node_ids.begin(), last));
+
+  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;
+}
+
+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())
+{}
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefFaceList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefFaceList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefFaceList&);
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefEdgeList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefEdgeList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefEdgeList&);
+
+template MeshNodeBoundary<1>::MeshNodeBoundary(const Mesh<Connectivity<1>>&, const RefNodeList&);
+template MeshNodeBoundary<2>::MeshNodeBoundary(const Mesh<Connectivity<2>>&, const RefNodeList&);
+template MeshNodeBoundary<3>::MeshNodeBoundary(const Mesh<Connectivity<3>>&, const RefNodeList&);
+
+template <size_t Dimension>
+MeshNodeBoundary<Dimension>
+getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const IBoundaryDescriptor& boundary_descriptor)
+{
+  for (size_t i_ref_node_list = 0; i_ref_node_list < mesh.connectivity().template numberOfRefItemList<ItemType::node>();
+       ++i_ref_node_list) {
+    const auto& ref_node_list = mesh.connectivity().template refItemList<ItemType::node>(i_ref_node_list);
+    const RefId& ref          = ref_node_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshNodeBoundary<Dimension>{mesh, ref_node_list};
+    }
+  }
+  for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<ItemType::edge>();
+       ++i_ref_edge_list) {
+    const auto& ref_edge_list = mesh.connectivity().template refItemList<ItemType::edge>(i_ref_edge_list);
+    const RefId& ref          = ref_edge_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshNodeBoundary<Dimension>{mesh, ref_edge_list};
+    }
+  }
+  for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<ItemType::face>();
+       ++i_ref_face_list) {
+    const auto& ref_face_list = mesh.connectivity().template refItemList<ItemType::face>(i_ref_face_list);
+    const RefId& ref          = ref_face_list.refId();
+    if (ref == boundary_descriptor) {
+      return MeshNodeBoundary<Dimension>{mesh, ref_face_list};
+    }
+  }
+
+  std::ostringstream ost;
+  ost << "cannot find surface with name " << rang::fgB::red << boundary_descriptor << rang::style::reset;
+
+  throw NormalError(ost.str());
+}
+
+template MeshNodeBoundary<1> getMeshNodeBoundary(const Mesh<Connectivity<1>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary<2> getMeshNodeBoundary(const Mesh<Connectivity<2>>&, const IBoundaryDescriptor&);
+template MeshNodeBoundary<3> getMeshNodeBoundary(const Mesh<Connectivity<3>>&, const IBoundaryDescriptor&);
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 0743d2d95dda05ca0cb990ade1c8d088c9ed2112..883a9d975682b55e0cad75c191e7f2be14e59190 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -1,16 +1,17 @@
 #ifndef MESH_NODE_BOUNDARY_HPP
 #define MESH_NODE_BOUNDARY_HPP
 
-#include <Kokkos_Vector.hpp>
-#include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
-#include <mesh/ConnectivityMatrix.hpp>
-#include <mesh/IConnectivity.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <mesh/ItemValue.hpp>
 #include <mesh/RefItemList.hpp>
 #include <utils/Array.hpp>
-#include <utils/Exceptions.hpp>
-#include <utils/Messenger.hpp>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
 
 template <size_t Dimension>
 class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
@@ -19,10 +20,14 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   Array<const NodeId> m_node_list;
   std::string m_boundary_name;
 
-  template <typename MeshType>
-  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const MeshType& mesh) const;
+  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(
+    const Mesh<Connectivity<Dimension>>& mesh) const;
 
  public:
+  template <size_t MeshDimension>
+  friend MeshNodeBoundary<MeshDimension> getMeshNodeBoundary(const Mesh<Connectivity<MeshDimension>>& mesh,
+                                                             const IBoundaryDescriptor& boundary_descriptor);
+
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
   MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
@@ -32,120 +37,12 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
     return m_node_list;
   }
 
-  template <typename MeshType>
-  MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : m_boundary_name(ref_face_list.refId().tagName())
-  {
-    static_assert(Dimension == MeshType::Dimension);
-    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
-
-    const Array<const FaceId>& face_list = ref_face_list.list();
-
-    bool is_bad = false;
-    parallel_for(face_list.size(), [=, &is_bad](int l) {
-      const auto& face_cells = face_to_cell_matrix[face_list[l]];
-      if (face_cells.size() > 1) {
-        is_bad = true;
-      }
-    });
-
-    if (parallel::allReduceOr(is_bad)) {
-      std::ostringstream ost;
-      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-          << ": inner faces cannot be used to define mesh boundaries";
-      throw NormalError(ost.str());
-    }
-
-    Kokkos::vector<unsigned int> node_ids;
-    // not enough but should reduce significantly the number of resizing
-    node_ids.reserve(Dimension * face_list.size());
-    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
-
-    for (size_t l = 0; l < face_list.size(); ++l) {
-      const FaceId face_number = face_list[l];
-      const auto& face_nodes   = face_to_node_matrix[face_number];
-
-      for (size_t r = 0; r < face_nodes.size(); ++r) {
-        node_ids.push_back(face_nodes[r]);
-      }
-    }
-    std::sort(node_ids.begin(), node_ids.end());
-    auto last = std::unique(node_ids.begin(), node_ids.end());
-    node_ids.resize(std::distance(node_ids.begin(), last));
-
-    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;
-  }
-
-  template <typename MeshType>
-  MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
-    : m_boundary_name(ref_edge_list.refId().tagName())
-  {
-    static_assert(Dimension == MeshType::Dimension);
-    const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
-    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
-
-    auto edge_is_owned = mesh.connectivity().edgeIsOwned();
-
-    const Array<const EdgeId>& edge_list = ref_edge_list.list();
-
-    bool is_bad = false;
-    parallel_for(edge_list.size(), [=, &is_bad](int l) {
-      const EdgeId edge_id = edge_list[l];
-      if (edge_is_owned[edge_id]) {
-        const auto& edge_faces             = edge_to_face_matrix[edge_id];
-        bool is_connected_to_boundary_face = false;
-        for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) {
-          const FaceId edge_face_id = edge_faces[i_edge_face];
-          if (face_to_cell_matrix[edge_face_id].size() == 1) {
-            is_connected_to_boundary_face = true;
-          }
-        }
-        if (not is_connected_to_boundary_face) {
-          is_bad = true;
-        }
-      }
-    });
-
-    if (parallel::allReduceOr(is_bad)) {
-      std::ostringstream ost;
-      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-          << ": inner edges cannot be used to define mesh boundaries";
-      throw NormalError(ost.str());
-    }
-
-    Kokkos::vector<unsigned int> node_ids;
-    node_ids.reserve(2 * edge_list.size());
-    const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
-
-    for (size_t l = 0; l < edge_list.size(); ++l) {
-      const EdgeId edge_number = edge_list[l];
-      const auto& edge_nodes   = edge_to_node_matrix[edge_number];
-
-      for (size_t r = 0; r < edge_nodes.size(); ++r) {
-        node_ids.push_back(edge_nodes[r]);
-      }
-    }
-    std::sort(node_ids.begin(), node_ids.end());
-    auto last = std::unique(node_ids.begin(), node_ids.end());
-    node_ids.resize(std::distance(node_ids.begin(), last));
-
-    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;
-  }
-
-  template <typename MeshType>
-  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
-    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
-
-  {
-    static_assert(Dimension == MeshType::Dimension);
-  }
+ protected:
+  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefFaceList& ref_face_list);
+  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh, const RefEdgeList& ref_edge_list);
+  MeshNodeBoundary(const Mesh<Connectivity<Dimension>>&, const RefNodeList& ref_node_list);
 
+ public:
   MeshNodeBoundary(const MeshNodeBoundary&) = default;
   MeshNodeBoundary(MeshNodeBoundary&&)      = default;
 
@@ -153,649 +50,8 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   virtual ~MeshNodeBoundary() = default;
 };
 
-template <>
-template <typename MeshType>
-std::array<TinyVector<2>, 2>
-MeshNodeBoundary<2>::_getBounds(const MeshType& mesh) const
-{
-  static_assert(MeshType::Dimension == 2);
-
-  using R2 = TinyVector<2, double>;
-
-  const NodeValue<const R2>& xr = mesh.xr();
-
-  std::array<R2, 2> bounds;
-  R2& xmin = bounds[0];
-  R2& xmax = bounds[1];
-
-  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
-
-  auto update_xmin = [](const R2& x, R2& xmin) {
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
-      xmin = x;
-    }
-  };
-
-  auto update_xmax = [](const R2& x, R2& xmax) {
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
-      xmax = x;
-    }
-  };
-
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R2& x = xr[m_node_list[r]];
-    update_xmin(x, xmin);
-    update_xmax(x, xmax);
-  }
-
-  if (parallel::size() > 1) {
-    Array<R2> xmin_array = parallel::allGather(xmin);
-    Array<R2> xmax_array = parallel::allGather(xmax);
-    for (size_t i = 0; i < xmin_array.size(); ++i) {
-      update_xmin(xmin_array[i], xmin);
-    }
-    for (size_t i = 0; i < xmax_array.size(); ++i) {
-      update_xmax(xmax_array[i], xmax);
-    }
-  }
-
-  return bounds;
-}
-
-template <>
-template <typename MeshType>
-std::array<TinyVector<3>, 6>
-MeshNodeBoundary<3>::_getBounds(const MeshType& mesh) const
-{
-  static_assert(MeshType::Dimension == 3);
-
-  using R3 = TinyVector<3, double>;
-
-  auto update_xmin = [](const R3& x, R3& xmin) {
-    // XMIN: X.xmin X.ymax X.zmax
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
-        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
-      xmin = x;
-    }
-  };
-
-  auto update_xmax = [](const R3& x, R3& xmax) {
-    // XMAX: X.xmax X.ymin X.zmin
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
-        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
-      xmax = x;
-    }
-  };
-
-  auto update_ymin = [](const R3& x, R3& ymin) {
-    // YMIN: X.ymin X.zmax X.xmin
-    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
-        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
-      ymin = x;
-    }
-  };
-
-  auto update_ymax = [](const R3& x, R3& ymax) {
-    // YMAX: X.ymax X.zmin X.xmax
-    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
-        ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
-      ymax = x;
-    }
-  };
-
-  auto update_zmin = [](const R3& x, R3& zmin) {
-    // ZMIN: X.zmin X.xmin X.ymin
-    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
-        ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
-      zmin = x;
-    }
-  };
-
-  auto update_zmax = [](const R3& x, R3& zmax) {
-    // ZMAX: X.zmax X.xmax X.ymax
-    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
-        ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
-      zmax = x;
-    }
-  };
-
-  std::array<R3, 6> bounds;
-  R3& xmin = bounds[0];
-  R3& ymin = bounds[1];
-  R3& zmin = bounds[2];
-  R3& xmax = bounds[3];
-  R3& ymax = bounds[4];
-  R3& zmax = bounds[5];
-
-  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
-  ymin = xmin;
-  zmin = xmin;
-
-  xmax = -xmin;
-  ymax = xmax;
-  zmax = xmax;
-
-  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]];
-    update_xmin(x, xmin);
-    update_xmax(x, xmax);
-    update_ymin(x, ymin);
-    update_ymax(x, ymax);
-    update_zmin(x, zmin);
-    update_zmax(x, zmax);
-  }
-
-  if (parallel::size() > 0) {
-    Array<const R3> xmin_array = parallel::allGather(xmin);
-    Array<const R3> xmax_array = parallel::allGather(xmax);
-    Array<const R3> ymin_array = parallel::allGather(ymin);
-    Array<const R3> ymax_array = parallel::allGather(ymax);
-    Array<const R3> zmin_array = parallel::allGather(zmin);
-    Array<const R3> zmax_array = parallel::allGather(zmax);
-
-    for (size_t i = 0; i < xmin_array.size(); ++i) {
-      update_xmin(xmin_array[i], xmin);
-    }
-    for (size_t i = 0; i < ymin_array.size(); ++i) {
-      update_ymin(ymin_array[i], ymin);
-    }
-    for (size_t i = 0; i < zmin_array.size(); ++i) {
-      update_zmin(zmin_array[i], zmin);
-    }
-    for (size_t i = 0; i < xmax_array.size(); ++i) {
-      update_xmax(xmax_array[i], xmax);
-    }
-    for (size_t i = 0; i < ymax_array.size(); ++i) {
-      update_ymax(ymax_array[i], ymax);
-    }
-    for (size_t i = 0; i < zmax_array.size(); ++i) {
-      update_zmax(zmax_array[i], zmax);
-    }
-  }
-
-  return bounds;
-}
-
-template <size_t Dimension>
-class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
-{
- public:
-  using Rd = TinyVector<Dimension, double>;
-
- private:
-  const Rd m_outgoing_normal;
-
-  template <typename MeshType>
-  PUGS_INLINE Rd _getNormal(const MeshType& mesh);
-
-  template <typename MeshType>
-  PUGS_INLINE void
-  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
-                       const TinyVector<Dimension, double>& origin,
-                       const double length,
-                       const MeshType& mesh) const
-  {
-    static_assert(MeshType::Dimension == Dimension);
-
-    const NodeValue<const Rd>& xr = mesh.xr();
-
-    bool is_bad = false;
-
-    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
-      const Rd& x = xr[this->m_node_list[r]];
-      if (dot(x - origin, normal) > 1E-13 * length) {
-        is_bad = true;
-      }
-    });
-
-    if (parallel::allReduceOr(is_bad)) {
-      std::ostringstream ost;
-      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-          << ": boundary is not flat!";
-      throw NormalError(ost.str());
-    }
-  }
-
-  template <typename MeshType>
-  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
-
- public:
-  const Rd&
-  outgoingNormal() const
-  {
-    return m_outgoing_normal;
-  }
-
-  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
-
-  template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
-  {
-    ;
-  }
-
-  template <typename MeshType>
-  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
-  {
-    ;
-  }
-
-  MeshFlatNodeBoundary()                            = default;
-  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
-  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&)      = default;
-  virtual ~MeshFlatNodeBoundary()                   = default;
-};
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 1);
-  using R = TinyVector<1, double>;
-
-  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]]);
-    }
-    return parallel::allReduceMax(number_of_bc_nodes);
-  }();
-
-  if (number_of_bc_nodes != 1) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
-        << ": node boundaries in 1D require to have exactly 1 node";
-    throw NormalError(ost.str());
-  }
-
-  return R{1};
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 2);
-  using R2 = TinyVector<2, double>;
-
-  std::array<R2, 2> bounds = this->_getBounds(mesh);
-
-  const R2& xmin = bounds[0];
-  const R2& xmax = bounds[1];
-
-  if (xmin == xmax) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-        << ": unable to compute normal";
-    throw NormalError(ost.str());
-  }
-
-  R2 dx = xmax - xmin;
-  dx *= 1. / l2Norm(dx);
-
-  const R2 normal(-dx[1], dx[0]);
-
-  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
-
-  return normal;
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 3);
-
-  using R3 = TinyVector<3, double>;
-
-  std::array<R3, 6> bounds = this->_getBounds(mesh);
-
-  const R3& xmin = bounds[0];
-  const R3& ymin = bounds[1];
-  const R3& zmin = bounds[2];
-  const R3& xmax = bounds[3];
-  const R3& ymax = bounds[4];
-  const R3& zmax = bounds[5];
-
-  const R3 u = xmax - xmin;
-  const R3 v = ymax - ymin;
-  const R3 w = zmax - zmin;
-
-  const R3 uv        = crossProduct(u, v);
-  const double uv_l2 = dot(uv, uv);
-
-  R3 normal        = uv;
-  double normal_l2 = uv_l2;
-
-  const R3 uw        = crossProduct(u, w);
-  const double uw_l2 = dot(uw, uw);
-
-  if (uw_l2 > uv_l2) {
-    normal    = uw;
-    normal_l2 = uw_l2;
-  }
-
-  const R3 vw        = crossProduct(v, w);
-  const double vw_l2 = dot(vw, vw);
-
-  if (vw_l2 > normal_l2) {
-    normal    = vw;
-    normal_l2 = vw_l2;
-  }
-
-  if (normal_l2 == 0) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-        << ": unable to compute normal";
-    throw NormalError(ost.str());
-  }
-
-  const double length = sqrt(normal_l2);
-
-  normal *= 1. / length;
-
-  this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
-
-  return normal;
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 1);
-  using R = TinyVector<1, double>;
-
-  const R normal = this->_getNormal(mesh);
-
-  double max_height = 0;
-
-  if (m_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 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) {
-      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
-      if (std::abs(height) > std::abs(max_height)) {
-        max_height = height;
-      }
-    }
-  }
-
-  Array<double> max_height_array = parallel::allGather(max_height);
-  for (size_t i = 0; i < max_height_array.size(); ++i) {
-    const double height = max_height_array[i];
-    if (std::abs(height) > std::abs(max_height)) {
-      max_height = height;
-    }
-  }
-
-  if (max_height > 0) {
-    return -normal;
-  } else {
-    return normal;
-  }
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 2);
-  using R2 = TinyVector<2, double>;
-
-  const R2 normal = this->_getNormal(mesh);
-
-  double max_height = 0;
-
-  if (m_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 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) {
-      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
-      if (std::abs(height) > std::abs(max_height)) {
-        max_height = height;
-      }
-    }
-  }
-
-  Array<double> max_height_array = parallel::allGather(max_height);
-  for (size_t i = 0; i < max_height_array.size(); ++i) {
-    const double height = max_height_array[i];
-    if (std::abs(height) > std::abs(max_height)) {
-      max_height = height;
-    }
-  }
-
-  if (max_height > 0) {
-    return -normal;
-  } else {
-    return normal;
-  }
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 3);
-  using R3 = TinyVector<3, double>;
-
-  const R3 normal = this->_getNormal(mesh);
-
-  double max_height = 0;
-
-  if (m_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 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) {
-      const double height = dot(xr[j0_nodes[r]] - xr[r0], normal);
-      if (std::abs(height) > std::abs(max_height)) {
-        max_height = height;
-      }
-    }
-  }
-
-  Array<double> max_height_array = parallel::allGather(max_height);
-  for (size_t i = 0; i < max_height_array.size(); ++i) {
-    const double height = max_height_array[i];
-    if (std::abs(height) > std::abs(max_height)) {
-      max_height = height;
-    }
-  }
-
-  if (max_height > 0) {
-    return -normal;
-  } else {
-    return normal;
-  }
-}
-
 template <size_t Dimension>
-class MeshLineNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
-{
- public:
-  using Rd = TinyVector<Dimension, double>;
-
- private:
-  const Rd m_direction;
-
-  template <typename MeshType>
-  PUGS_INLINE TinyVector<Dimension> _getDirection(const MeshType&);
-
-  template <typename MeshType>
-  PUGS_INLINE void
-  _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
-                       const TinyVector<Dimension, double>& origin,
-                       const double length,
-                       const MeshType& mesh) const
-  {
-    static_assert(MeshType::Dimension == Dimension);
-    using Rdxd = TinyMatrix<Dimension>;
-
-    const NodeValue<const Rd>& xr = mesh.xr();
-
-    Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
-
-    bool is_bad = false;
-    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
-      const Rd& x    = xr[this->m_node_list[r]];
-      const Rd delta = x - origin;
-      if (dot(P * delta, direction) > 1E-13 * length) {
-        is_bad = true;
-      }
-    });
-
-    if (parallel::allReduceOr(is_bad)) {
-      std::ostringstream ost;
-      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-          << ": boundary is not a line!";
-      throw NormalError(ost.str());
-    }
-  }
-
- public:
-  const Rd&
-  direction() const
-  {
-    return m_direction;
-  }
-
-  MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
-  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
-
-  template <typename MeshType>
-  MeshLineNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
-  {
-    ;
-  }
-
-  template <typename MeshType>
-  MeshLineNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
-  {
-    ;
-  }
-
-  template <typename MeshType>
-  MeshLineNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
-    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
-  {
-    ;
-  }
-
-  MeshLineNodeBoundary()                            = default;
-  MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
-  MeshLineNodeBoundary(MeshLineNodeBoundary&&)      = default;
-  virtual ~MeshLineNodeBoundary()                   = default;
-};
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<1>
-MeshLineNodeBoundary<1>::_getDirection(const MeshType&)
-{
-  throw UnexpectedError("MeshLineNodeBoundary::_getDirection make no sense in dimension 1");
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<2>
-MeshLineNodeBoundary<2>::_getDirection(const MeshType& mesh)
-{
-  using R2 = TinyVector<2, double>;
-
-  std::array<R2, 2> bounds = this->_getBounds(mesh);
-
-  const R2& xmin = bounds[0];
-  const R2& xmax = bounds[1];
-
-  if (xmin == xmax) {
-    std::ostringstream ost;
-    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
-        << ": unable to compute direction";
-    throw NormalError(ost.str());
-  }
-
-  R2 direction        = xmax - xmin;
-  const double length = l2Norm(direction);
-  direction *= 1. / length;
-
-  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
-
-  return direction;
-}
-
-template <>
-template <typename MeshType>
-PUGS_INLINE TinyVector<3>
-MeshLineNodeBoundary<3>::_getDirection(const MeshType& mesh)
-{
-  static_assert(MeshType::Dimension == 3);
-
-  using R3 = TinyVector<3, double>;
-
-  std::array<R3, 6> bounds = this->_getBounds(mesh);
-
-  const R3& xmin = bounds[0];
-  const R3& ymin = bounds[1];
-  const R3& zmin = bounds[2];
-  const R3& xmax = bounds[3];
-  const R3& ymax = bounds[4];
-  const R3& zmax = bounds[5];
-
-  const R3 u = xmax - xmin;
-  const R3 v = ymax - ymin;
-  const R3 w = zmax - zmin;
-
-  R3 direction = u;
-  if (dot(v, v) > dot(direction, direction)) {
-    direction = v;
-  }
-
-  if (dot(w, w) > dot(direction, direction)) {
-    direction = w;
-  }
-
-  const double length = l2Norm(direction);
-  direction *= 1. / length;
-
-  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
-
-  return direction;
-}
+MeshNodeBoundary<Dimension> getMeshNodeBoundary(const Mesh<Connectivity<Dimension>>& mesh,
+                                                const IBoundaryDescriptor& boundary_descriptor);
 
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
index b6bfe625bf6a01ab8e5e04c3b970ebc77d8ed438..fb29175a8b6254931e055d773bcf345b940de054 100644
--- a/src/mesh/MeshRandomizer.hpp
+++ b/src/mesh/MeshRandomizer.hpp
@@ -7,6 +7,8 @@
 #include <mesh/Connectivity.hpp>
 #include <mesh/ItemValueUtils.hpp>
 #include <mesh/Mesh.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <scheme/AxisBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
@@ -43,82 +45,32 @@ class MeshRandomizer
   {
     BoundaryConditionList bc_list;
 
-    constexpr ItemType EdgeType = [] {
-      if constexpr (Dimension == 3) {
-        return ItemType::edge;
-      } else if constexpr (Dimension == 2) {
-        return ItemType::face;
-      } else {
-        return ItemType::node;
-      }
-    }();
-
-    constexpr ItemType FaceType = [] {
-      if constexpr (Dimension > 1) {
-        return ItemType::face;
-      } else {
-        return ItemType::node;
-      }
-    }();
-
     for (const auto& bc_descriptor : bc_descriptor_list) {
-      bool is_valid_boundary_condition = true;
-
       switch (bc_descriptor->type()) {
       case IBoundaryConditionDescriptor::Type::axis: {
-        const AxisBoundaryConditionDescriptor& axis_bc_descriptor =
-          dynamic_cast<const AxisBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<EdgeType>();
-             ++i_ref_edge_list) {
-          const auto& ref_edge_list = mesh.connectivity().template refItemList<EdgeType>(i_ref_edge_list);
-          const RefId& ref          = ref_edge_list.refId();
-          if (ref == axis_bc_descriptor.boundaryDescriptor()) {
-            if constexpr (Dimension == 1) {
-              bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_edge_list}});
-            } else {
-              bc_list.emplace_back(AxisBoundaryCondition{MeshLineNodeBoundary<Dimension>(mesh, ref_edge_list)});
-            }
-          }
+        if constexpr (Dimension == 1) {
+          bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        } else {
+          bc_list.emplace_back(
+            AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
         }
-        is_valid_boundary_condition = true;
         break;
       }
       case IBoundaryConditionDescriptor::Type::symmetry: {
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
-             ++i_ref_face_list) {
-          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
-          const RefId& ref          = ref_face_list.refId();
-          if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-            bc_list.emplace_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
-          }
-        }
-        is_valid_boundary_condition = true;
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
         break;
       }
       case IBoundaryConditionDescriptor::Type::fixed: {
-        const FixedBoundaryConditionDescriptor& fixed_bc_descriptor =
-          dynamic_cast<const FixedBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
-             ++i_ref_face_list) {
-          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
-          const RefId& ref          = ref_face_list.refId();
-          if (ref == fixed_bc_descriptor.boundaryDescriptor()) {
-            bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_face_list}});
-          }
-        }
+        bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
         break;
       }
       default: {
-        is_valid_boundary_condition = false;
-      }
-      }
-      if (not is_valid_boundary_condition) {
         std::ostringstream error_msg;
         error_msg << *bc_descriptor << " is an invalid boundary condition for mesh randomizer";
         throw NormalError(error_msg.str());
       }
+      }
     }
 
     return bc_list;
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 0ad9b7901ceb6dd39cab12d817201e307d035308..5b4d9cb57cdea7a374b68447fc7a79ba81562e84 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -2,10 +2,13 @@
 
 #include <language/utils/InterpolateItemValue.hpp>
 #include <mesh/ItemValueUtils.hpp>
+#include <mesh/MeshFaceBoundary.hpp>
+#include <mesh/MeshFlatNodeBoundary.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
@@ -68,12 +71,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
   using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
 
+  class FixedBoundaryCondition;
   class PressureBoundaryCondition;
   class SymmetryBoundaryCondition;
   class VelocityBoundaryCondition;
 
-  using BoundaryCondition =
-    std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+  using BoundaryCondition = std::
+    variant<FixedBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
 
@@ -83,27 +87,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   NodeValue<const Rd> m_ur;
   NodeValuePerCell<const Rd> m_Fjr;
 
-  CellValue<const double>
-  _getRhoC(const DiscreteScalarFunction& rho, const DiscreteScalarFunction& c)
-  {
-    Assert(rho.mesh() == c.mesh());
-
-    const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
-
-    CellValue<double> rhoc{mesh.connectivity()};
-    parallel_for(
-      mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
-
-    return rhoc;
-  }
-
   NodeValuePerCell<const Rdxd>
-  _computeGlaceAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
-    const NodeValuePerCell<const TinyVector<Dimension>> Cjr = mesh_data.Cjr();
-    const NodeValuePerCell<const TinyVector<Dimension>> njr = mesh_data.njr();
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const NodeValuePerCell<const Rd> njr = mesh_data.njr();
 
     NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
 
@@ -120,7 +110,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   }
 
   NodeValuePerCell<const Rdxd>
-  _computeEucclhydAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  _computeEucclhydAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
   {
     MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
@@ -168,7 +158,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   }
 
   NodeValuePerCell<const Rdxd>
-  _computeAjr(const MeshType& mesh, const CellValue<const double>& rhoc)
+  _computeAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoc)
   {
     if constexpr (Dimension == 1) {
       return _computeGlaceAjr(mesh, rhoc);
@@ -264,67 +254,50 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
       switch (bc_descriptor->type()) {
       case IBoundaryConditionDescriptor::Type::symmetry: {
-        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
-          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
-             ++i_ref_face_list) {
-          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
-          const RefId& ref          = ref_face_list.refId();
-          if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-            SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
-            bc_list.push_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
-          }
-        }
-        is_valid_boundary_condition = true;
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fixed: {
+        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
         break;
       }
-
       case IBoundaryConditionDescriptor::Type::dirichlet: {
         const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
         if (dirichlet_bc_descriptor.name() == "velocity") {
-          for (size_t i_ref_face_list = 0;
-               i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
-            const RefId& ref          = ref_face_list.refId();
-            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-              MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+          MeshNodeBoundary<Dimension> mesh_node_boundary =
+            getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
-              const FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
+          Array<const Rd> value_list =
+            InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
+                                                                               mesh.xr(),
+                                                                               mesh_node_boundary.nodeList());
 
-              const auto& node_list = mesh_node_boundary.nodeList();
+          bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
+        } else if (dirichlet_bc_descriptor.name() == "pressure") {
+          const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
 
-              Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-                TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh.xr(), node_list);
+          if constexpr (Dimension == 1) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary =
+              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
 
-              bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
-            }
-          }
-        } else if (dirichlet_bc_descriptor.name() == "pressure") {
-          for (size_t i_ref_face_list = 0;
-               i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
-            const RefId& ref          = ref_face_list.refId();
-            if (ref == dirichlet_bc_descriptor.boundaryDescriptor()) {
-              const auto& face_list = ref_face_list.list();
-
-              const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
-
-              Array<const double> face_values = [&] {
-                if constexpr (Dimension == 1) {
-                  return InterpolateItemValue<double(
-                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh.xr(), face_list);
-                } else {
-                  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
-
-                  return InterpolateItemValue<double(
-                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
-                }
-              }();
-
-              bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
-            }
+            Array<const double> node_values =
+              InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
+                                                                                     mesh_node_boundary.nodeList());
+
+            bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
+          } else {
+            MeshFaceBoundary<Dimension> mesh_face_boundary =
+              getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
+            Array<const double> face_values =
+              InterpolateItemValue<double(Rd)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(),
+                                                                               mesh_face_boundary.faceList());
+            bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
           }
+
         } else {
           is_valid_boundary_condition = false;
         }
@@ -401,8 +374,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
     : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(mesh, bc_descriptor_list)}
   {
-    CellValue<const double> rhoc     = this->_getRhoC(rho, c);
-    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc);
+    NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rho * c);
 
     NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
@@ -485,10 +457,9 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
       throw NormalError("acoustic solver expects P0 functions");
     }
 
-    return this->apply(dt, dynamic_cast<const MeshType&>(*i_mesh),
-                       *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
-                       *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
-                       *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
+    return this->apply(dt, dynamic_cast<const MeshType&>(*i_mesh), dynamic_cast<const DiscreteScalarFunction&>(*rho),
+                       dynamic_cast<const DiscreteVectorFunction&>(*u),
+                       dynamic_cast<const DiscreteScalarFunction&>(*E));
   }
 
   AcousticSolver(SolverType solver_type,
@@ -500,10 +471,10 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
     : AcousticSolver{solver_type,
                      dynamic_cast<const MeshType&>(*mesh),
-                     *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
-                     *std::dynamic_pointer_cast<const DiscreteScalarFunction>(c),
-                     *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
-                     *std::dynamic_pointer_cast<const DiscreteScalarFunction>(p),
+                     dynamic_cast<const DiscreteScalarFunction&>(*rho),
+                     dynamic_cast<const DiscreteScalarFunction&>(*c),
+                     dynamic_cast<const DiscreteVectorFunction&>(*u),
+                     dynamic_cast<const DiscreteScalarFunction&>(*p),
                      bc_descriptor_list}
   {}
 
@@ -528,7 +499,7 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyPressureBC(const MeshTyp
             const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
             const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
 
-            const auto& node_list  = bc.faceList();
+            const auto& node_list  = bc.nodeList();
             const auto& value_list = bc.valueList();
             parallel_for(
               node_list.size(), PUGS_LAMBDA(size_t i_node) {
@@ -624,24 +595,53 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdx
               Ar[node_id] = identity;
               br[node_id] = value;
             });
+        } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
+          const auto& node_list = bc.nodeList();
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id = node_list[i_node];
+
+              Ar[node_id] = identity;
+              br[node_id] = zero;
+            });
         }
       },
       boundary_condition);
   }
 }
 
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::FixedBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
+    : m_mesh_node_boundary{mesh_node_boundary}
+  {}
+
+  ~FixedBoundaryCondition() = default;
+};
+
 template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryCondition
 {
  private:
+  const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
   const Array<const double> m_value_list;
-  const Array<const FaceId> m_face_list;
 
  public:
   const Array<const FaceId>&
   faceList() const
   {
-    return m_face_list;
+    return m_mesh_face_boundary.faceList();
   }
 
   const Array<const double>&
@@ -650,8 +650,9 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::PressureBoundaryConditio
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
-    : m_value_list{value_list}, m_face_list{face_list}
+  PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
+                            const Array<const double>& value_list)
+    : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
   {}
 
   ~PressureBoundaryCondition() = default;
@@ -661,14 +662,14 @@ template <>
 class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
 {
  private:
+  const MeshNodeBoundary<1> m_mesh_node_boundary;
   const Array<const double> m_value_list;
-  const Array<const NodeId> m_face_list;
 
  public:
   const Array<const NodeId>&
-  faceList() const
+  nodeList() const
   {
-    return m_face_list;
+    return m_mesh_node_boundary.nodeList();
   }
 
   const Array<const double>&
@@ -677,8 +678,8 @@ class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
-    : m_value_list{value_list}, m_face_list{face_list}
+  PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
 
   ~PressureBoundaryCondition() = default;
@@ -688,14 +689,15 @@ template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition
 {
  private:
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
   const Array<const TinyVector<Dimension>> m_value_list;
-  const Array<const NodeId> m_node_list;
 
  public:
   const Array<const NodeId>&
   nodeList() const
   {
-    return m_node_list;
+    return m_mesh_node_boundary.nodeList();
   }
 
   const Array<const TinyVector<Dimension>>&
@@ -704,8 +706,9 @@ class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryConditio
     return m_value_list;
   }
 
-  VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
-    : m_value_list{value_list}, m_node_list{node_list}
+  VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+                            const Array<const TinyVector<Dimension>>& value_list)
+    : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
   {}
 
   ~VelocityBoundaryCondition() = default;
diff --git a/src/scheme/AxisBoundaryConditionDescriptor.hpp b/src/scheme/AxisBoundaryConditionDescriptor.hpp
index a03a4bcdd5c2bfbcb23824d3da4e693d4bdbc345..07ab3be487f30af021fb51ceae84892091e3ae8b 100644
--- a/src/scheme/AxisBoundaryConditionDescriptor.hpp
+++ b/src/scheme/AxisBoundaryConditionDescriptor.hpp
@@ -1,8 +1,8 @@
 #ifndef AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 #define AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -20,7 +20,7 @@ class AxisBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
 
  public:
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/DirichletBoundaryConditionDescriptor.hpp b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
index 5fa6c42f9d4b5333ad1ca5ef2fdb46a4f06b4242..36f95df0d0ae02ea9c5b49fcd33fa8640dc7ebf1 100644
--- a/src/scheme/DirichletBoundaryConditionDescriptor.hpp
+++ b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
@@ -2,8 +2,8 @@
 #define DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
 #include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -36,7 +36,7 @@ class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   }
 
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/FixedBoundaryConditionDescriptor.hpp b/src/scheme/FixedBoundaryConditionDescriptor.hpp
index 9167c24beb7a91d1810cff528d23a4f72e39a25a..b69648e38364f80ef34552f10f7de41ceebf8f57 100644
--- a/src/scheme/FixedBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FixedBoundaryConditionDescriptor.hpp
@@ -2,8 +2,8 @@
 #define FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
 #include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -23,7 +23,7 @@ class FixedBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
 
  public:
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/FourierBoundaryConditionDescriptor.hpp b/src/scheme/FourierBoundaryConditionDescriptor.hpp
index e7aeb830c3447bcf9783060210bc97e59bf2ebdf..c9d75883dc9c8a57685db4f3477f2de67c4098ad 100644
--- a/src/scheme/FourierBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FourierBoundaryConditionDescriptor.hpp
@@ -2,8 +2,8 @@
 #define FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
 #include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -43,7 +43,7 @@ class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   }
 
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/FreeBoundaryConditionDescriptor.hpp b/src/scheme/FreeBoundaryConditionDescriptor.hpp
index b0eb2ceb218d77cae347a44219f13d278a25aee5..2d46f2a36a04cbaca0abc5b83c9ad1983e844841 100644
--- a/src/scheme/FreeBoundaryConditionDescriptor.hpp
+++ b/src/scheme/FreeBoundaryConditionDescriptor.hpp
@@ -1,8 +1,8 @@
 #ifndef FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 #define FREE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -20,7 +20,7 @@ class FreeBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
 
  public:
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index bf92cd450f000db208dfc87c340f2520b9671757..c0ea161db5d5f869d8ac8517b85ae2385af1d94a 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -3,6 +3,8 @@
 
 #include <iostream>
 
+class IBoundaryDescriptor;
+
 class IBoundaryConditionDescriptor
 {
  public:
@@ -27,6 +29,8 @@ class IBoundaryConditionDescriptor
     return bcd._write(os);
   }
 
+  virtual const IBoundaryDescriptor& boundaryDescriptor() const = 0;
+
   virtual Type type() const = 0;
 
   IBoundaryConditionDescriptor()                                    = default;
diff --git a/src/scheme/NamedBoundaryDescriptor.hpp b/src/scheme/NamedBoundaryDescriptor.hpp
index db88ee4a0dd9ac39c65ff7b1a181d3c6aad262b4..5d8959261ef37459a7ab449d3834c4f2e2050c79 100644
--- a/src/scheme/NamedBoundaryDescriptor.hpp
+++ b/src/scheme/NamedBoundaryDescriptor.hpp
@@ -1,7 +1,7 @@
 #ifndef NAMED_BOUNDARY_DESCRIPTOR_HPP
 #define NAMED_BOUNDARY_DESCRIPTOR_HPP
 
-#include <scheme/IBoundaryDescriptor.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 
 #include <iostream>
 #include <string>
diff --git a/src/scheme/NeumannBoundaryConditionDescriptor.hpp b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
index edc1b1fb35b3f530c47adcd2f0073c9727707965..597ab95e7eb5e1f61f637d840c873ba8aeb30f84 100644
--- a/src/scheme/NeumannBoundaryConditionDescriptor.hpp
+++ b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
@@ -2,8 +2,8 @@
 #define NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
 #include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -36,7 +36,7 @@ class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   }
 
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }
diff --git a/src/scheme/NumberedBoundaryDescriptor.hpp b/src/scheme/NumberedBoundaryDescriptor.hpp
index 4ed9a262859762a4103582fef57b642c82820d11..47aa7fa537f6f1070e48077cdee62aa28868782b 100644
--- a/src/scheme/NumberedBoundaryDescriptor.hpp
+++ b/src/scheme/NumberedBoundaryDescriptor.hpp
@@ -1,7 +1,7 @@
 #ifndef NUMBERED_BOUNDARY_DESCRIPTOR_HPP
 #define NUMBERED_BOUNDARY_DESCRIPTOR_HPP
 
-#include <scheme/IBoundaryDescriptor.hpp>
+#include <mesh/IBoundaryDescriptor.hpp>
 
 #include <iostream>
 
diff --git a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
index f5657ad58856f38c5d645501a14b62e77039b4c1..9364090dde18490a4803662c7eb770d9f6354b1c 100644
--- a/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
+++ b/src/scheme/SymmetryBoundaryConditionDescriptor.hpp
@@ -1,8 +1,8 @@
 #ifndef SYMMETRY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 #define SYMMETRY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
 
+#include <mesh/IBoundaryDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
 
 #include <memory>
 
@@ -20,7 +20,7 @@ class SymmetryBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
 
  public:
   const IBoundaryDescriptor&
-  boundaryDescriptor() const
+  boundaryDescriptor() const final
   {
     return *m_boundary_descriptor;
   }