diff --git a/CMakeLists.txt b/CMakeLists.txt
index 58fc97c458dbb7fc65fd0ebbf8d519a19f254f5e..8b5a8d153925bbe9c57dd04c817ad25ad67e01c6 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -425,6 +425,7 @@ target_link_libraries(
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsMesh
   PugsUtils
   PugsLanguageUtils
diff --git a/src/algebra/SparseMatrixDescriptor.hpp b/src/algebra/SparseMatrixDescriptor.hpp
index 3533a4017011b0c5c1f5d7f04ae2b893fa8a7d0b..100000366d3bd0e5a072084864e82cdb4b73ba29 100644
--- a/src/algebra/SparseMatrixDescriptor.hpp
+++ b/src/algebra/SparseMatrixDescriptor.hpp
@@ -28,14 +28,16 @@ class SparseMatrixDescriptor
       return m_id_value_map.size();
     }
 
-    const DataType& operator[](const IndexType& j) const
+    const DataType&
+    operator[](const IndexType& j) const
     {
       auto i_value = m_id_value_map.find(j);
       Assert(i_value != m_id_value_map.end());
       return i_value->second;
     }
 
-    DataType& operator[](const IndexType& j)
+    DataType&
+    operator[](const IndexType& j)
     {
       auto i_value = m_id_value_map.find(j);
       if (i_value != m_id_value_map.end()) {
@@ -51,7 +53,7 @@ class SparseMatrixDescriptor
     operator<<(std::ostream& os, const SparseRowDescriptor& row)
     {
       for (auto [j, value] : row.m_id_value_map) {
-        os << ' ' << j << ':' << value;
+        os << ' ' << static_cast<size_t>(j) << ':' << value;
       }
       return os;
     }
@@ -73,24 +75,30 @@ class SparseMatrixDescriptor
   SparseRowDescriptor&
   row(const IndexType i)
   {
+    Assert(i < m_row_array.size());
     return m_row_array[i];
   }
 
   const SparseRowDescriptor&
   row(const IndexType i) const
   {
+    Assert(i < m_row_array.size());
     return m_row_array[i];
   }
 
   DataType&
   operator()(const IndexType& i, const IndexType& j)
   {
+    Assert(i < m_row_array.size());
+    Assert(j < m_row_array.size());
     return m_row_array[i][j];
   }
 
   const DataType&
   operator()(const IndexType& i, const IndexType& j) const
   {
+    Assert(i < m_row_array.size());
+    Assert(j < m_row_array.size());
     const auto& r = m_row_array[i];   // split to ensure const-ness of call
     return r[j];
   }
@@ -99,7 +107,7 @@ class SparseMatrixDescriptor
   operator<<(std::ostream& os, const SparseMatrixDescriptor& M)
   {
     for (IndexType i = 0; i < M.m_row_array.size(); ++i) {
-      os << i << " |" << M.m_row_array[i] << '\n';
+      os << static_cast<size_t>(i) << " |" << M.m_row_array[i] << '\n';
     }
     return os;
   }
diff --git a/src/language/CMakeLists.txt b/src/language/CMakeLists.txt
index b9c6094e4403cb1d2289d8cfcae58de9e2e3c4e5..3565082161700da019088f7bfca70b0d9dfe6fd8 100644
--- a/src/language/CMakeLists.txt
+++ b/src/language/CMakeLists.txt
@@ -1,5 +1,6 @@
 # ------------------- Source files --------------------
 
+add_subdirectory(algorithms)
 add_subdirectory(ast)
 add_subdirectory(modules)
 add_subdirectory(node_processor)
@@ -12,6 +13,7 @@ add_library(
 # Additional dependencies
 add_dependencies(PugsLanguage
   PugsLanguage
+  PugsLanguageAlgorithms
   PugsLanguageAST
   PugsLanguageModules
   PugsLanguageUtils
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.cpp b/src/language/algorithms/AcousticSolverAlgorithm.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9811c90f0b0f2c52d5f72369efea4a4398f007e3
--- /dev/null
+++ b/src/language/algorithms/AcousticSolverAlgorithm.cpp
@@ -0,0 +1,245 @@
+#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,
+  std::shared_ptr<const IMesh> i_mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& rho_id,
+  const FunctionSymbolId& u_id,
+  const FunctionSymbolId& p_id)
+{
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+  using MeshDataType     = MeshData<Dimension>;
+  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
+
+  std::shared_ptr<const MeshType> mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+
+  typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
+  {
+    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::symmetry: {
+        using SymmetryBoundaryCondition = typename AcousticSolver<MeshType>::SymmetryBoundaryCondition;
+
+        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.push_back(
+              SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(mesh, ref_face_list)});
+          }
+        }
+        is_valid_boundary_condition = true;
+        break;
+      }
+
+      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;
+
+          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 FunctionSymbolId velocity_id = dirichlet_bc_descriptor.rhsSymbolId();
+
+              const auto& node_list = mesh_node_boundary.nodeList();
+
+              Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+                TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, mesh->xr(), node_list);
+
+              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());
+      }
+    }
+  }
+
+  UnknownsType unknowns(*mesh);
+
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    unknowns.rhoj() =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id, mesh_data.xj());
+
+    unknowns.pj() =
+      InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
+
+    unknowns.uj() =
+      InterpolateItemValue<TinyVector<Dimension>(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id,
+                                                                                                               mesh_data
+                                                                                                                 .xj());
+  }
+  unknowns.gammaj().fill(1.4);
+
+  AcousticSolver acoustic_solver(mesh, bc_list, solver_type);
+
+  const double tmax = 0.2;
+  double t          = 0;
+
+  int itermax   = std::numeric_limits<int>::max();
+  int iteration = 0;
+
+  CellValue<double>& rhoj              = unknowns.rhoj();
+  CellValue<double>& ej                = unknowns.ej();
+  CellValue<double>& pj                = unknowns.pj();
+  CellValue<double>& gammaj            = unknowns.gammaj();
+  CellValue<double>& cj                = unknowns.cj();
+  CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
+  CellValue<double>& Ej                = unknowns.Ej();
+  CellValue<double>& mj                = unknowns.mj();
+  CellValue<double>& inv_mj            = unknowns.invMj();
+
+  BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
+  block_eos.updateEandCFromRhoP();
+
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    const CellValue<const double> Vj = mesh_data.Vj();
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
+
+    parallel_for(
+      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
+  }
+
+  VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
+
+  while ((t < tmax) and (iteration < itermax)) {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    vtk_writer.write(mesh,
+                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                     t);
+
+    double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
+    if (t + dt > tmax) {
+      dt = tmax - t;
+    }
+
+    std::cout.setf(std::cout.scientific);
+    std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
+              << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
+              << rang::style::reset << '\n';
+
+    mesh = acoustic_solver.computeNextStep(dt, unknowns);
+
+    block_eos.updatePandCFromRhoE();
+
+    t += dt;
+    ++iteration;
+  }
+  std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
+            << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
+            << rang::style::reset << '\n';
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+
+    vtk_writer.write(mesh,
+                     {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
+                      NamedItemValue{"coords", mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
+                      NamedItemValue{"cell_owner", mesh->connectivity().cellOwner()},
+                      NamedItemValue{"node_owner", mesh->connectivity().nodeOwner()}},
+                     t, true);   // forces last output
+  }
+}
+
+template AcousticSolverAlgorithm<1>::AcousticSolverAlgorithm(
+  const AcousticSolverType&,
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);
+
+template AcousticSolverAlgorithm<2>::AcousticSolverAlgorithm(
+  const AcousticSolverType&,
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);
+
+template AcousticSolverAlgorithm<3>::AcousticSolverAlgorithm(
+  const AcousticSolverType&,
+  std::shared_ptr<const IMesh>,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&,
+  const FunctionSymbolId&);
diff --git a/src/language/algorithms/AcousticSolverAlgorithm.hpp b/src/language/algorithms/AcousticSolverAlgorithm.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a6973480c4cd913c6463d39012f6eea906540003
--- /dev/null
+++ b/src/language/algorithms/AcousticSolverAlgorithm.hpp
@@ -0,0 +1,20 @@
+#ifndef ACOUSTIC_SOLVER_ALGORITHM_HPP
+#define ACOUSTIC_SOLVER_ALGORITHM_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <mesh/IMesh.hpp>
+#include <scheme/AcousticSolverType.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+
+template <size_t Dimension>
+struct AcousticSolverAlgorithm
+{
+  AcousticSolverAlgorithm(const AcousticSolverType& solver_type,
+                          std::shared_ptr<const IMesh> i_mesh,
+                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                          const FunctionSymbolId& rho_id,
+                          const FunctionSymbolId& u_id,
+                          const FunctionSymbolId& p_id);
+};
+
+#endif   // ACOUSTIC_SOLVER_ALGORITHM_HPP
diff --git a/src/language/algorithms/CMakeLists.txt b/src/language/algorithms/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a8599f575058185b5ceeed0814b242ae6ec70b0a
--- /dev/null
+++ b/src/language/algorithms/CMakeLists.txt
@@ -0,0 +1,10 @@
+# ------------------- Source files --------------------
+
+add_library(PugsLanguageAlgorithms
+  AcousticSolverAlgorithm.cpp
+)
+
+
+add_dependencies(PugsLanguageAlgorithms
+  PugsUtils
+  PugsMesh)
diff --git a/src/language/modules/CMakeLists.txt b/src/language/modules/CMakeLists.txt
index 703f19857cfb99c0053aecf1af0984dab371d72e..f8cf2f0b445882a672d54e60dbc75858bbb246c7 100644
--- a/src/language/modules/CMakeLists.txt
+++ b/src/language/modules/CMakeLists.txt
@@ -11,5 +11,6 @@ add_library(PugsLanguageModules
 
 
 add_dependencies(PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsUtils
   PugsMesh)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index c04f29ad41df122df252cdc9e7009b49241e8ed2..a6dcec9c8d510e36a8b70547eafbf11775acef8a 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1,215 +1,20 @@
 #include <language/modules/SchemeModule.hpp>
 
+#include <language/algorithms/AcousticSolverAlgorithm.hpp>
 #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>
 
-/////////// TEMPORARY
-
-#include <language/utils/PugsFunctionAdapter.hpp>
-#include <output/VTKWriter.hpp>
-
-template <typename T>
-class InterpolateItemValue;
-template <typename OutputType, typename InputType>
-class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
-{
-  static constexpr size_t Dimension = OutputType::Dimension;
-  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
-
- public:
-  template <ItemType item_type>
-  static inline ItemValue<OutputType, item_type>
-  interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
-  {
-    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
-    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
-
-    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
-
-    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
-    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
-    const IConnectivity& connectivity = *position.connectivity_ptr();
-
-    ItemValue<OutputType, item_type> value(connectivity);
-    using ItemId = ItemIdT<item_type>;
-
-    parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
-      const int32_t t = tokens.acquire();
-
-      auto& execution_policy = context_list[t];
-
-      Adapter::convertArgs(execution_policy.currentContext(), position[i]);
-      auto result = expression.execute(execution_policy);
-      value[i]    = convert_result(std::move(result));
-
-      tokens.release(t);
-    });
-
-    return value;
-  }
-};
-
-template <size_t Dimension>
-struct GlaceScheme
-{
-  using ConnectivityType = Connectivity<Dimension>;
-  using MeshType         = Mesh<ConnectivityType>;
-  using MeshDataType     = MeshData<Dimension>;
-  using UnknownsType     = FiniteVolumesEulerUnknowns<MeshType>;
-
-  std::shared_ptr<const MeshType> m_mesh;
-
-  GlaceScheme(std::shared_ptr<const IMesh> i_mesh,
-              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
-              const FunctionSymbolId& rho_id,
-              const FunctionSymbolId& u_id,
-              const FunctionSymbolId& p_id)
-    : m_mesh{std::dynamic_pointer_cast<const MeshType>(i_mesh)}
-  {
-    std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
-
-    std::vector<BoundaryConditionHandler> bc_list;
-    {
-      constexpr ItemType FaceType = [] {
-        if constexpr (Dimension > 1) {
-          return ItemType::face;
-        } else {
-          return ItemType::node;
-        }
-      }();
-
-      for (const auto& bc_descriptor : bc_descriptor_list) {
-        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 < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
-            const RefId& ref          = ref_face_list.refId();
-            if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-              SymmetryBoundaryCondition<MeshType::Dimension>* sym_bc =
-                new SymmetryBoundaryCondition<MeshType::Dimension>(
-                  MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list));
-              std::shared_ptr<SymmetryBoundaryCondition<MeshType::Dimension>> bc(sym_bc);
-              bc_list.push_back(BoundaryConditionHandler(bc));
-            }
-          }
-          break;
-        }
-        default: {
-          throw UnexpectedError("Unknown BCDescription\n");
-        }
-        }
-      }
-    }
-
-    UnknownsType unknowns(*m_mesh);
-
-    {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-      unknowns.rhoj() =
-        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(rho_id,
-                                                                                                  mesh_data.xj());
-
-      unknowns.pj() =
-        InterpolateItemValue<double(TinyVector<Dimension>)>::template interpolate<ItemType::cell>(p_id, mesh_data.xj());
-
-      unknowns.uj() = InterpolateItemValue<TinyVector<Dimension>(
-        TinyVector<Dimension>)>::template interpolate<ItemType::cell>(u_id, mesh_data.xj());
-    }
-    unknowns.gammaj().fill(1.4);
-
-    AcousticSolver acoustic_solver(m_mesh, bc_list);
-
-    const double tmax = 0.2;
-    double t          = 0;
-
-    int itermax   = std::numeric_limits<int>::max();
-    int iteration = 0;
-
-    CellValue<double>& rhoj              = unknowns.rhoj();
-    CellValue<double>& ej                = unknowns.ej();
-    CellValue<double>& pj                = unknowns.pj();
-    CellValue<double>& gammaj            = unknowns.gammaj();
-    CellValue<double>& cj                = unknowns.cj();
-    CellValue<TinyVector<Dimension>>& uj = unknowns.uj();
-    CellValue<double>& Ej                = unknowns.Ej();
-    CellValue<double>& mj                = unknowns.mj();
-    CellValue<double>& inv_mj            = unknowns.invMj();
-
-    BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-    block_eos.updateEandCFromRhoP();
-
-    {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-      const CellValue<const double> Vj = mesh_data.Vj();
-
-      parallel_for(
-        m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { Ej[j] = ej[j] + 0.5 * (uj[j], uj[j]); });
-
-      parallel_for(
-        m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { mj[j] = rhoj[j] * Vj[j]; });
-
-      parallel_for(
-        m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { inv_mj[j] = 1. / mj[j]; });
-    }
-
-    VTKWriter vtk_writer("mesh_" + std::to_string(Dimension), 0.01);
-
-    while ((t < tmax) and (iteration < itermax)) {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-      vtk_writer.write(m_mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                        NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
-                       t);
-      double dt = 0.4 * acoustic_solver.acoustic_dt(mesh_data.Vj(), cj);
-      if (t + dt > tmax) {
-        dt = tmax - t;
-      }
-
-      std::cout.setf(std::cout.scientific);
-      std::cout << "iteration " << rang::fg::cyan << std::setw(4) << iteration << rang::style::reset
-                << " time=" << rang::fg::green << t << rang::style::reset << " dt=" << rang::fgB::blue << dt
-                << rang::style::reset << '\n';
-
-      m_mesh = acoustic_solver.computeNextStep(dt, unknowns);
-
-      block_eos.updatePandCFromRhoE();
-
-      t += dt;
-      ++iteration;
-    }
-    std::cout << rang::style::bold << "Final time=" << rang::fgB::green << t << rang::style::reset << " reached after "
-              << rang::fgB::cyan << iteration << rang::style::reset << rang::style::bold << " iterations"
-              << rang::style::reset << '\n';
-    {
-      MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
-
-      vtk_writer.write(m_mesh,
-                       {NamedItemValue{"density", rhoj}, NamedItemValue{"velocity", unknowns.uj()},
-                        NamedItemValue{"coords", m_mesh->xr()}, NamedItemValue{"xj", mesh_data.xj()},
-                        NamedItemValue{"cell_owner", m_mesh->connectivity().cellOwner()},
-                        NamedItemValue{"node_owner", m_mesh->connectivity().nodeOwner()}},
-                       t, true);   // forces last output
-    }
-  }
-};
-
 SchemeModule::SchemeModule()
 {
   this->_addTypeDescriptor(ast_node_data_type_from<std::shared_ptr<const IBoundaryDescriptor>>);
@@ -247,6 +52,34 @@ SchemeModule::SchemeModule()
 
                             ));
 
+  this->_addBuiltinFunction("pressure",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+                                                                  const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const FunctionSymbolId& pressure_id)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<DirichletBoundaryConditionDescriptor>("pressure", boundary,
+                                                                                              pressure_id);
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("velocity",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+                                                                  const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const FunctionSymbolId& velocity_id)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<DirichletBoundaryConditionDescriptor>("velocity", boundary,
+                                                                                              velocity_id);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("glace",
                             std::make_shared<BuiltinFunctionEmbedder<
                               void(std::shared_ptr<const IMesh>,
@@ -260,15 +93,77 @@ SchemeModule::SchemeModule()
                                  const FunctionSymbolId& p_id) -> void {
                                 switch (p_mesh->dimension()) {
                                 case 1: {
-                                  GlaceScheme<1>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 2: {
+                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                case 3: {
+                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Glace,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
+                                  break;
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("eucclhyd",
+                            std::make_shared<BuiltinFunctionEmbedder<
+                              void(std::shared_ptr<const IMesh>,
+                                   const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                   const FunctionSymbolId&, const FunctionSymbolId&, const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& rho_id, const FunctionSymbolId& u_id,
+                                 const FunctionSymbolId& p_id) -> void {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  AcousticSolverAlgorithm<1>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
                                   break;
                                 }
                                 case 2: {
-                                  GlaceScheme<2>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  AcousticSolverAlgorithm<2>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
                                   break;
                                 }
                                 case 3: {
-                                  GlaceScheme<3>{p_mesh, bc_descriptor_list, rho_id, u_id, p_id};
+                                  AcousticSolverAlgorithm<3>{AcousticSolverType::Eucclhyd,
+                                                             p_mesh,
+                                                             bc_descriptor_list,
+                                                             rho_id,
+                                                             u_id,
+                                                             p_id};
                                   break;
                                 }
                                 default: {
diff --git a/src/language/utils/InterpolateItemValue.hpp b/src/language/utils/InterpolateItemValue.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..16a87b4952c7cb0c2907833705e8142d9c0a7178
--- /dev/null
+++ b/src/language/utils/InterpolateItemValue.hpp
@@ -0,0 +1,82 @@
+#ifndef INTERPOLATE_ITEM_VALUE_HPP
+#define INTERPOLATE_ITEM_VALUE_HPP
+
+#include <language/utils/PugsFunctionAdapter.hpp>
+#include <mesh/ItemType.hpp>
+#include <mesh/ItemValue.hpp>
+
+template <typename T>
+class InterpolateItemValue;
+template <typename OutputType, typename InputType>
+class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<OutputType(InputType)>
+{
+  static constexpr size_t Dimension = OutputType::Dimension;
+  using Adapter                     = PugsFunctionAdapter<OutputType(InputType)>;
+
+ public:
+  template <ItemType item_type>
+  static inline ItemValue<OutputType, item_type>
+  interpolate(const FunctionSymbolId& function_symbol_id, const ItemValue<const InputType, item_type>& position)
+  {
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+    const IConnectivity& connectivity = *position.connectivity_ptr();
+
+    ItemValue<OutputType, item_type> value(connectivity);
+    using ItemId = ItemIdT<item_type>;
+
+    parallel_for(connectivity.template numberOf<item_type>(), [=, &expression, &tokens](ItemId i) {
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), position[i]);
+      auto result = expression.execute(execution_policy);
+      value[i]    = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return value;
+  }
+
+  template <ItemType item_type>
+  static inline Array<OutputType>
+  interpolate(const FunctionSymbolId& function_symbol_id,
+              const ItemValue<const InputType, item_type>& position,
+              const Array<const ItemIdT<item_type>>& list_of_items)
+  {
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    Array<OutputType> value{list_of_items.size()};
+    using ItemId = ItemIdT<item_type>;
+
+    parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) {
+      ItemId item_id  = list_of_items[i_item];
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), position[item_id]);
+      auto result   = expression.execute(execution_policy);
+      value[i_item] = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return value;
+  }
+};
+
+#endif   // INTERPOLATE_ITEM_VALUE_HPP
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
index e8d7071af5ae2c09b60c10478560087f5b2bdc86..1cb20d3ee595cf0963b95b4cc67f201cf2cddb5c 100644
--- a/src/mesh/ItemToItemMatrix.hpp
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -29,8 +29,10 @@ class ItemToItemMatrix
     }
 
     PUGS_INLINE
-    TargetItemId operator[](size_t j) const
+    TargetItemId
+    operator[](size_t j) const
     {
+      Assert(j < m_row.length);
       return m_row(j);
     }
 
@@ -73,15 +75,19 @@ class ItemToItemMatrix
   }
 
   PUGS_INLINE
-  auto operator[](const SourceItemId& source_id) const
+  auto
+  operator[](const SourceItemId& source_id) const
   {
+    Assert(source_id < m_connectivity_matrix.numRows());
     using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
     return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
   }
 
   template <typename IndexType>
-  PUGS_INLINE const auto& operator[](const IndexType& source_id) const
+  PUGS_INLINE const auto&
+  operator[](const IndexType& source_id) const
   {
+    Assert(source_id < m_connectivity_matrix.numRows());
     static_assert(std::is_same_v<IndexType, SourceItemId>, "ItemToItemMatrix must be indexed using correct ItemId");
     using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
     return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 86e78f5f7fde25656882535258f6fb3779855d01..aeddc582f09332543df40fbcb11e368f1d7e929e 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -32,16 +32,90 @@ class MeshData : public IMeshData
 
  private:
   const MeshType& m_mesh;
-  std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
-  std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
-  std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
-  std::shared_ptr<CellValue<const Rd>> m_xj;
-  std::shared_ptr<CellValue<const double>> m_Vj;
-  std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr;
+  NodeValuePerFace<const Rd> m_Nlr;
+  NodeValuePerFace<const Rd> m_nlr;
+  NodeValuePerCell<const Rd> m_Cjr;
+  NodeValuePerCell<const double> m_ljr;
+  NodeValuePerCell<const Rd> m_njr;
+  FaceValue<const Rd> m_xl;
+  CellValue<const Rd> m_xj;
+  CellValue<const double> m_Vj;
+  FaceValue<const double> m_ll;
 
   PUGS_INLINE
   void
-  _computeIsobarycenter()
+  _compute_ll()
+  {
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "ll does not make sense in 1d");
+    } else {
+      const auto& Nlr = this->Nlr();
+
+      FaceValue<double> ll{m_mesh.connectivity()};
+
+      const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+      parallel_for(
+        m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+          const auto& face_nodes = face_to_node_matrix[face_id];
+
+          double lenght = 0;
+          for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+            lenght += l2Norm(Nlr(face_id, i_node));
+          }
+
+          ll[face_id] = lenght;
+        });
+
+      m_ll = ll;
+    }
+  }
+
+  PUGS_INLINE
+  void
+  _compute_nlr()
+  {
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "nlr does not make sense in 1d");
+    } else {
+      const auto& Nlr = this->Nlr();
+
+      NodeValuePerFace<Rd> nlr{m_mesh.connectivity()};
+
+      parallel_for(
+        Nlr.numberOfValues(), PUGS_LAMBDA(size_t lr) {
+          double length = l2Norm(Nlr[lr]);
+          nlr[lr]       = 1. / length * Nlr[lr];
+        });
+
+      m_nlr = nlr;
+    }
+  }
+
+  PUGS_INLINE void
+  _computeFaceIsobarycenter()
+  {   // Computes vertices isobarycenter
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "xl does not make sense in 1d");
+    } else {
+      const NodeValue<const Rd>& xr = m_mesh.xr();
+
+      const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+      FaceValue<Rd> xl(m_mesh.connectivity());
+      parallel_for(
+        m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+          Rd X                   = zero;
+          const auto& face_nodes = face_to_node_matrix[face_id];
+          for (size_t R = 0; R < face_nodes.size(); ++R) {
+            X += xr[face_nodes[R]];
+          }
+          xl[face_id] = 1. / face_nodes.size() * X;
+        });
+      m_xl = xl;
+    }
+  }
+
+  void
+  _computeCellIsobarycenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
       const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -54,7 +128,7 @@ class MeshData : public IMeshData
           const auto& cell_nodes = cell_to_node_matrix[j];
           xj[j]                  = 0.5 * (xr[cell_nodes[0]] + xr[cell_nodes[1]]);
         });
-      m_xj = std::make_shared<CellValue<const Rd>>(xj);
+      m_xj = xj;
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
@@ -71,7 +145,7 @@ class MeshData : public IMeshData
           }
           xj[j] = inv_cell_nb_nodes[j] * X;
         });
-      m_xj = std::make_shared<CellValue<const Rd>>(xj);
+      m_xj = xj;
     }
   }
 
@@ -95,7 +169,7 @@ class MeshData : public IMeshData
         }
         Vj[j] = inv_Dimension * sum_cjr_xr;
       });
-    m_Vj = std::make_shared<CellValue<const double>>(Vj);
+    m_Vj = Vj;
   }
 
   PUGS_INLINE
@@ -123,7 +197,7 @@ class MeshData : public IMeshData
           Nlr(l, 0) = Nr;
           Nlr(l, 1) = Nr;
         });
-      m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+      m_Nlr = Nlr;
     } else {
       const NodeValue<const Rd>& xr = m_mesh.xr();
 
@@ -150,7 +224,7 @@ class MeshData : public IMeshData
             Nlr(l, r) = Nr;
           }
         });
-      m_Nlr = std::make_shared<NodeValuePerFace<const Rd>>(Nlr);
+      m_Nlr = Nlr;
     }
   }
 
@@ -165,7 +239,7 @@ class MeshData : public IMeshData
           Cjr(j, 0) = -1;
           Cjr(j, 1) = 1;
         });
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     } else if constexpr (Dimension == 2) {
       const NodeValue<const Rd>& xr   = m_mesh.xr();
       const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
@@ -181,7 +255,7 @@ class MeshData : public IMeshData
             Cjr(j, R)       = Rd{half_xrp_xrm[1], -half_xrp_xrm[0]};
           }
         });
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     } else if (Dimension == 3) {
       auto Nlr = this->Nlr();
 
@@ -228,7 +302,7 @@ class MeshData : public IMeshData
           }
         });
 
-      m_Cjr = std::make_shared<NodeValuePerCell<const Rd>>(Cjr);
+      m_Cjr = Cjr;
     }
   }
 
@@ -241,13 +315,13 @@ class MeshData : public IMeshData
       NodeValuePerCell<double> ljr(m_mesh.connectivity());
       parallel_for(
         ljr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = 1; });
-      m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
+      m_ljr = ljr;
 
     } else {
       NodeValuePerCell<double> ljr(m_mesh.connectivity());
       parallel_for(
         Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { ljr[jr] = l2Norm(Cjr[jr]); });
-      m_ljr = std::make_shared<NodeValuePerCell<const double>>(ljr);
+      m_ljr = ljr;
     }
   }
 
@@ -265,19 +339,18 @@ class MeshData : public IMeshData
       NodeValuePerCell<Rd> njr(m_mesh.connectivity());
       parallel_for(
         Cjr.numberOfValues(), PUGS_LAMBDA(size_t jr) { njr[jr] = (1. / ljr[jr]) * Cjr[jr]; });
-      m_njr = std::make_shared<NodeValuePerCell<const Rd>>(njr);
+      m_njr = njr;
     }
   }
 
   void
   _checkCellVolume()
   {
-    Assert(m_Vj);
-    auto& Vj = *m_Vj;
+    Assert(m_Vj.isBuilt());
 
     bool is_valid = [&] {
       for (CellId j = 0; j < m_mesh.numberOfCells(); ++j) {
-        if (Vj[j] <= 0) {
+        if (m_Vj[j] <= 0) {
           return false;
         }
       }
@@ -297,65 +370,95 @@ class MeshData : public IMeshData
     return m_mesh;
   }
 
+  PUGS_INLINE
+  FaceValue<const double>
+  ll()
+  {
+    if (not m_ll.isBuilt()) {
+      this->_compute_ll();
+    }
+    return m_ll;
+  }
+
   PUGS_INLINE
   NodeValuePerFace<const Rd>
   Nlr()
   {
-    if (not m_Nlr) {
+    if (not m_Nlr.isBuilt()) {
       this->_computeNlr();
     }
-    return *m_Nlr;
+    return m_Nlr;
+  }
+
+  PUGS_INLINE
+  NodeValuePerFace<const Rd>
+  nlr()
+  {
+    if (not m_nlr.isBuilt()) {
+      this->_compute_nlr();
+    }
+    return m_nlr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const Rd>
   Cjr()
   {
-    if (not m_Cjr) {
+    if (not m_Cjr.isBuilt()) {
       this->_computeCjr();
     }
-    return *m_Cjr;
+    return m_Cjr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const double>
   ljr()
   {
-    if (not m_ljr) {
+    if (not m_ljr.isBuilt()) {
       this->_compute_ljr();
     }
-    return *m_ljr;
+    return m_ljr;
   }
 
   PUGS_INLINE
   NodeValuePerCell<const Rd>
   njr()
   {
-    if (not m_njr) {
+    if (not m_njr.isBuilt()) {
       this->_compute_njr();
     }
-    return *m_njr;
+    return m_njr;
+  }
+
+  PUGS_INLINE
+  FaceValue<const Rd>
+  xl()
+  {
+    if (not m_xl.isBuilt()) {
+      this->_computeFaceIsobarycenter();
+    }
+    return m_xl;
   }
 
   PUGS_INLINE
   CellValue<const Rd>
   xj()
   {
-    if (not m_xj) {
-      this->_computeIsobarycenter();
+    if (not m_xj.isBuilt()) {
+      this->_computeCellIsobarycenter();
     }
-    return *m_xj;
+    return m_xj;
   }
 
   PUGS_INLINE
   CellValue<const double>
   Vj()
   {
-    if (not m_Vj) {
+    if (not m_Vj.isBuilt()) {
       this->_computeCellVolume();
       this->_checkCellVolume();
     }
-    return *m_Vj;
+    return m_Vj;
   }
 
  private:
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 d0eb4d0b699709e0cddff74035e6b2263b7fffd2..3e3151d25429593955cb281c2dd814690f67bfa5 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,47 +1,170 @@
 #ifndef ACOUSTIC_SOLVER_HPP
 #define ACOUSTIC_SOLVER_HPP
 
-#include <rang.hpp>
-
-#include <utils/ArrayUtils.hpp>
-
-#include <scheme/BlockPerfectGas.hpp>
-#include <utils/PugsAssert.hpp>
-
-#include <scheme/BoundaryCondition.hpp>
-#include <scheme/FiniteVolumesEulerUnknowns.hpp>
-
-#include <mesh/Mesh.hpp>
-#include <mesh/MeshData.hpp>
-#include <mesh/MeshDataManager.hpp>
-
 #include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
-
 #include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshData.hpp>
+#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>
 #include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
+
+#include <rang.hpp>
 
 #include <iostream>
 
 template <typename MeshType>
 class AcousticSolver
 {
+ public:
+  class PressureBoundaryCondition;
+  class SymmetryBoundaryCondition;
+  class VelocityBoundaryCondition;
+
+  using BoundaryCondition =
+    std::variant<PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+
   constexpr static size_t Dimension = MeshType::Dimension;
 
+ private:
   using MeshDataType = MeshData<Dimension>;
   using UnknownsType = FiniteVolumesEulerUnknowns<MeshType>;
 
   std::shared_ptr<const MeshType> m_mesh;
   const typename MeshType::Connectivity& m_connectivity;
-  const std::vector<BoundaryConditionHandler>& m_boundary_condition_list;
 
   using Rd  = TinyVector<Dimension>;
   using Rdd = TinyMatrix<Dimension>;
 
- private:
+  const BoundaryConditionList m_boundary_condition_list;
+
+  const AcousticSolverType m_solver_type;
+
+  void
+  _applyPressureBC()
+  {
+    for (const auto& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
+            MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+            if constexpr (Dimension == 1) {
+              const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+
+              const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
+              const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
+
+              const auto& node_list  = bc.faceList();
+              const auto& value_list = bc.valueList();
+              parallel_for(
+                node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                  const NodeId node_id       = node_list[i_node];
+                  const auto& node_cell_list = node_to_cell_matrix[node_id];
+                  Assert(node_cell_list.size() == 1);
+
+                  CellId node_cell_id              = node_cell_list[0];
+                  size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
+
+                  m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
+                });
+            } else {
+              const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+              const auto& face_to_cell_matrix               = m_connectivity.faceToCellMatrix();
+              const auto& face_to_node_matrix               = m_connectivity.faceToNodeMatrix();
+              const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
+              const auto& face_cell_is_reversed             = m_connectivity.cellFaceIsReversed();
+
+              const auto& face_list  = bc.faceList();
+              const auto& value_list = bc.valueList();
+              for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                const FaceId face_id       = face_list[i_face];
+                const auto& face_cell_list = face_to_cell_matrix[face_id];
+                Assert(face_cell_list.size() == 1);
+
+                CellId face_cell_id              = face_cell_list[0];
+                size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+                const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+                const auto& face_nodes = face_to_node_matrix[face_id];
+
+                for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                  NodeId node_id = face_nodes[i_node];
+                  m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+                }
+              }
+            }
+          }
+        },
+        boundary_condition);
+    }
+  }
+
+  void
+  _applySymmetryBC()
+  {
+    for (const auto& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
+            const Rd& n = bc.outgoingNormal();
+
+            const Rdd I   = identity;
+            const Rdd nxn = tensorProduct(n, n);
+            const Rdd P   = I - nxn;
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+                const NodeId r = node_list[r_number];
+
+                m_Ar[r] = P * m_Ar[r] * P + nxn;
+                m_br[r] = P * m_br[r];
+              });
+          }
+        },
+        boundary_condition);
+    }
+  }
+
+  void
+  _applyVelocityBC()
+  {
+    for (const auto& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using T = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
+            const auto& node_list  = bc.nodeList();
+            const auto& value_list = bc.valueList();
+
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                NodeId node_id    = node_list[i_node];
+                const auto& value = value_list[i_node];
+
+                m_Ar[node_id] = identity;
+                m_br[node_id] = value;
+              });
+          }
+        },
+        boundary_condition);
+    }
+  }
+
   PUGS_INLINE
   const CellValue<const double>
   computeRhoCj(const CellValue<const double>& rhoj, const CellValue<const double>& cj)
@@ -53,11 +176,13 @@ class AcousticSolver
 
   PUGS_INLINE
   void
-  computeAjr(const CellValue<const double>& rhocj,
-             const NodeValuePerCell<const Rd>& Cjr,
-             const NodeValuePerCell<const double>& /* ljr */,
-             const NodeValuePerCell<const Rd>& njr)
+  _computeGlaceAjr(const CellValue<const double>& rhocj)
   {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
+    const NodeValuePerCell<const Rd> njr = mesh_data.njr();
+
     parallel_for(
       m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
         const size_t& nb_nodes = m_Ajr.numberOfSubValues(j);
@@ -68,6 +193,71 @@ class AcousticSolver
       });
   }
 
+  PUGS_INLINE
+  void
+  _computeEucclhydAjr(const CellValue<const double>& rhocj)
+  {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+    const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
+
+    const auto& face_to_node_matrix = m_connectivity.faceToNodeMatrix();
+    const auto& cell_to_node_matrix = m_connectivity.cellToNodeMatrix();
+    const auto& cell_to_face_matrix = m_connectivity.cellToFaceMatrix();
+
+    parallel_for(
+      m_Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { m_Ajr[jr] = zero; });
+
+    parallel_for(
+      m_mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+        const auto& cell_nodes = cell_to_node_matrix[j];
+
+        const auto& cell_faces = cell_to_face_matrix[j];
+
+        const double& rho_c = rhocj[j];
+
+        for (size_t L = 0; L < cell_faces.size(); ++L) {
+          const FaceId& l        = cell_faces[L];
+          const auto& face_nodes = face_to_node_matrix[l];
+
+          auto local_node_number_in_cell = [&](NodeId node_number) {
+            for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
+              if (node_number == cell_nodes[i_node]) {
+                return i_node;
+              }
+            }
+            return std::numeric_limits<size_t>::max();
+          };
+
+          for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
+            const size_t R = local_node_number_in_cell(face_nodes[rl]);
+            m_Ajr(j, R) += tensorProduct(rho_c * Nlr(l, rl), nlr(l, rl));
+          }
+        }
+      });
+  }
+
+  PUGS_INLINE
+  void
+  computeAjr(const CellValue<const double>& rhocj)
+  {
+    switch (m_solver_type) {
+    case AcousticSolverType::Glace: {
+      this->_computeGlaceAjr(rhocj);
+      break;
+    }
+    case AcousticSolverType::Eucclhyd: {
+      if constexpr (Dimension > 1) {
+        this->_computeEucclhydAjr(rhocj);
+      } else {
+        this->_computeGlaceAjr(rhocj);
+      }
+      break;
+    }
+    }
+  }
+
   PUGS_INLINE
   const NodeValue<const Rdd>
   computeAr(const NodeValuePerCell<const Rdd>& Ajr)
@@ -94,11 +284,13 @@ class AcousticSolver
 
   PUGS_INLINE
   const NodeValue<const Rd>
-  computeBr(const NodeValuePerCell<Rdd>& Ajr,
-            const NodeValuePerCell<const Rd>& Cjr,
-            const CellValue<const Rd>& uj,
-            const CellValue<const double>& pj)
+  computeBr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
   {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValuePerCell<const Rd>& Cjr  = mesh_data.Cjr();
+    const NodeValuePerCell<const Rdd>& Ajr = m_Ajr;
+
     const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
     const auto& node_local_numbers_in_their_cells = m_connectivity.nodeLocalNumbersInTheirCells();
 
@@ -121,58 +313,33 @@ class AcousticSolver
   void
   applyBoundaryConditions()
   {
-    for (const auto& handler : m_boundary_condition_list) {
-      switch (handler.boundaryCondition().type()) {
-      case BoundaryCondition::normal_velocity: {
-        throw NotImplementedError("normal_velocity BC");
-      }
-      case BoundaryCondition::velocity: {
-        throw NotImplementedError("velocity BC");
-      }
-      case BoundaryCondition::pressure: {
-        throw NotImplementedError("pressure BC");
-      }
-      case BoundaryCondition::symmetry: {
-        const SymmetryBoundaryCondition<Dimension>& symmetry_bc =
-          dynamic_cast<const SymmetryBoundaryCondition<Dimension>&>(handler.boundaryCondition());
-        const Rd& n = symmetry_bc.outgoingNormal();
-
-        const Rdd I   = identity;
-        const Rdd nxn = tensorProduct(n, n);
-        const Rdd P   = I - nxn;
-
-        const Array<const NodeId>& node_list = symmetry_bc.nodeList();
-        parallel_for(
-          symmetry_bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
-            const NodeId r = node_list[r_number];
-
-            m_Ar[r] = P * m_Ar[r] * P + nxn;
-            m_br[r] = P * m_br[r];
-          });
-        break;
-      }
-      }
-    }
+    this->_applyPressureBC();
+    this->_applySymmetryBC();
+    this->_applyVelocityBC();
   }
 
-  NodeValue<Rd>
-  computeUr(const NodeValue<const Rdd>& Ar, const NodeValue<const Rd>& br)
+  void
+  computeUr()
   {
+    const NodeValue<const Rdd> Ar = m_Ar;
+    const NodeValue<const Rd> br  = m_br;
+
     inverse(Ar, m_inv_Ar);
     const NodeValue<const Rdd> invAr = m_inv_Ar;
     parallel_for(
       m_mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { m_ur[r] = invAr[r] * br[r]; });
-
-    return m_ur;
   }
 
   void
-  computeFjr(const NodeValuePerCell<Rdd>& Ajr,
-             const NodeValue<const Rd>& ur,
-             const NodeValuePerCell<const Rd>& Cjr,
-             const CellValue<const Rd>& uj,
-             const CellValue<const double>& pj)
+  computeFjr(const CellValue<const Rd>& uj, const CellValue<const double>& pj)
   {
+    MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+    const NodeValuePerCell<const Rd> Cjr  = mesh_data.Cjr();
+    const NodeValuePerCell<const Rdd> Ajr = m_Ajr;
+
+    const NodeValue<const Rd> ur = m_ur;
+
     const auto& cell_to_node_matrix = m_mesh->connectivity().cellToNodeMatrix();
 
     parallel_for(
@@ -197,26 +364,22 @@ class AcousticSolver
   computeExplicitFluxes(const CellValue<const double>& rhoj,
                         const CellValue<const Rd>& uj,
                         const CellValue<const double>& pj,
-                        const CellValue<const double>& cj,
-                        const NodeValuePerCell<const Rd>& Cjr,
-                        const NodeValuePerCell<const double>& ljr,
-                        const NodeValuePerCell<const Rd>& njr)
+                        const CellValue<const double>& cj)
   {
     const CellValue<const double> rhocj = computeRhoCj(rhoj, cj);
-    computeAjr(rhocj, Cjr, ljr, njr);
+    computeAjr(rhocj);
 
     NodeValuePerCell<const Rdd> Ajr = m_Ajr;
     this->computeAr(Ajr);
-    this->computeBr(m_Ajr, Cjr, uj, pj);
+    this->computeBr(uj, pj);
 
     this->applyBoundaryConditions();
 
     synchronize(m_Ar);
     synchronize(m_br);
 
-    NodeValue<Rd>& ur = m_ur;
-    ur                = computeUr(m_Ar, m_br);
-    computeFjr(m_Ajr, ur, Cjr, uj, pj);
+    computeUr();
+    computeFjr(uj, pj);
   }
 
   NodeValue<Rd> m_br;
@@ -229,10 +392,13 @@ class AcousticSolver
   CellValue<double> m_Vj_over_cj;
 
  public:
-  AcousticSolver(std::shared_ptr<const MeshType> p_mesh, const std::vector<BoundaryConditionHandler>& bc_list)
+  AcousticSolver(std::shared_ptr<const MeshType> p_mesh,
+                 const BoundaryConditionList& bc_list,
+                 const AcousticSolverType solver_type)
     : m_mesh(p_mesh),
       m_connectivity(m_mesh->connectivity()),
       m_boundary_condition_list(bc_list),
+      m_solver_type(solver_type),
       m_br(m_connectivity),
       m_Ajr(m_connectivity),
       m_Ar(m_connectivity),
@@ -281,12 +447,9 @@ class AcousticSolver
 
     MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
 
-    const CellValue<const double> Vj         = mesh_data.Vj();
-    const NodeValuePerCell<const Rd> Cjr     = mesh_data.Cjr();
-    const NodeValuePerCell<const double> ljr = mesh_data.ljr();
-    const NodeValuePerCell<const Rd> njr     = mesh_data.njr();
+    const CellValue<const double> Vj = mesh_data.Vj();
 
-    computeExplicitFluxes(rhoj, uj, pj, cj, Cjr, ljr, njr);
+    computeExplicitFluxes(rhoj, uj, pj, cj);
 
     const NodeValuePerCell<Rd>& Fjr = m_Fjr;
     const NodeValue<const Rd> ur    = m_ur;
@@ -326,4 +489,124 @@ class AcousticSolver
   }
 };
 
+template <typename MeshType>
+class AcousticSolver<MeshType>::PressureBoundaryCondition
+{
+ private:
+  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;
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    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() = default;
+};
+
+template <>
+class AcousticSolver<Mesh<Connectivity<1>>>::PressureBoundaryCondition
+{
+ private:
+  const Array<const double> m_value_list;
+  const Array<const NodeId> m_face_list;
+
+ public:
+  const Array<const NodeId>&
+  faceList() const
+  {
+    return m_face_list;
+  }
+
+  const Array<const double>&
+  valueList() const
+  {
+    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() = default;
+};
+
+template <typename MeshType>
+class AcousticSolver<MeshType>::VelocityBoundaryCondition
+{
+ private:
+  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;
+  }
+
+  const Array<const TinyVector<Dimension>>&
+  valueList() const
+  {
+    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() = default;
+};
+
+template <typename MeshType>
+class AcousticSolver<MeshType>::SymmetryBoundaryCondition
+{
+ public:
+  static constexpr size_t Dimension = MeshType::Dimension;
+
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
+
+ public:
+  const Rd&
+  outgoingNormal() const
+  {
+    return m_mesh_flat_node_boundary.outgoingNormal();
+  }
+
+  size_t
+  numberOfNodes() const
+  {
+    return m_mesh_flat_node_boundary.nodeList().size();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
 #endif   // ACOUSTIC_SOLVER_HPP
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/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
deleted file mode 100644
index 3d3332fc92cfa88792ecf556cb97781de007a104..0000000000000000000000000000000000000000
--- a/src/scheme/BoundaryCondition.hpp
+++ /dev/null
@@ -1,146 +0,0 @@
-#ifndef BOUNDARY_CONDITION_HANDLER_HPP
-#define BOUNDARY_CONDITION_HANDLER_HPP
-
-#include <utils/Array.hpp>
-
-#include <mesh/MeshNodeBoundary.hpp>
-#include <mesh/RefItemList.hpp>
-
-#include <memory>
-#include <vector>
-
-class BoundaryCondition   // clazy:exclude=copyable-polymorphic
-{
- public:
-  enum Type
-  {
-    velocity,
-    normal_velocity,
-    pressure,
-    symmetry
-  };
-
-  const Type m_type;
-
- protected:
-  BoundaryCondition(Type type) : m_type(type)
-  {
-    ;
-  }
-
- public:
-  const Type&
-  type() const
-  {
-    return m_type;
-  }
-
-  virtual ~BoundaryCondition() = default;
-};
-
-class PressureBoundaryCondition final : public BoundaryCondition   // clazy:exclude=copyable-polymorphic
-{
- private:
-  const double m_value;
-
-  const size_t m_number_of_faces;
-  Array<const unsigned int> m_face_list;
-
- public:
-  double
-  value() const
-  {
-    return m_value;
-  }
-
-  const size_t&
-  numberOfFaces() const
-  {
-    return m_number_of_faces;
-  }
-
-  const Array<const unsigned int>&
-  faceList() const
-  {
-    return m_face_list;
-  }
-
-  PressureBoundaryCondition(double value, const std::vector<unsigned int>& faces)
-    : BoundaryCondition(BoundaryCondition::pressure), m_value(value), m_number_of_faces(faces.size())
-  {
-    Array<unsigned int> face_list(faces.size());
-    parallel_for(
-      m_number_of_faces, PUGS_LAMBDA(int f) { face_list[f] = faces[f]; });
-    m_face_list = face_list;
-  }
-
-  ~PressureBoundaryCondition() = default;
-};
-
-template <size_t dimension>
-class SymmetryBoundaryCondition : public BoundaryCondition
-{
- public:
-  using Rd = TinyVector<dimension, double>;
-
- private:
-  const MeshFlatNodeBoundary<dimension> m_mesh_flat_node_boundary;
-
- public:
-  const Rd&
-  outgoingNormal() const
-  {
-    return m_mesh_flat_node_boundary.outgoingNormal();
-  }
-
-  size_t
-  numberOfNodes() const
-  {
-    return m_mesh_flat_node_boundary.nodeList().size();
-  }
-
-  const Array<const NodeId>&
-  nodeList() const
-  {
-    return m_mesh_flat_node_boundary.nodeList();
-  }
-
-  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<dimension>& mesh_flat_node_boundary)
-    : BoundaryCondition(BoundaryCondition::symmetry), m_mesh_flat_node_boundary(mesh_flat_node_boundary)
-  {
-    ;
-  }
-
-  ~SymmetryBoundaryCondition() = default;
-};
-
-class BoundaryConditionHandler
-{
- private:
-  std::shared_ptr<BoundaryCondition> m_boundary_condition;
-
- public:
-  const BoundaryCondition&
-  boundaryCondition() const
-  {
-    return *m_boundary_condition;
-  }
-
-  PUGS_INLINE
-  BoundaryConditionHandler& operator=(BoundaryConditionHandler&) = default;
-
-  PUGS_INLINE
-  BoundaryConditionHandler(const BoundaryConditionHandler&) = default;
-
-  PUGS_INLINE
-  BoundaryConditionHandler(BoundaryConditionHandler&&) = default;
-
-  PUGS_INLINE
-  BoundaryConditionHandler(std::shared_ptr<BoundaryCondition> boundary_condition)
-    : m_boundary_condition(boundary_condition)
-  {}
-
-  ~BoundaryConditionHandler() = default;
-};
-
-#endif   // BOUNDARY_CONDITION_HANDLER_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 0814c9e8580f6c36ad1152c3143d5877bff45dbc..be25ba73180f001ee264e9087379dc15ee03f378 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,6 +8,9 @@ class IBoundaryConditionDescriptor
  public:
   enum class Type
   {
+    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/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2dfd3234486b352e5b9f27e13e10001846008bc1..7cd1ed89023e619af478625b58441002ecae0af2 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -90,6 +90,7 @@ target_link_libraries (unit_tests
   PugsLanguage
   PugsLanguageAST
   PugsLanguageModules
+  PugsLanguageAlgorithms
   PugsLanguageUtils
   PugsMesh
   PugsUtils
diff --git a/tests/test_SparseMatrixDescriptor.cpp b/tests/test_SparseMatrixDescriptor.cpp
index 4e8b90134060c786cc7f97bc8fc8cd00bfba3257..cd66380b1d1beb6bcb9ba1a9c0d9c495947a3235 100644
--- a/tests/test_SparseMatrixDescriptor.cpp
+++ b/tests/test_SparseMatrixDescriptor.cpp
@@ -168,4 +168,30 @@ TEST_CASE("SparseMatrixDescriptor", "[algebra]")
     REQUIRE(value_array[6] == 1);
     REQUIRE(value_array[7] == -2);
   }
+
+  SECTION("output")
+  {
+    SparseMatrixDescriptor<int, uint8_t> S{5};
+    S(0, 2) = 3;
+    S(1, 2) = 11;
+    S(1, 1) = 1;
+    S(0, 2) += 2;
+    S(3, 3) = 5;
+    S(3, 1) = -3;
+    S(4, 1) = 1;
+    S(2, 2) = 4;
+    S(4, 4) = -2;
+
+    std::ostringstream output;
+    output << '\n' << S;
+
+    std::string expected_output = R"(
+0 | 2:5
+1 | 1:1 2:11
+2 | 2:4
+3 | 1:-3 3:5
+4 | 1:1 4:-2
+)";
+    REQUIRE(output.str() == expected_output);
+  }
 }