From 96e80be056339c041f83ffc4045e608d892d7d03 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 28 Jul 2020 16:14:10 +0200
Subject: [PATCH] Finalize boundary condition rewriting for acoustic solver
 (Glace)

Now boundary condition treatment is embedded within the scheme which
makes lots more sense
---
 src/language/modules/SchemeModule.cpp |  91 ++++---
 src/scheme/AcousticSolver.hpp         | 375 ++++++++++++++++++--------
 src/scheme/BoundaryCondition.hpp      | 187 -------------
 3 files changed, 311 insertions(+), 342 deletions(-)
 delete mode 100644 src/scheme/BoundaryCondition.hpp

diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 4e19eb473..2491767f4 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -112,7 +112,7 @@ struct GlaceScheme
   {
     std::cout << "number of bc descr = " << bc_descriptor_list.size() << '\n';
 
-    std::vector<BoundaryConditionHandler> bc_list;
+    typename AcousticSolver<MeshType>::BoundaryConditionList bc_list;
     {
       constexpr ItemType FaceType = [] {
         if constexpr (Dimension > 1) {
@@ -125,6 +125,8 @@ struct GlaceScheme
       for (const auto& bc_descriptor : bc_descriptor_list) {
         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;
@@ -132,16 +134,15 @@ struct GlaceScheme
             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));
+              bc_list.push_back(
+                SymmetryBoundaryCondition{MeshFlatNodeBoundary<MeshType::Dimension>(m_mesh, ref_face_list)});
             }
           }
           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;
@@ -149,43 +150,51 @@ struct GlaceScheme
             const auto& ref_face_list = m_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{m_mesh, ref_face_list};
+
               const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
 
-              if constexpr (Dimension == 1) {
-                const auto& node_list = ref_face_list.list();
-
-                Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-                  TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
-
-                std::shared_ptr bc =
-                  std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
-                bc_list.push_back(BoundaryConditionHandler(bc));
-              } else {
-                const auto& face_list           = ref_face_list.list();
-                const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
-                std::set<NodeId> node_set;
-                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
-                  FaceId face_id         = face_list[i_face];
-                  const auto& face_nodes = face_to_node_matrix[face_id];
-                  for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
-                    node_set.insert(face_nodes[i_node]);
-                  }
-                }
+              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, m_mesh->xr(), node_list);
+
+              bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+
+              // if constexpr (Dimension == 1) {
+              //   const auto& node_list = ref_face_list.list();
+
+              //   Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+              //     TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
 
-                Array<NodeId> node_list = convert_to_array(node_set);
+              //   bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+              // } else {
+              //   const auto& face_list           = ref_face_list.list();
+              //   const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
+              //   std::set<NodeId> node_set;
+              //   for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+              //     FaceId face_id         = face_list[i_face];
+              //     const auto& face_nodes = face_to_node_matrix[face_id];
+              //     for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+              //       node_set.insert(face_nodes[i_node]);
+              //     }
+              //   }
 
-                Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
-                  TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
+              //   Array<NodeId> node_list = convert_to_array(node_set);
 
-                std::shared_ptr bc =
-                  std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
-                bc_list.push_back(BoundaryConditionHandler(bc));
-              }
+              //   Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+              //     TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(),
+              //     node_list);
+
+              //   bc_list.push_back(VelocityBoundaryCondition{node_list, value_list});
+              // }
             }
           }
           break;
         }
         case IBoundaryConditionDescriptor::Type::pressure: {
+          using PressureBoundaryCondition = typename AcousticSolver<MeshType>::PressureBoundaryCondition;
+
           const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
             dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor);
           for (size_t i_ref_face_list = 0;
@@ -209,9 +218,7 @@ struct GlaceScheme
                 }
               }();
 
-              std::shared_ptr bc =
-                std::make_shared<PressureBoundaryCondition<MeshType::Dimension>>(face_list, face_values);
-              bc_list.push_back(BoundaryConditionHandler(bc));
+              bc_list.push_back(PressureBoundaryCondition{face_list, face_values});
             }
           }
           break;
@@ -283,12 +290,12 @@ struct GlaceScheme
     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);
+      // 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;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 1fcfa105e..88007ac47 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -1,47 +1,167 @@
 #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/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;
+
+  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)
@@ -121,100 +241,9 @@ class AcousticSolver
   void
   applyBoundaryConditions()
   {
-    for (const auto& handler : m_boundary_condition_list) {
-      switch (handler.boundaryCondition().type()) {
-      case BoundaryCondition::velocity: {
-        const VelocityBoundaryCondition<Dimension>& velocity_bc =
-          dynamic_cast<const VelocityBoundaryCondition<Dimension>&>(handler.boundaryCondition());
-
-        const auto& node_list  = velocity_bc.nodeList();
-        const auto& value_list = velocity_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;
-          });
-        break;
-      }
-      case BoundaryCondition::pressure: {
-        const PressureBoundaryCondition<Dimension>& pressure_bc =
-          dynamic_cast<const PressureBoundaryCondition<Dimension>&>(handler.boundaryCondition());
-
-        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  = pressure_bc.faceList();
-          const auto& value_list = pressure_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  = pressure_bc.faceList();
-          const auto& value_list = pressure_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);
-            }
-          }
-        }
-        break;
-      }
-      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>
@@ -291,7 +320,7 @@ 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)
     : m_mesh(p_mesh),
       m_connectivity(m_mesh->connectivity()),
       m_boundary_condition_list(bc_list),
@@ -388,4 +417,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/BoundaryCondition.hpp b/src/scheme/BoundaryCondition.hpp
deleted file mode 100644
index 2bd924974..000000000
--- a/src/scheme/BoundaryCondition.hpp
+++ /dev/null
@@ -1,187 +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,
-    pressure,
-    symmetry
-  };
-
-  const Type m_type;
-
- protected:
-  BoundaryCondition(Type type) : m_type(type)
-  {
-    ;
-  }
-
- public:
-  const Type&
-  type() const
-  {
-    return m_type;
-  }
-
-  virtual ~BoundaryCondition() = default;
-};
-
-template <size_t Dimension>
-class PressureBoundaryCondition final : public BoundaryCondition   // clazy:exclude=copyable-polymorphic
-{
- 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)
-    : BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, m_face_list{face_list}
-  {}
-
-  ~PressureBoundaryCondition() = default;
-};
-
-template <>
-class PressureBoundaryCondition<1> final : public BoundaryCondition   // clazy:exclude=copyable-polymorphic
-{
- 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)
-    : BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, 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;
-};
-
-template <size_t Dimension>
-class VelocityBoundaryCondition final : public BoundaryCondition   // clazy:exclude=copyable-polymorphic
-{
- 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)
-    : BoundaryCondition{BoundaryCondition::velocity}, m_value_list{value_list}, m_node_list{node_list}
-  {}
-
-  ~VelocityBoundaryCondition() = 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
-- 
GitLab