diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
index 1a11de171d8172f9f276718d0624cd31b38312f4..9811c90f0b0f2c52d5f72369efea4a4398f007e3 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.cpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -1,5 +1,12 @@
 #include <language/algorithms/AcousticSolverAlgorithm.hpp>
 
+#include <language/utils/InterpolateItemValue.hpp>
+#include <output/VTKWriter.hpp>
+#include <scheme/AcousticSolver.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+
 template <size_t Dimension>
 AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
   const AcousticSolverType& solver_type,
@@ -9,9 +16,12 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
   const FunctionSymbolId& u_id,
   const FunctionSymbolId& p_id)
 {
-  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+  using MeshDataType     = MeshData<Dimension>;
+  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
 
-  std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
+  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
 
   typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
   {
@@ -24,6 +34,8 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
     }();
 
     for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
       switch (bc_descriptor->type()) {
       case IBoundaryConditionDescriptor::Type::symmetry: {
         using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
@@ -39,72 +51,74 @@ AcousticSolverAlgorithm<Dimension>::AcousticSolverAlgorithm(
               SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list)});
           }
         }
+        is_valid_boundary_condition = true;
         break;
       }
-      case IBoundaryConditionDescriptor::Type::velocity: {
-        using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
 
-        const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
-          dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*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 == velocity_bc_descriptor.boundaryDescriptor()) {
-            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+      case IBoundaryConditionDescriptor::Type::dirichlet: {
+        const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
+          dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
+        if (dirichlet_bc_descriptor.name() == "velocity") {
+          using VelocityBoundaryCondition = typename AcousticSolver<MeshType>::VelocityBoundaryCondition;
 
-            const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
+          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};
 
-            const auto& node_list = mesh_node_boundary.nodeList();
+              const FunctionSymbolId velocity_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);
+              const auto& node_list = mesh_node_boundary.nodeList();
 
-            bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
-          }
-        }
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::pressure: {
-        using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+              Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+                TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
 
-        const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
-          dynamic_cast<const PressureBoundaryConditionDescriptor&>(*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 == pressure_bc_descriptor.boundaryDescriptor()) {
-            const auto& face_list = ref_face_list.list();
-
-            const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId();
-
-            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});
+              bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+            }
+          }
+        } else if (dirichlet_bc_descriptor.name() == "pressure") {
+          using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+
+          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});
+            }
           }
+        } else {
+          is_valid_boundary_condition = false;
         }
         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 acoustic solver";
         throw NormalError(error_msg.str());
       }
-      }
     }
   }
 
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp
index bb9d9ef66d176274d305a02067e4f299c71d2c78..a6973480c4cd913c6463d39012f6eea906540003 100644
--- a/src/language/algorithms/AcousticSolverAlgorithm.hpp
+++ b/src/language/algorithms/AcousticSolverAlgorithm.hpp
@@ -2,23 +2,13 @@
 #define ACOUSTIC_SOLVER_ALGORITHM_HPP
 
 #include <language/utils/FunctionSymbolId.hpp>
-#include <language/utils/InterpolateItemValue.hpp>
-#include <output/VTKWriter.hpp>
-#include <scheme/AcousticSolver.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/AcousticSolverType.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
-#include <scheme/PressureBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <scheme/VelocityBoundaryConditionDescriptor.hpp>
 
 template <size_t Dimension>
 struct AcousticSolverAlgorithm
 {
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
-
   AcousticSolverAlgorithm(const AcousticSolverType& solver_type,
                           std::shared_ptr<const IMesh> i_mesh,
                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 0bdf7de8ec290730ff10595ba1f4866978bbb333..a6dcec9c8d510e36a8b70547eafbf11775acef8a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -4,11 +4,14 @@
 #include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/TypeDescriptor.hpp>
 #include <mesh/Mesh.hpp>
-#include <scheme/AcousticSolver.hpp>
+#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
+#include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryDescriptor.hpp>
 #include <scheme/NamedBoundaryDescriptor.hpp>
+#include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/NumberedBoundaryDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
 
 #include <memory>
 
@@ -57,7 +60,8 @@ SchemeModule::SchemeModule()
                               [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                  const FunctionSymbolId& pressure_id)
                                 -> std::shared_ptr<const IBoundaryConditionDescriptor> {
-                                return std::make_shared<PressureBoundaryConditionDescriptor>(boundary, pressure_id);
+                                return std::make_shared<DirichletBoundaryConditionDescriptor>("pressure", boundary,
+                                                                                              pressure_id);
                               }
 
                               ));
@@ -70,7 +74,8 @@ SchemeModule::SchemeModule()
                               [](std::shared_ptr<const IBoundaryDescriptor> boundary,
                                  const FunctionSymbolId& velocity_id)
                                 -> std::shared_ptr<const IBoundaryConditionDescriptor> {
-                                return std::make_shared<VelocityBoundaryConditionDescriptor>(boundary, velocity_id);
+                                return std::make_shared<DirichletBoundaryConditionDescriptor>("velocity", boundary,
+                                                                                              velocity_id);
                               }
 
                               ));
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 9ce6369e07b63537ac698149703547aa107e83cd..8189fc13a8b6a6b9948735da21456ff3296602d7 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -54,16 +54,20 @@ class SubItemValuePerItem
     const size_t m_size;
 
    public:
-    PUGS_INLINE
-    const DataType& operator[](size_t i) const noexcept(NO_ASSERT)
+    template <typename IndexType>
+    PUGS_INLINE const DataType&
+    operator[](IndexType i) const noexcept(NO_ASSERT)
     {
+      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
       Assert(i < m_size);
       return m_sub_values[i];
     }
 
-    PUGS_FORCEINLINE
-    DataType& operator[](size_t i) noexcept(NO_ASSERT)
+    template <typename IndexType>
+    PUGS_FORCEINLINE DataType&
+    operator[](IndexType i) noexcept(NO_ASSERT)
     {
+      static_assert(std::is_integral_v<IndexType>, "SubView is indexed by integral values");
       Assert(i < m_size);
       return m_sub_values[i];
     }
@@ -96,23 +100,15 @@ class SubItemValuePerItem
     return m_connectivity_ptr.use_count() != 0;
   }
 
-  // Following Kokkos logic, these classes are view and const view does allow
-  // changes in data
-  PUGS_FORCEINLINE
-  DataType&
-  operator()(ItemId j, size_t r) const noexcept(NO_ASSERT)
-  {
-    Assert(this->isBuilt());
-    return m_values[m_host_row_map(size_t{j}) + r];
-  }
-
-  template <typename IndexType>
+  template <typename IndexType, typename SubIndexType>
   PUGS_FORCEINLINE DataType&
-  operator()(IndexType j, size_t r) const noexcept(NO_ASSERT)
+  operator()(IndexType item_id, SubIndexType i) const noexcept(NO_ASSERT)
   {
-    static_assert(std::is_same_v<IndexType, ItemId>, "SubItemValuePerItem indexed by ItemId");
+    static_assert(std::is_same_v<IndexType, ItemId>, "first index must be an ItemId");
+    static_assert(std::is_integral_v<SubIndexType>, "second index must be an integral type");
     Assert(this->isBuilt());
-    return m_values[m_host_row_map(size_t{j}) + r];
+    Assert(i < m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id}));
+    return m_values[m_host_row_map(size_t{item_id}) + i];
   }
 
   PUGS_INLINE
@@ -125,19 +121,13 @@ class SubItemValuePerItem
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  PUGS_FORCEINLINE
-  DataType& operator[](size_t i) const noexcept(NO_ASSERT)
-  {
-    Assert(this->isBuilt());
-    return m_values[i];
-  }
-
-  template <typename IndexType>
-  DataType& operator[](const IndexType& i) const noexcept(NO_ASSERT)
+  template <typename ArrayIndexType>
+  DataType&
+  operator[](const ArrayIndexType& i) const noexcept(NO_ASSERT)
   {
-    static_assert(std::is_same_v<IndexType, size_t>, "Access to SubItemValuePerItem's array must be indexed by "
-                                                     "size_t");
+    static_assert(std::is_integral_v<ArrayIndexType>, "index must be an integral type");
     Assert(this->isBuilt());
+    Assert(i < m_values.size());
     return m_values[i];
   }
 
@@ -150,34 +140,40 @@ class SubItemValuePerItem
     return m_host_row_map.extent(0);
   }
 
-  PUGS_INLINE
-  size_t
-  numberOfSubValues(size_t i_cell) const noexcept(NO_ASSERT)
+  template <typename IndexType>
+  PUGS_INLINE size_t
+  numberOfSubValues(IndexType item_id) const noexcept(NO_ASSERT)
   {
+    static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
-    return m_host_row_map(i_cell + 1) - m_host_row_map(i_cell);
+    Assert(item_id < m_host_row_map.extent(0));
+    return m_host_row_map(size_t{item_id} + 1) - m_host_row_map(size_t{item_id});
   }
 
-  PUGS_INLINE
-  SubView
-  itemValues(size_t i_cell) noexcept(NO_ASSERT)
+  template <typename IndexType>
+  PUGS_INLINE SubView
+  itemValues(IndexType item_id) noexcept(NO_ASSERT)
   {
+    static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
-    const auto& cell_begin = m_host_row_map(i_cell);
-    const auto& cell_end   = m_host_row_map(i_cell + 1);
-    return SubView(m_values, cell_begin, cell_end);
+    Assert(item_id < m_host_row_map.extent(0));
+    const auto& item_begin = m_host_row_map(size_t{item_id});
+    const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
+    return SubView(m_values, item_begin, item_end);
   }
 
   // Following Kokkos logic, these classes are view and const view does allow
   // changes in data
-  PUGS_INLINE
-  SubView
-  itemValues(size_t i_cell) const noexcept(NO_ASSERT)
+  template <typename IndexType>
+  PUGS_INLINE SubView
+  itemValues(IndexType item_id) const noexcept(NO_ASSERT)
   {
+    static_assert(std::is_same_v<IndexType, ItemId>, "index must be an ItemId");
     Assert(this->isBuilt());
-    const auto& cell_begin = m_host_row_map(i_cell);
-    const auto& cell_end   = m_host_row_map(i_cell + 1);
-    return SubView(m_values, cell_begin, cell_end);
+    Assert(item_id < m_host_row_map.extent(0));
+    const auto& item_begin = m_host_row_map(size_t{item_id});
+    const auto& item_end   = m_host_row_map(size_t{item_id} + 1);
+    return SubView(m_values, item_begin, item_end);
   }
 
   template <typename DataType2, typename ConnectivityPtr2>
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 1e246f2dbf054715d67e80819078ed84b9441785..3e3151d25429593955cb281c2dd814690f67bfa5 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -9,6 +9,7 @@
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/MeshNodeBoundary.hpp>
 #include <mesh/SubItemValuePerItem.hpp>
+#include <scheme/AcousticSolverType.hpp>
 #include <scheme/BlockPerfectGas.hpp>
 #include <scheme/FiniteVolumesEulerUnknowns.hpp>
 #include <utils/ArrayUtils.hpp>
@@ -20,12 +21,6 @@
 
 #include <iostream>
 
-enum class AcousticSolverType
-{
-  Eucclhyd,
-  Glace
-};
-
 template <typename MeshType>
 class AcousticSolver
 {
diff --git a/src/scheme/AcousticSolverType.hpp b/src/scheme/AcousticSolverType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..3095dcc864d5598e985290302ad94106adbd6b2e
--- /dev/null
+++ b/src/scheme/AcousticSolverType.hpp
@@ -0,0 +1,10 @@
+#ifndef ACOUSTIC_SOLVER_TYPE_HPP
+#define ACOUSTIC_SOLVER_TYPE_HPP
+
+enum class AcousticSolverType
+{
+  Eucclhyd,
+  Glace
+};
+
+#endif   // ACOUSTIC_SOLVER_TYPE_HPP
diff --git a/src/scheme/DirichletBoundaryConditionDescriptor.hpp b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5fa6c42f9d4b5333ad1ca5ef2fdb46a4f06b4242
--- /dev/null
+++ b/src/scheme/DirichletBoundaryConditionDescriptor.hpp
@@ -0,0 +1,64 @@
+#ifndef DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class DirichletBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "dirichlet(" << m_name << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::string_view m_name;
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+  std::string_view
+  name() const
+  {
+    return m_name;
+  }
+
+  FunctionSymbolId
+  rhsSymbolId() const
+  {
+    return m_rhs_symbol_id;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::dirichlet;
+  }
+
+  DirichletBoundaryConditionDescriptor(const std::string_view name,
+                                       std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                       const FunctionSymbolId& rhs_symbol_id)
+    : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+  {
+    ;
+  }
+
+  DirichletBoundaryConditionDescriptor(const DirichletBoundaryConditionDescriptor&) = delete;
+  DirichletBoundaryConditionDescriptor(DirichletBoundaryConditionDescriptor&&)      = delete;
+
+  ~DirichletBoundaryConditionDescriptor() = default;
+};
+
+#endif   // DIRICHLET_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/FourierBoundaryConditionDescriptor.hpp b/src/scheme/FourierBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..e7aeb830c3447bcf9783060210bc97e59bf2ebdf
--- /dev/null
+++ b/src/scheme/FourierBoundaryConditionDescriptor.hpp
@@ -0,0 +1,75 @@
+#ifndef FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class FourierBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "fourier(" << m_name << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::string_view m_name;
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const FunctionSymbolId m_mass_symbol_id;
+  const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+  std::string_view
+  name() const
+  {
+    return m_name;
+  }
+
+  FunctionSymbolId
+  massSymbolId() const
+  {
+    return m_mass_symbol_id;
+  }
+
+  FunctionSymbolId
+  rhsSymbolId() const
+  {
+    return m_rhs_symbol_id;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::fourier;
+  }
+
+  FourierBoundaryConditionDescriptor(const std::string_view name,
+                                     std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                     const FunctionSymbolId& mass_symbol_id,
+                                     const FunctionSymbolId& rhs_symbol_id)
+    : m_name{name},
+      m_boundary_descriptor(boundary_descriptor),
+      m_mass_symbol_id{mass_symbol_id},
+      m_rhs_symbol_id{rhs_symbol_id}
+  {
+    ;
+  }
+
+  FourierBoundaryConditionDescriptor(const FourierBoundaryConditionDescriptor&) = delete;
+  FourierBoundaryConditionDescriptor(FourierBoundaryConditionDescriptor&&)      = delete;
+
+  ~FourierBoundaryConditionDescriptor() = default;
+};
+
+#endif   // FOURIER_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 95f68bd2575db1ce87b6657e32e1d7863297a8ae..be25ba73180f001ee264e9087379dc15ee03f378 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,8 +8,9 @@ class IBoundaryConditionDescriptor
  public:
   enum class Type
   {
-    velocity,
-    pressure,
+    dirichlet,
+    fourier,
+    neumann,
     symmetry
   };
 
diff --git a/src/scheme/NeumannBoundaryConditionDescriptor.hpp b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..edc1b1fb35b3f530c47adcd2f0073c9727707965
--- /dev/null
+++ b/src/scheme/NeumannBoundaryConditionDescriptor.hpp
@@ -0,0 +1,64 @@
+#ifndef NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class NeumannBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "neumann(" << m_name << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::string_view m_name;
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const FunctionSymbolId m_rhs_symbol_id;
+
+ public:
+  std::string_view
+  name() const
+  {
+    return m_name;
+  }
+
+  FunctionSymbolId
+  rhsSymbolId() const
+  {
+    return m_rhs_symbol_id;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::neumann;
+  }
+
+  NeumannBoundaryConditionDescriptor(const std::string_view name,
+                                     std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                     const FunctionSymbolId& rhs_symbol_id)
+    : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_rhs_symbol_id{rhs_symbol_id}
+  {
+    ;
+  }
+
+  NeumannBoundaryConditionDescriptor(const NeumannBoundaryConditionDescriptor&) = delete;
+  NeumannBoundaryConditionDescriptor(NeumannBoundaryConditionDescriptor&&)      = delete;
+
+  ~NeumannBoundaryConditionDescriptor() = default;
+};
+
+#endif   // NEUMANN_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/PressureBoundaryConditionDescriptor.hpp b/src/scheme/PressureBoundaryConditionDescriptor.hpp
deleted file mode 100644
index 8a3d39f5138e05fa41e856ebd676a0cc0c729767..0000000000000000000000000000000000000000
--- a/src/scheme/PressureBoundaryConditionDescriptor.hpp
+++ /dev/null
@@ -1,55 +0,0 @@
-#ifndef PRESSURE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-#define PRESSURE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
-
-#include <memory>
-
-class PressureBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
-{
- private:
-  std::ostream&
-  _write(std::ostream& os) const final
-  {
-    os << "pressure(" << *m_boundary_descriptor << ")";
-    return os;
-  }
-
-  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-  const FunctionSymbolId m_function_symbol_id;
-
- public:
-  FunctionSymbolId
-  functionSymbolId() const
-  {
-    return m_function_symbol_id;
-  }
-
-  const IBoundaryDescriptor&
-  boundaryDescriptor() const
-  {
-    return *m_boundary_descriptor;
-  }
-
-  Type
-  type() const final
-  {
-    return Type::pressure;
-  }
-
-  PressureBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
-                                      const FunctionSymbolId& function_symbol_id)
-    : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
-  {
-    ;
-  }
-
-  PressureBoundaryConditionDescriptor(const PressureBoundaryConditionDescriptor&) = delete;
-  PressureBoundaryConditionDescriptor(PressureBoundaryConditionDescriptor&&)      = delete;
-
-  ~PressureBoundaryConditionDescriptor() = default;
-};
-
-#endif   // PRESSURE_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/VelocityBoundaryConditionDescriptor.hpp b/src/scheme/VelocityBoundaryConditionDescriptor.hpp
deleted file mode 100644
index 5ce10e6a5ace389db0526ad3c5bd9fac31e80c9d..0000000000000000000000000000000000000000
--- a/src/scheme/VelocityBoundaryConditionDescriptor.hpp
+++ /dev/null
@@ -1,55 +0,0 @@
-#ifndef VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-#define VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
-
-#include <language/utils/FunctionSymbolId.hpp>
-#include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/IBoundaryDescriptor.hpp>
-
-#include <memory>
-
-class VelocityBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
-{
- private:
-  std::ostream&
-  _write(std::ostream& os) const final
-  {
-    os << "velocity(" << *m_boundary_descriptor << ")";
-    return os;
-  }
-
-  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
-  const FunctionSymbolId m_function_symbol_id;
-
- public:
-  FunctionSymbolId
-  functionSymbolId() const
-  {
-    return m_function_symbol_id;
-  }
-
-  const IBoundaryDescriptor&
-  boundaryDescriptor() const
-  {
-    return *m_boundary_descriptor;
-  }
-
-  Type
-  type() const final
-  {
-    return Type::velocity;
-  }
-
-  VelocityBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
-                                      const FunctionSymbolId& function_symbol_id)
-    : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
-  {
-    ;
-  }
-
-  VelocityBoundaryConditionDescriptor(const VelocityBoundaryConditionDescriptor&) = delete;
-  VelocityBoundaryConditionDescriptor(VelocityBoundaryConditionDescriptor&&)      = delete;
-
-  ~VelocityBoundaryConditionDescriptor() = default;
-};
-
-#endif   // VELOCITY!_BOUNDARY_CONDITION_DESCRIPTOR_HPP