diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 2d429282e7ee7963121b1ee827a3ae3903a00449..97b03eda3487e6506c2d18f77a0ac5fd13d9b780 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -17,6 +17,7 @@ add_library(
   LogicalConnectivityBuilder.cpp
   MeshBuilderBase.cpp
   MeshDataManager.cpp
+  MeshFaceBoundary.cpp
   MeshFlatNodeBoundary.cpp
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
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/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index dedfbd2af950547841e13b5c953ee069df878ffe..5b4d9cb57cdea7a374b68447fc7a79ba81562e84 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -2,6 +2,7 @@
 
 #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>
@@ -91,8 +92,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   {
     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()};
 
@@ -253,12 +254,12 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
 
       switch (bc_descriptor->type()) {
       case IBoundaryConditionDescriptor::Type::symmetry: {
-        bc_list.push_back(
+        bc_list.emplace_back(
           SymmetryBoundaryCondition(getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
         break;
       }
       case IBoundaryConditionDescriptor::Type::fixed: {
-        bc_list.push_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+        bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
         break;
       }
       case IBoundaryConditionDescriptor::Type::dirichlet: {
@@ -268,36 +269,35 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
           MeshNodeBoundary<Dimension> mesh_node_boundary =
             getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
 
-          Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-            TinyVector<Dimension>)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
-                                                                          mesh.xr(), mesh_node_boundary.nodeList());
+          Array<const Rd> value_list =
+            InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
+                                                                               mesh.xr(),
+                                                                               mesh_node_boundary.nodeList());
 
-          bc_list.push_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
+          bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, 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});
-            }
+          const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+          if constexpr (Dimension == 1) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary =
+              getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
+
+            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;
         }
@@ -499,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) {
@@ -634,14 +634,14 @@ 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;