diff --git a/src/language/modules/CoreModule.cpp b/src/language/modules/CoreModule.cpp
index 0c4f2b354adb6a059fb92cff28e196178a08effb..1cc699455e23c123fbd1fee059f5296917a46b9b 100644
--- a/src/language/modules/CoreModule.cpp
+++ b/src/language/modules/CoreModule.cpp
@@ -1,5 +1,8 @@
 #include <language/modules/CoreModule.hpp>
 
+#include <language/modules/CoreModule.hpp>
+#include <language/modules/ModuleRepository.hpp>
+#include <language/utils/ASTExecutionInfo.hpp>
 #include <language/utils/AffectationProcessorBuilder.hpp>
 #include <language/utils/AffectationRegisterForB.hpp>
 #include <language/utils/AffectationRegisterForN.hpp>
@@ -8,7 +11,6 @@
 #include <language/utils/AffectationRegisterForRnxn.hpp>
 #include <language/utils/AffectationRegisterForString.hpp>
 #include <language/utils/AffectationRegisterForZ.hpp>
-
 #include <language/utils/BinaryOperatorRegisterForB.hpp>
 #include <language/utils/BinaryOperatorRegisterForN.hpp>
 #include <language/utils/BinaryOperatorRegisterForR.hpp>
@@ -16,23 +18,21 @@
 #include <language/utils/BinaryOperatorRegisterForRnxn.hpp>
 #include <language/utils/BinaryOperatorRegisterForString.hpp>
 #include <language/utils/BinaryOperatorRegisterForZ.hpp>
-
+#include <language/utils/BuiltinFunctionEmbedder.hpp>
 #include <language/utils/IncDecOperatorRegisterForN.hpp>
 #include <language/utils/IncDecOperatorRegisterForR.hpp>
 #include <language/utils/IncDecOperatorRegisterForZ.hpp>
-
 #include <language/utils/UnaryOperatorRegisterForB.hpp>
 #include <language/utils/UnaryOperatorRegisterForN.hpp>
 #include <language/utils/UnaryOperatorRegisterForR.hpp>
 #include <language/utils/UnaryOperatorRegisterForRn.hpp>
 #include <language/utils/UnaryOperatorRegisterForRnxn.hpp>
 #include <language/utils/UnaryOperatorRegisterForZ.hpp>
-
-#include <language/modules/CoreModule.hpp>
-#include <language/modules/ModuleRepository.hpp>
-#include <language/utils/ASTExecutionInfo.hpp>
-#include <language/utils/BuiltinFunctionEmbedder.hpp>
+#include <utils/Messenger.hpp>
 #include <utils/PugsUtils.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <random>
 
 CoreModule::CoreModule() : BuiltinModule(true)
 {
@@ -69,6 +69,20 @@ CoreModule::CoreModule() : BuiltinModule(true)
                                                }
 
                                                ));
+
+  this->_addBuiltinFunction("setRandomSeed", std::make_shared<BuiltinFunctionEmbedder<void(const int64_t&)>>(
+
+                                               [](const int64_t& random_seed) -> void {
+                                                 RandomEngine::instance().setRandomSeed(random_seed);
+                                               }
+
+                                               ));
+
+  this->_addBuiltinFunction("resetRandomSeed", std::make_shared<BuiltinFunctionEmbedder<void(void)>>(
+
+                                                 []() { RandomEngine::instance().resetRandomSeed(); }
+
+                                                 ));
 }
 
 void
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b4022d660475cf93eb9f9e1cb291622de22dba22..634ffe04fe88817042f028319ea1d22e8104d98c 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -10,7 +10,9 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshRandomizer.hpp>
 #include <scheme/AcousticSolver.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0Vector.hpp>
@@ -92,6 +94,84 @@ SchemeModule::SchemeModule()
 
       ));
 
+  this->_addBuiltinFunction("randomizeMesh",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IMesh>(std::shared_ptr<const IMesh>,
+                                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list) -> std::shared_ptr<const IMesh> {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  constexpr size_t Dimension = 1;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh();
+                                }
+                                case 2: {
+                                  constexpr size_t Dimension = 2;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh();
+                                }
+                                case 3: {
+                                  constexpr size_t Dimension = 3;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh();
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
+  this->_addBuiltinFunction("randomizeMesh",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IMesh>(std::shared_ptr<const IMesh>,
+                                           const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&,
+                                           const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IMesh> p_mesh,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const FunctionSymbolId& function_symbol_id) -> std::shared_ptr<const IMesh> {
+                                switch (p_mesh->dimension()) {
+                                case 1: {
+                                  constexpr size_t Dimension = 1;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh(function_symbol_id);
+                                }
+                                case 2: {
+                                  constexpr size_t Dimension = 2;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh(function_symbol_id);
+                                }
+                                case 3: {
+                                  constexpr size_t Dimension = 3;
+                                  using MeshType             = Mesh<Connectivity<Dimension>>;
+                                  const MeshType& mesh       = dynamic_cast<const MeshType&>(*p_mesh);
+                                  MeshRandomizer randomizer(mesh, bc_descriptor_list);
+                                  return randomizer.getRandomizedMesh(function_symbol_id);
+                                }
+                                default: {
+                                  throw UnexpectedError("invalid mesh dimension");
+                                }
+                                }
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("boundaryName",
                             std::make_shared<
                               BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>(const std::string&)>>(
@@ -112,6 +192,30 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this
+    ->_addBuiltinFunction("fixed",
+                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
+                            std::shared_ptr<const IBoundaryDescriptor>)>>(
+
+                            [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+                              -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                              return std::make_shared<FixedBoundaryConditionDescriptor>(boundary);
+                            }
+
+                            ));
+
+  this
+    ->_addBuiltinFunction("axis",
+                          std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
+                            std::shared_ptr<const IBoundaryDescriptor>)>>(
+
+                            [](std::shared_ptr<const IBoundaryDescriptor> boundary)
+                              -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                              return std::make_shared<AxisBoundaryConditionDescriptor>(boundary);
+                            }
+
+                            ));
+
   this
     ->_addBuiltinFunction("symmetry",
                           std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryConditionDescriptor>(
diff --git a/src/main.cpp b/src/main.cpp
index aad67d58a5d228cdc502e9b024cfc497fef23291..cf382d69b491e85968d977e9023ac094d81c3bc2 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -1,10 +1,10 @@
-#include <utils/PugsUtils.hpp>
-
 #include <language/PugsParser.hpp>
 #include <mesh/DiamondDualConnectivityManager.hpp>
 #include <mesh/DiamondDualMeshManager.hpp>
 #include <mesh/MeshDataManager.hpp>
 #include <mesh/SynchronizerManager.hpp>
+#include <utils/PugsUtils.hpp>
+#include <utils/RandomEngine.hpp>
 
 int
 main(int argc, char* argv[])
@@ -12,6 +12,7 @@ main(int argc, char* argv[])
   std::string filename = initialize(argc, argv);
 
   SynchronizerManager::create();
+  RandomEngine::create();
   MeshDataManager::create();
   DiamondDualConnectivityManager::create();
   DiamondDualMeshManager::create();
@@ -21,7 +22,10 @@ main(int argc, char* argv[])
   DiamondDualMeshManager::destroy();
   DiamondDualConnectivityManager::destroy();
   MeshDataManager::destroy();
+  RandomEngine::destroy();
   SynchronizerManager::destroy();
+
   finalize();
+
   return 0;
 }
diff --git a/src/mesh/ItemArrayUtils.hpp b/src/mesh/ItemArrayUtils.hpp
index 618dc8c8e061e55f9db3950b7cff71b3b281b111..6025bd0f3b17763be3f4b6646ab6730be0d6de18 100644
--- a/src/mesh/ItemArrayUtils.hpp
+++ b/src/mesh/ItemArrayUtils.hpp
@@ -23,4 +23,27 @@ synchronize(ItemArray<DataType, item_type, ConnectivityPtr>& item_array)
   }
 }
 
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+bool
+isSynchronized(ItemArray<const DataType, item_type, ConnectivityPtr> item_array)
+{
+  bool is_synchronized = true;
+
+  if (parallel::size() > 1) {
+    ItemArray<std::remove_const_t<DataType>, item_type> item_array_copy = copy(item_array);
+    synchronize(item_array_copy);
+    for (ItemIdT<item_type> item_id = 0; item_id < item_array_copy.numberOfItems(); ++item_id) {
+      for (size_t i = 0; i < item_array.sizeOfArrays(); ++i) {
+        if (item_array_copy[item_id][i] != item_array[item_id][i]) {
+          is_synchronized = false;
+        }
+      }
+    }
+
+    is_synchronized = parallel::allReduceAnd(is_synchronized);
+  }
+
+  return is_synchronized;
+}
+
 #endif   // ITEM_ARRAY_UTILS_HPP
diff --git a/src/mesh/ItemValueUtils.hpp b/src/mesh/ItemValueUtils.hpp
index ab39d5c5c6423dac3d8c7d3fd2c1cbeffd1a9809..b5151953c8bff1501124b027ef82c259c7bc17b7 100644
--- a/src/mesh/ItemValueUtils.hpp
+++ b/src/mesh/ItemValueUtils.hpp
@@ -304,4 +304,26 @@ synchronize(ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
   }
 }
 
+template <typename DataType, ItemType item_type, typename ConnectivityPtr>
+bool
+isSynchronized(const ItemValue<DataType, item_type, ConnectivityPtr>& item_value)
+{
+  bool is_synchronized = true;
+
+  if (parallel::size() > 1) {
+    ItemValue<std::remove_const_t<DataType>, item_type> item_value_copy = copy(item_value);
+    synchronize(item_value_copy);
+    for (ItemIdT<item_type> item_id = 0; item_id < item_value_copy.numberOfItems(); ++item_id) {
+      if (item_value_copy[item_id] != item_value[item_id]) {
+        is_synchronized = false;
+        break;
+      }
+    }
+
+    is_synchronized = parallel::allReduceAnd(is_synchronized);
+  }
+
+  return is_synchronized;
+}
+
 #endif   // ITEM_VALUE_UTILS_HPP
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 92d354fc52ddfd3e35ab69a8d892edb32f6432de..0743d2d95dda05ca0cb990ade1c8d088c9ed2112 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -1,28 +1,26 @@
 #ifndef MESH_NODE_BOUNDARY_HPP
 #define MESH_NODE_BOUNDARY_HPP
 
-#include <utils/Array.hpp>
-
+#include <Kokkos_Vector.hpp>
+#include <algebra/TinyMatrix.hpp>
 #include <algebra/TinyVector.hpp>
-
-#include <mesh/ItemValue.hpp>
-#include <mesh/RefItemList.hpp>
-
 #include <mesh/ConnectivityMatrix.hpp>
 #include <mesh/IConnectivity.hpp>
-
+#include <mesh/ItemValue.hpp>
+#include <mesh/RefItemList.hpp>
+#include <utils/Array.hpp>
 #include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
 
-#include <Kokkos_Vector.hpp>
-
-#include <iostream>
-
 template <size_t Dimension>
 class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
 {
  protected:
   Array<const NodeId> m_node_list;
+  std::string m_boundary_name;
+
+  template <typename MeshType>
+  std::array<TinyVector<Dimension>, Dimension*(Dimension - 1)> _getBounds(const MeshType& mesh) const;
 
  public:
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
@@ -35,24 +33,33 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
+  MeshNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : m_boundary_name(ref_face_list.refId().tagName())
   {
     static_assert(Dimension == MeshType::Dimension);
-    const auto& face_to_cell_matrix = mesh->connectivity().faceToCellMatrix();
+    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
 
     const Array<const FaceId>& face_list = ref_face_list.list();
-    parallel_for(
-      face_list.size(), PUGS_LAMBDA(int l) {
-        const auto& face_cells = face_to_cell_matrix[face_list[l]];
-        if (face_cells.size() > 1) {
-          throw NormalError("internal faces cannot be used to define mesh boundaries");
-        }
-      });
+
+    bool is_bad = false;
+    parallel_for(face_list.size(), [=, &is_bad](int l) {
+      const auto& face_cells = face_to_cell_matrix[face_list[l]];
+      if (face_cells.size() > 1) {
+        is_bad = true;
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": inner faces cannot be used to define mesh boundaries";
+      throw NormalError(ost.str());
+    }
 
     Kokkos::vector<unsigned int> node_ids;
     // not enough but should reduce significantly the number of resizing
     node_ids.reserve(Dimension * face_list.size());
-    const auto& face_to_node_matrix = mesh->connectivity().faceToNodeMatrix();
+    const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
 
     for (size_t l = 0; l < face_list.size(); ++l) {
       const FaceId face_number = face_list[l];
@@ -73,20 +80,245 @@ class MeshNodeBoundary   // clazy:exclude=copyable-polymorphic
   }
 
   template <typename MeshType>
-  MeshNodeBoundary(const std::shared_ptr<const MeshType>&, const RefNodeList& ref_node_list)
-    : m_node_list(ref_node_list.list())
+  MeshNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
+    : m_boundary_name(ref_edge_list.refId().tagName())
   {
     static_assert(Dimension == MeshType::Dimension);
+    const auto& edge_to_face_matrix = mesh.connectivity().edgeToFaceMatrix();
+    const auto& face_to_cell_matrix = mesh.connectivity().faceToCellMatrix();
+
+    auto edge_is_owned = mesh.connectivity().edgeIsOwned();
+
+    const Array<const EdgeId>& edge_list = ref_edge_list.list();
+
+    bool is_bad = false;
+    parallel_for(edge_list.size(), [=, &is_bad](int l) {
+      const EdgeId edge_id = edge_list[l];
+      if (edge_is_owned[edge_id]) {
+        const auto& edge_faces             = edge_to_face_matrix[edge_id];
+        bool is_connected_to_boundary_face = false;
+        for (size_t i_edge_face = 0; i_edge_face < edge_faces.size(); ++i_edge_face) {
+          const FaceId edge_face_id = edge_faces[i_edge_face];
+          if (face_to_cell_matrix[edge_face_id].size() == 1) {
+            is_connected_to_boundary_face = true;
+          }
+        }
+        if (not is_connected_to_boundary_face) {
+          is_bad = true;
+        }
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": inner edges cannot be used to define mesh boundaries";
+      throw NormalError(ost.str());
+    }
+
+    Kokkos::vector<unsigned int> node_ids;
+    node_ids.reserve(2 * edge_list.size());
+    const auto& edge_to_node_matrix = mesh.connectivity().edgeToNodeMatrix();
+
+    for (size_t l = 0; l < edge_list.size(); ++l) {
+      const EdgeId edge_number = edge_list[l];
+      const auto& edge_nodes   = edge_to_node_matrix[edge_number];
+
+      for (size_t r = 0; r < edge_nodes.size(); ++r) {
+        node_ids.push_back(edge_nodes[r]);
+      }
+    }
+    std::sort(node_ids.begin(), node_ids.end());
+    auto last = std::unique(node_ids.begin(), node_ids.end());
+    node_ids.resize(std::distance(node_ids.begin(), last));
+
+    Array<NodeId> node_list(node_ids.size());
+    parallel_for(
+      node_ids.size(), PUGS_LAMBDA(int r) { node_list[r] = node_ids[r]; });
+    m_node_list = node_list;
   }
 
-  MeshNodeBoundary()          = default;
-  virtual ~MeshNodeBoundary() = default;
+  template <typename MeshType>
+  MeshNodeBoundary(const MeshType&, const RefNodeList& ref_node_list)
+    : m_node_list(ref_node_list.list()), m_boundary_name(ref_node_list.refId().tagName())
+
+  {
+    static_assert(Dimension == MeshType::Dimension);
+  }
 
- protected:
   MeshNodeBoundary(const MeshNodeBoundary&) = default;
   MeshNodeBoundary(MeshNodeBoundary&&)      = default;
+
+  MeshNodeBoundary()          = default;
+  virtual ~MeshNodeBoundary() = default;
 };
 
+template <>
+template <typename MeshType>
+std::array<TinyVector<2>, 2>
+MeshNodeBoundary<2>::_getBounds(const MeshType& mesh) const
+{
+  static_assert(MeshType::Dimension == 2);
+
+  using R2 = TinyVector<2, double>;
+
+  const NodeValue<const R2>& xr = mesh.xr();
+
+  std::array<R2, 2> bounds;
+  R2& xmin = bounds[0];
+  R2& xmax = bounds[1];
+
+  xmin = R2{std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  xmax = R2{-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max()};
+
+  auto update_xmin = [](const R2& x, R2& xmin) {
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R2& x, R2& xmax) {
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
+      xmax = x;
+    }
+  };
+
+  for (size_t r = 0; r < m_node_list.size(); ++r) {
+    const R2& x = xr[m_node_list[r]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+  }
+
+  if (parallel::size() > 1) {
+    Array<R2> xmin_array = parallel::allGather(xmin);
+    Array<R2> xmax_array = parallel::allGather(xmax);
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+  }
+
+  return bounds;
+}
+
+template <>
+template <typename MeshType>
+std::array<TinyVector<3>, 6>
+MeshNodeBoundary<3>::_getBounds(const MeshType& mesh) const
+{
+  static_assert(MeshType::Dimension == 3);
+
+  using R3 = TinyVector<3, double>;
+
+  auto update_xmin = [](const R3& x, R3& xmin) {
+    // XMIN: X.xmin X.ymax X.zmax
+    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] > xmin[1])) or
+        ((x[0] == xmin[0]) and (x[1] == xmin[1]) and (x[2] > xmin[2]))) {
+      xmin = x;
+    }
+  };
+
+  auto update_xmax = [](const R3& x, R3& xmax) {
+    // XMAX: X.xmax X.ymin X.zmin
+    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] < xmax[1])) or
+        ((x[0] == xmax[0]) and (x[1] == xmax[1]) and (x[2] < xmax[2]))) {
+      xmax = x;
+    }
+  };
+
+  auto update_ymin = [](const R3& x, R3& ymin) {
+    // YMIN: X.ymin X.zmax X.xmin
+    if ((x[1] < ymin[1]) or ((x[1] == ymin[1]) and (x[2] > ymin[2])) or
+        ((x[1] == ymin[1]) and (x[2] == ymin[2]) and (x[0] < ymin[0]))) {
+      ymin = x;
+    }
+  };
+
+  auto update_ymax = [](const R3& x, R3& ymax) {
+    // YMAX: X.ymax X.zmin X.xmax
+    if ((x[1] > ymax[1]) or ((x[1] == ymax[1]) and (x[2] < ymax[2])) or
+        ((x[1] == ymax[1]) and (x[2] == ymax[1]) and (x[0] > ymax[0]))) {
+      ymax = x;
+    }
+  };
+
+  auto update_zmin = [](const R3& x, R3& zmin) {
+    // ZMIN: X.zmin X.xmin X.ymin
+    if ((x[2] < zmin[2]) or ((x[2] == zmin[2]) and (x[2] < zmin[2])) or
+        ((x[1] == zmin[1]) and (x[2] == zmin[1]) and (x[0] < zmin[0]))) {
+      zmin = x;
+    }
+  };
+
+  auto update_zmax = [](const R3& x, R3& zmax) {
+    // ZMAX: X.zmax X.xmax X.ymax
+    if ((x[2] > zmax[2]) or ((x[2] == zmax[2]) and (x[2] > zmax[2])) or
+        ((x[1] == zmax[1]) and (x[2] == zmax[1]) and (x[0] > zmax[0]))) {
+      zmax = x;
+    }
+  };
+
+  std::array<R3, 6> bounds;
+  R3& xmin = bounds[0];
+  R3& ymin = bounds[1];
+  R3& zmin = bounds[2];
+  R3& xmax = bounds[3];
+  R3& ymax = bounds[4];
+  R3& zmax = bounds[5];
+
+  xmin = R3{std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max()};
+  ymin = xmin;
+  zmin = xmin;
+
+  xmax = -xmin;
+  ymax = xmax;
+  zmax = xmax;
+
+  const NodeValue<const R3>& xr = mesh.xr();
+
+  for (size_t r = 0; r < m_node_list.size(); ++r) {
+    const R3& x = xr[m_node_list[r]];
+    update_xmin(x, xmin);
+    update_xmax(x, xmax);
+    update_ymin(x, ymin);
+    update_ymax(x, ymax);
+    update_zmin(x, zmin);
+    update_zmax(x, zmax);
+  }
+
+  if (parallel::size() > 0) {
+    Array<const R3> xmin_array = parallel::allGather(xmin);
+    Array<const R3> xmax_array = parallel::allGather(xmax);
+    Array<const R3> ymin_array = parallel::allGather(ymin);
+    Array<const R3> ymax_array = parallel::allGather(ymax);
+    Array<const R3> zmin_array = parallel::allGather(zmin);
+    Array<const R3> zmax_array = parallel::allGather(zmax);
+
+    for (size_t i = 0; i < xmin_array.size(); ++i) {
+      update_xmin(xmin_array[i], xmin);
+    }
+    for (size_t i = 0; i < ymin_array.size(); ++i) {
+      update_ymin(ymin_array[i], ymin);
+    }
+    for (size_t i = 0; i < zmin_array.size(); ++i) {
+      update_zmin(zmin_array[i], zmin);
+    }
+    for (size_t i = 0; i < xmax_array.size(); ++i) {
+      update_xmax(xmax_array[i], xmax);
+    }
+    for (size_t i = 0; i < ymax_array.size(); ++i) {
+      update_ymax(ymax_array[i], ymax);
+    }
+    for (size_t i = 0; i < zmax_array.size(); ++i) {
+      update_zmax(zmax_array[i], zmax);
+    }
+  }
+
+  return bounds;
+}
+
 template <size_t Dimension>
 class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
 {
@@ -97,16 +329,38 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   const Rd m_outgoing_normal;
 
   template <typename MeshType>
-  PUGS_INLINE Rd _getNormal(const std::shared_ptr<const MeshType>& mesh);
+  PUGS_INLINE Rd _getNormal(const MeshType& mesh);
 
   template <typename MeshType>
-  PUGS_INLINE void _checkBoundaryIsFlat(TinyVector<2, double> normal,
-                                        TinyVector<2, double> xmin,
-                                        TinyVector<2, double> xmax,
-                                        const std::shared_ptr<const MeshType>& mesh) const;
+  PUGS_INLINE void
+  _checkBoundaryIsFlat(const TinyVector<Dimension, double>& normal,
+                       const TinyVector<Dimension, double>& origin,
+                       const double length,
+                       const MeshType& mesh) const
+  {
+    static_assert(MeshType::Dimension == Dimension);
+
+    const NodeValue<const Rd>& xr = mesh.xr();
+
+    bool is_bad = false;
+
+    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+      const Rd& x = xr[this->m_node_list[r]];
+      if (dot(x - origin, normal) > 1E-13 * length) {
+        is_bad = true;
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": boundary is not flat!";
+      throw NormalError(ost.str());
+    }
+  }
 
   template <typename MeshType>
-  PUGS_INLINE Rd _getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh);
+  PUGS_INLINE Rd _getOutgoingNormal(const MeshType& mesh);
 
  public:
   const Rd&
@@ -119,14 +373,14 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
   template <typename MeshType>
-  MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefFaceList& ref_face_list)
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
     : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
   }
 
   template <typename MeshType>
-  MeshFlatNodeBoundary(const std::shared_ptr<const MeshType>& mesh, const RefNodeList& ref_node_list)
+  MeshFlatNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
     : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
@@ -138,42 +392,17 @@ class MeshFlatNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclu
   virtual ~MeshFlatNodeBoundary()                   = default;
 };
 
-template <>
-template <typename MeshType>
-void
-MeshFlatNodeBoundary<2>::_checkBoundaryIsFlat(TinyVector<2, double> normal,
-                                              TinyVector<2, double> xmin,
-                                              TinyVector<2, double> xmax,
-                                              const std::shared_ptr<const MeshType>& mesh) const
-{
-  static_assert(MeshType::Dimension == 2);
-  using R2 = TinyVector<2, double>;
-
-  const R2 origin     = 0.5 * (xmin + xmax);
-  const double length = l2Norm(xmax - xmin);
-
-  const NodeValue<const R2>& xr = mesh->xr();
-
-  parallel_for(
-    m_node_list.size(), PUGS_LAMBDA(size_t r) {
-      const R2& x = xr[m_node_list[r]];
-      if (dot(x - origin, normal) > 1E-13 * length) {
-        throw NormalError("this FlatBoundary is not flat!");
-      }
-    });
-}
-
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<1>::_getNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 1);
   using R = TinyVector<1, double>;
 
   const size_t number_of_bc_nodes = [&]() {
     size_t number_of_bc_nodes = 0;
-    auto node_is_owned        = mesh->connectivity().nodeIsOwned();
+    auto node_is_owned        = mesh.connectivity().nodeIsOwned();
     for (size_t i_node = 0; i_node < m_node_list.size(); ++i_node) {
       number_of_bc_nodes += (node_is_owned[m_node_list[i_node]]);
     }
@@ -181,7 +410,10 @@ MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   }();
 
   if (number_of_bc_nodes != 1) {
-    throw NormalError("Node boundaries in 1D require to have exactly 1 node");
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << m_boundary_name << rang::style::reset
+        << ": node boundaries in 1D require to have exactly 1 node";
+    throw NormalError(ost.str());
   }
 
   return R{1};
@@ -190,54 +422,29 @@ MeshFlatNodeBoundary<1>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<2>::_getNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 2);
   using R2 = TinyVector<2, double>;
 
-  const NodeValue<const R2>& xr = mesh->xr();
-
-  R2 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
 
-  R2 xmax(-std::numeric_limits<double>::max(), -std::numeric_limits<double>::max());
-
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R2& x = xr[m_node_list[r]];
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
-      xmin = x;
-    }
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
-      xmax = x;
-    }
-  }
-
-  Array<R2> xmin_array = parallel::allGather(xmin);
-  Array<R2> xmax_array = parallel::allGather(xmax);
-  for (size_t i = 0; i < xmin_array.size(); ++i) {
-    const R2& x = xmin_array[i];
-    if ((x[0] < xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
-      xmin = x;
-    }
-  }
-  for (size_t i = 0; i < xmax_array.size(); ++i) {
-    const R2& x = xmax_array[i];
-    if ((x[0] > xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
-      xmax = x;
-    }
-  }
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
 
   if (xmin == xmax) {
-    std::stringstream os;
-    os << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
-    throw NormalError(os.str());
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
   }
 
   R2 dx = xmax - xmin;
   dx *= 1. / l2Norm(dx);
 
-  R2 normal(-dx[1], dx[0]);
+  const R2 normal(-dx[1], dx[0]);
 
-  this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+  this->_checkBoundaryIsFlat(normal, 0.5 * (xmin + xmax), l2Norm(xmax - xmin), mesh);
 
   return normal;
 }
@@ -245,85 +452,20 @@ MeshFlatNodeBoundary<2>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<3>::_getNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 3);
-  using R3 = TinyVector<3, double>;
 
-  R3 xmin(std::numeric_limits<double>::max(), std::numeric_limits<double>::max(), std::numeric_limits<double>::max());
-  R3 ymin = xmin;
-  R3 zmin = xmin;
+  using R3 = TinyVector<3, double>;
 
-  R3 xmax = -xmin;
-  R3 ymax = xmax;
-  R3 zmax = xmax;
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
 
-  const NodeValue<const R3>& xr = mesh->xr();
-
-  for (size_t r = 0; r < m_node_list.size(); ++r) {
-    const R3& x = xr[m_node_list[r]];
-    if (x[0] < xmin[0]) {
-      xmin = x;
-    }
-    if (x[1] < ymin[1]) {
-      ymin = x;
-    }
-    if (x[2] < zmin[2]) {
-      zmin = x;
-    }
-    if (x[0] > xmax[0]) {
-      xmax = x;
-    }
-    if (x[1] > ymax[1]) {
-      ymax = x;
-    }
-    if (x[2] > zmax[2]) {
-      zmax = x;
-    }
-  }
-  Array<R3> xmin_array = parallel::allGather(xmin);
-  Array<R3> xmax_array = parallel::allGather(xmax);
-  Array<R3> ymin_array = parallel::allGather(ymin);
-  Array<R3> ymax_array = parallel::allGather(ymax);
-  Array<R3> zmin_array = parallel::allGather(zmin);
-  Array<R3> zmax_array = parallel::allGather(zmax);
-
-  for (size_t i = 0; i < xmin_array.size(); ++i) {
-    const R3& x = xmin_array[i];
-    if (x[0] < xmin[0]) {
-      xmin = x;
-    }
-  }
-  for (size_t i = 0; i < ymin_array.size(); ++i) {
-    const R3& x = ymin_array[i];
-    if (x[1] < ymin[1]) {
-      ymin = x;
-    }
-  }
-  for (size_t i = 0; i < zmin_array.size(); ++i) {
-    const R3& x = zmin_array[i];
-    if (x[2] < zmin[2]) {
-      zmin = x;
-    }
-  }
-  for (size_t i = 0; i < xmax_array.size(); ++i) {
-    const R3& x = xmax_array[i];
-    if (x[0] > xmax[0]) {
-      xmax = x;
-    }
-  }
-  for (size_t i = 0; i < ymax_array.size(); ++i) {
-    const R3& x = ymax_array[i];
-    if (x[1] > ymax[1]) {
-      ymax = x;
-    }
-  }
-  for (size_t i = 0; i < zmax_array.size(); ++i) {
-    const R3& x = zmax_array[i];
-    if (x[2] > zmax[2]) {
-      zmax = x;
-    }
-  }
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
 
   const R3 u = xmax - xmin;
   const R3 v = ymax - ymin;
@@ -352,12 +494,17 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
   }
 
   if (normal_l2 == 0) {
-    throw NormalError("cannot to compute normal!");
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute normal";
+    throw NormalError(ost.str());
   }
 
-  normal *= 1. / sqrt(normal_l2);
+  const double length = sqrt(normal_l2);
+
+  normal *= 1. / length;
 
-  // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+  this->_checkBoundaryIsFlat(normal, 1. / 6. * (xmin + xmax + ymin + ymax + zmin + zmax), length, mesh);
 
   return normal;
 }
@@ -365,7 +512,7 @@ MeshFlatNodeBoundary<3>::_getNormal(const std::shared_ptr<const MeshType>& mesh)
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<1, double>
-MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<1>::_getOutgoingNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 1);
   using R = TinyVector<1, double>;
@@ -375,10 +522,10 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R>& xr    = mesh->xr();
-    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+    const NodeValue<const R>& xr    = mesh.xr();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
@@ -410,7 +557,7 @@ MeshFlatNodeBoundary<1>::_getOutgoingNormal(const std::shared_ptr<const MeshType
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<2, double>
-MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<2>::_getOutgoingNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 2);
   using R2 = TinyVector<2, double>;
@@ -420,10 +567,10 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R2>& xr   = mesh->xr();
-    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+    const NodeValue<const R2>& xr   = mesh.xr();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
@@ -454,7 +601,7 @@ MeshFlatNodeBoundary<2>::_getOutgoingNormal(const std::shared_ptr<const MeshType
 template <>
 template <typename MeshType>
 PUGS_INLINE TinyVector<3, double>
-MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType>& mesh)
+MeshFlatNodeBoundary<3>::_getOutgoingNormal(const MeshType& mesh)
 {
   static_assert(MeshType::Dimension == 3);
   using R3 = TinyVector<3, double>;
@@ -464,10 +611,10 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType
   double max_height = 0;
 
   if (m_node_list.size() > 0) {
-    const NodeValue<const R3>& xr   = mesh->xr();
-    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+    const NodeValue<const R3>& xr   = mesh.xr();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
-    const auto& node_to_cell_matrix = mesh->connectivity().nodeToCellMatrix();
+    const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
 
     const NodeId r0      = m_node_list[0];
     const CellId j0      = node_to_cell_matrix[r0][0];
@@ -496,4 +643,159 @@ MeshFlatNodeBoundary<3>::_getOutgoingNormal(const std::shared_ptr<const MeshType
   }
 }
 
+template <size_t Dimension>
+class MeshLineNodeBoundary : public MeshNodeBoundary<Dimension>   // clazy:exclude=copyable-polymorphic
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const Rd m_direction;
+
+  template <typename MeshType>
+  PUGS_INLINE TinyVector<Dimension> _getDirection(const MeshType&);
+
+  template <typename MeshType>
+  PUGS_INLINE void
+  _checkBoundaryIsLine(const TinyVector<Dimension, double>& direction,
+                       const TinyVector<Dimension, double>& origin,
+                       const double length,
+                       const MeshType& mesh) const
+  {
+    static_assert(MeshType::Dimension == Dimension);
+    using Rdxd = TinyMatrix<Dimension>;
+
+    const NodeValue<const Rd>& xr = mesh.xr();
+
+    Rdxd P = Rdxd{identity} - tensorProduct(direction, direction);
+
+    bool is_bad = false;
+    parallel_for(this->m_node_list.size(), [=, &is_bad](int r) {
+      const Rd& x    = xr[this->m_node_list[r]];
+      const Rd delta = x - origin;
+      if (dot(P * delta, direction) > 1E-13 * length) {
+        is_bad = true;
+      }
+    });
+
+    if (parallel::allReduceOr(is_bad)) {
+      std::ostringstream ost;
+      ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+          << ": boundary is not a line!";
+      throw NormalError(ost.str());
+    }
+  }
+
+ public:
+  const Rd&
+  direction() const
+  {
+    return m_direction;
+  }
+
+  MeshLineNodeBoundary& operator=(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary& operator=(MeshLineNodeBoundary&&) = default;
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefEdgeList& ref_edge_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_edge_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefFaceList& ref_face_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_face_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  template <typename MeshType>
+  MeshLineNodeBoundary(const MeshType& mesh, const RefNodeList& ref_node_list)
+    : MeshNodeBoundary<Dimension>(mesh, ref_node_list), m_direction(_getDirection(mesh))
+  {
+    ;
+  }
+
+  MeshLineNodeBoundary()                            = default;
+  MeshLineNodeBoundary(const MeshLineNodeBoundary&) = default;
+  MeshLineNodeBoundary(MeshLineNodeBoundary&&)      = default;
+  virtual ~MeshLineNodeBoundary()                   = default;
+};
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<1>
+MeshLineNodeBoundary<1>::_getDirection(const MeshType&)
+{
+  throw UnexpectedError("MeshLineNodeBoundary::_getDirection make no sense in dimension 1");
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<2>
+MeshLineNodeBoundary<2>::_getDirection(const MeshType& mesh)
+{
+  using R2 = TinyVector<2, double>;
+
+  std::array<R2, 2> bounds = this->_getBounds(mesh);
+
+  const R2& xmin = bounds[0];
+  const R2& xmax = bounds[1];
+
+  if (xmin == xmax) {
+    std::ostringstream ost;
+    ost << "invalid boundary " << rang::fgB::yellow << this->m_boundary_name << rang::style::reset
+        << ": unable to compute direction";
+    throw NormalError(ost.str());
+  }
+
+  R2 direction        = xmax - xmin;
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
+template <>
+template <typename MeshType>
+PUGS_INLINE TinyVector<3>
+MeshLineNodeBoundary<3>::_getDirection(const MeshType& mesh)
+{
+  static_assert(MeshType::Dimension == 3);
+
+  using R3 = TinyVector<3, double>;
+
+  std::array<R3, 6> bounds = this->_getBounds(mesh);
+
+  const R3& xmin = bounds[0];
+  const R3& ymin = bounds[1];
+  const R3& zmin = bounds[2];
+  const R3& xmax = bounds[3];
+  const R3& ymax = bounds[4];
+  const R3& zmax = bounds[5];
+
+  const R3 u = xmax - xmin;
+  const R3 v = ymax - ymin;
+  const R3 w = zmax - zmin;
+
+  R3 direction = u;
+  if (dot(v, v) > dot(direction, direction)) {
+    direction = v;
+  }
+
+  if (dot(w, w) > dot(direction, direction)) {
+    direction = w;
+  }
+
+  const double length = l2Norm(direction);
+  direction *= 1. / length;
+
+  this->_checkBoundaryIsLine(direction, xmin, length, mesh);
+
+  return direction;
+}
+
 #endif   // MESH_NODE_BOUNDARY_HPP
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1daaff6a04dfbfd4856c8f2d0875bb4cc4549ca5
--- /dev/null
+++ b/src/mesh/MeshRandomizer.hpp
@@ -0,0 +1,403 @@
+#ifndef MESH_RANDOMIZER_HPP
+#define MESH_RANDOMIZER_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <language/utils/InterpolateItemValue.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <variant>
+#include <vector>
+
+template <size_t Dimension>
+class MeshRandomizer
+{
+ private:
+  using Rd               = TinyVector<Dimension>;
+  using Rdxd             = TinyMatrix<Dimension>;
+  using ConnectivityType = Connectivity<Dimension>;
+  using MeshType         = Mesh<ConnectivityType>;
+
+  const MeshType& m_given_mesh;
+
+  class AxisBoundaryCondition;
+  class FixedBoundaryCondition;
+  class SymmetryBoundaryCondition;
+
+  using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
+
+  using BoundaryConditionList = std::vector<BoundaryCondition>;
+  BoundaryConditionList m_boundary_condition_list;
+
+  BoundaryConditionList
+  _getBCList(const MeshType& mesh,
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+  {
+    BoundaryConditionList bc_list;
+
+    constexpr ItemType EdgeType = [] {
+      if constexpr (Dimension == 3) {
+        return ItemType::edge;
+      } else if constexpr (Dimension == 2) {
+        return ItemType::face;
+      } else {
+        return ItemType::node;
+      }
+    }();
+
+    constexpr ItemType FaceType = [] {
+      if constexpr (Dimension > 1) {
+        return ItemType::face;
+      } else {
+        return ItemType::node;
+      }
+    }();
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::axis: {
+        const AxisBoundaryConditionDescriptor& axis_bc_descriptor =
+          dynamic_cast<const AxisBoundaryConditionDescriptor&>(*bc_descriptor);
+        for (size_t i_ref_edge_list = 0; i_ref_edge_list < mesh.connectivity().template numberOfRefItemList<EdgeType>();
+             ++i_ref_edge_list) {
+          const auto& ref_edge_list = mesh.connectivity().template refItemList<EdgeType>(i_ref_edge_list);
+          const RefId& ref          = ref_edge_list.refId();
+          if (ref == axis_bc_descriptor.boundaryDescriptor()) {
+            if constexpr (Dimension == 1) {
+              bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_edge_list}});
+            } else {
+              bc_list.emplace_back(AxisBoundaryCondition{MeshLineNodeBoundary<Dimension>(mesh, ref_edge_list)});
+            }
+          }
+        }
+        is_valid_boundary_condition = true;
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
+          dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
+             ++i_ref_face_list) {
+          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
+          const RefId& ref          = ref_face_list.refId();
+          if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+            bc_list.emplace_back(SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)});
+          }
+        }
+        is_valid_boundary_condition = true;
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fixed: {
+        const FixedBoundaryConditionDescriptor& fixed_bc_descriptor =
+          dynamic_cast<const FixedBoundaryConditionDescriptor&>(*bc_descriptor);
+        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
+             ++i_ref_face_list) {
+          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
+          const RefId& ref          = ref_face_list.refId();
+          if (ref == fixed_bc_descriptor.boundaryDescriptor()) {
+            bc_list.emplace_back(FixedBoundaryCondition{MeshNodeBoundary<Dimension>{mesh, ref_face_list}});
+          }
+        }
+        break;
+      }
+      default: {
+        is_valid_boundary_condition = false;
+      }
+      }
+      if (not is_valid_boundary_condition) {
+        std::ostringstream error_msg;
+        error_msg << *bc_descriptor << " is an invalid boundary condition for mesh randomizer";
+        throw NormalError(error_msg.str());
+      }
+    }
+
+    return bc_list;
+  }
+
+  void
+  _applyBC(NodeValue<Rd>& shift) const
+  {
+    for (auto&& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+            const Rd& n = bc.outgoingNormal();
+
+            const Rdxd I   = identity;
+            const Rdxd nxn = tensorProduct(n, n);
+            const Rdxd P   = I - nxn;
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                shift[node_id] = P * shift[node_id];
+              });
+
+          } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
+            const Rd& t = bc.direction();
+
+            const Rdxd txt = tensorProduct(t, t);
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+
+                shift[node_id] = txt * shift[node_id];
+              });
+
+          } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(const size_t i_node) {
+                const NodeId node_id = node_list[i_node];
+                shift[node_id]       = zero;
+              });
+
+          } else {
+            throw UnexpectedError("invalid boundary condition type");
+          }
+        },
+        boundary_condition);
+    }
+  }
+
+  NodeValue<Rd>
+  _getDisplacement() const
+  {
+    const ConnectivityType& connectivity = m_given_mesh.connectivity();
+    NodeValue<const Rd> given_xr         = m_given_mesh.xr();
+
+    auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
+    auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
+    auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+    NodeValue<double> max_delta_xr{connectivity};
+    parallel_for(
+      connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        const Rd& x0 = given_xr[node_id];
+
+        const auto& node_cell_list = node_to_cell_matrix[node_id];
+        double min_distance_2      = std::numeric_limits<double>::max();
+
+        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
+
+          const CellId cell_id       = node_cell_list[i_cell];
+          const auto& cell_node_list = cell_to_node_matrix[cell_id];
+
+          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+            if (i_node != i_cell_node) {
+              const NodeId cell_node_id = cell_node_list[i_node];
+              const Rd delta            = x0 - given_xr[cell_node_id];
+              min_distance_2            = std::min(min_distance_2, dot(delta, delta));
+            }
+          }
+        }
+        double max_delta = std::sqrt(min_distance_2);
+
+        max_delta_xr[node_id] = max_delta;
+      });
+
+    synchronize(max_delta_xr);
+
+    std::uniform_real_distribution<> dis(-0.5, 0.5);
+
+    NodeValue<const int> node_numbers = connectivity.nodeNumber();
+    using IdCorrespondance            = std::pair<int, NodeId>;
+    Array<IdCorrespondance> node_numbers_to_node_id{node_numbers.numberOfItems()};
+    parallel_for(
+      node_numbers.numberOfItems(), PUGS_LAMBDA(const NodeId node_id) {
+        node_numbers_to_node_id[node_id] = std::make_pair(node_numbers[node_id], node_id);
+      });
+
+    std::sort(&node_numbers_to_node_id[0], &node_numbers_to_node_id[0] + node_numbers_to_node_id.size(),
+              [](IdCorrespondance a, IdCorrespondance b) { return a.first < b.first; });
+
+    RandomEngine& random_engine = RandomEngine::instance();
+
+    Assert(isSynchronized(random_engine), "seed is not synchronized when entering mesh randomization");
+
+    NodeValue<Rd> shift_r{connectivity};
+
+    int i_node_number = 0;
+    for (size_t i = 0; i < node_numbers_to_node_id.size(); ++i) {
+      const auto [node_number, node_id] = node_numbers_to_node_id[i];
+      while (i_node_number < node_number) {
+        for (size_t j = 0; j < Dimension; ++j) {
+          dis(random_engine.engine());
+        }
+        ++i_node_number;
+      }
+
+      double max_delta = max_delta_xr[node_id];
+
+      Rd shift;
+      for (size_t i_component = 0; i_component < Dimension; ++i_component) {
+        shift[i_component] = max_delta * dis(random_engine.engine());
+      }
+
+      shift_r[node_id] = shift;
+
+      ++i_node_number;
+    }
+
+    const int max_node_number =
+      parallel::allReduceMax(node_numbers_to_node_id[node_numbers_to_node_id.size() - 1].first);
+
+    // Advances random engine to preserve CPU random number generators synchronization
+    for (; i_node_number <= max_node_number; ++i_node_number) {
+      for (size_t j = 0; j < Dimension; ++j) {
+        dis(random_engine.engine());
+      }
+    }
+
+    this->_applyBC(shift_r);
+
+#ifndef NDEBUG
+    if (not isSynchronized(shift_r)) {
+      throw UnexpectedError("randomized mesh coordinates are not synchronized");
+    }
+#endif   // NDEBUG
+
+    Assert(isSynchronized(random_engine), "seed is not synchronized after mesh randomization");
+
+    return shift_r;
+  }
+
+ public:
+  std::shared_ptr<const MeshType>
+  getRandomizedMesh() const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { xr[node_id] += given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  std::shared_ptr<const MeshType>
+  getRandomizedMesh(const FunctionSymbolId& function_symbol_id) const
+  {
+    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    NodeValue<const bool> is_displaced =
+      InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
+
+    NodeValue<Rd> xr = this->_getDisplacement();
+
+    parallel_for(
+      m_given_mesh.numberOfNodes(),
+      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+  }
+
+  MeshRandomizer(const MeshRandomizer&) = delete;
+  MeshRandomizer(MeshRandomizer&&)      = delete;
+
+  MeshRandomizer(const MeshType& given_mesh,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+    : m_given_mesh(given_mesh), m_boundary_condition_list(this->_getBCList(given_mesh, bc_descriptor_list))
+  {}
+
+  ~MeshRandomizer() = default;
+};
+
+template <size_t Dimension>
+class MeshRandomizer<Dimension>::AxisBoundaryCondition
+{
+ public:
+  using Rd = TinyVector<Dimension, double>;
+
+ private:
+  const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
+
+ public:
+  const Rd&
+  direction() const
+  {
+    return m_mesh_line_node_boundary.direction();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_line_node_boundary.nodeList();
+  }
+
+  AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
+    : m_mesh_line_node_boundary(mesh_line_node_boundary)
+  {
+    ;
+  }
+
+  ~AxisBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshRandomizer<Dimension>::FixedBoundaryCondition
+{
+ private:
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
+
+  ~FixedBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshRandomizer<Dimension>::SymmetryBoundaryCondition
+{
+ 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();
+  }
+
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_flat_node_boundary.nodeList();
+  }
+
+  SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+#endif   // MESH_RANDOMIZER_HPP
diff --git a/src/scheme/AcousticSolver.cpp b/src/scheme/AcousticSolver.cpp
index 6973037571b23f861cec5c70ed28fbe67c036fcf..0ad9b7901ceb6dd39cab12d817201e307d035308 100644
--- a/src/scheme/AcousticSolver.cpp
+++ b/src/scheme/AcousticSolver.cpp
@@ -18,15 +18,14 @@ template <size_t Dimension>
 double
 acoustic_dt(const DiscreteFunctionP0<Dimension, double>& c)
 {
-  const std::shared_ptr<const Mesh<Connectivity<Dimension>>> mesh =
-    std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(c.mesh());
+  const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh());
 
-  const auto Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
-  const auto Sj = MeshDataManager::instance().getMeshData(*mesh).sumOverRLjr();
+  const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+  const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
 
-  CellValue<double> local_dt{mesh->connectivity()};
+  CellValue<double> local_dt{mesh.connectivity()};
   parallel_for(
-    mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
+    mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
 
   return min(local_dt);
 }
@@ -89,11 +88,11 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   {
     Assert(rho.mesh() == c.mesh());
 
-    std::shared_ptr mesh = std::dynamic_pointer_cast<const MeshType>(rho.mesh());
+    const MeshType& mesh = dynamic_cast<const MeshType&>(*rho.mesh());
 
-    CellValue<double> rhoc{mesh->connectivity()};
+    CellValue<double> rhoc{mesh.connectivity()};
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId cell_id) { rhoc[cell_id] = rho[cell_id] * c[cell_id]; });
 
     return rhoc;
   }
@@ -247,7 +246,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   }
 
   BoundaryConditionList
-  _getBCList(const std::shared_ptr<const MeshType>& mesh,
+  _getBCList(const MeshType& mesh,
              const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
   {
     BoundaryConditionList bc_list;
@@ -267,9 +266,9 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
       case IBoundaryConditionDescriptor::Type::symmetry: {
         const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor =
           dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
-        for (size_t i_ref_face_list = 0;
-             i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-          const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+        for (size_t i_ref_face_list = 0; i_ref_face_list < mesh.connectivity().template numberOfRefItemList<FaceType>();
+             ++i_ref_face_list) {
+          const auto& ref_face_list = mesh.connectivity().template refItemList<FaceType>(i_ref_face_list);
           const RefId& ref          = ref_face_list.refId();
           if (ref == sym_bc_descriptor.boundaryDescriptor()) {
             SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
@@ -285,8 +284,8 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
         if (dirichlet_bc_descriptor.name() == "velocity") {
           for (size_t i_ref_face_list = 0;
-               i_ref_face_list < mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
-            const auto& ref_face_list = mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+               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};
@@ -296,15 +295,15 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
               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);
+                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") {
           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);
+               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();
@@ -314,9 +313,9 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
               Array<const double> face_values = [&] {
                 if constexpr (Dimension == 1) {
                   return InterpolateItemValue<double(
-                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh->xr(), face_list);
+                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh.xr(), face_list);
                 } else {
-                  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*mesh);
+                  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
 
                   return InterpolateItemValue<double(
                     TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
@@ -394,16 +393,14 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   }
 
   AcousticSolver(SolverType solver_type,
-                 const std::shared_ptr<const MeshType>& p_mesh,
+                 const MeshType& mesh,
                  const DiscreteScalarFunction& rho,
                  const DiscreteScalarFunction& c,
                  const DiscreteVectorFunction& u,
                  const DiscreteScalarFunction& p,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(p_mesh, bc_descriptor_list)}
+    : m_solver_type{solver_type}, m_boundary_condition_list{this->_getBCList(mesh, bc_descriptor_list)}
   {
-    const MeshType& mesh = *p_mesh;
-
     CellValue<const double> rhoc     = this->_getRhoC(rho, c);
     NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(mesh, rhoc);
 
@@ -425,27 +422,27 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
              std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
              std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>>
   apply(const double& dt,
-        const std::shared_ptr<const MeshType>& mesh,
+        const MeshType& mesh,
         const DiscreteFunctionP0<Dimension, double>& rho,
         const DiscreteFunctionP0<Dimension, Rd>& u,
         const DiscreteFunctionP0<Dimension, double>& E) const
   {
-    const auto& cell_to_node_matrix = mesh->connectivity().cellToNodeMatrix();
+    const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
-    NodeValue<Rd> new_xr = copy(mesh->xr());
+    NodeValue<Rd> new_xr = copy(mesh.xr());
     parallel_for(
-      mesh->numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
+      mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * m_ur[r]; });
 
-    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh->shared_connectivity(), new_xr);
+    std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
 
-    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(*mesh).Vj();
+    CellValue<const double> Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
 
     CellValue<double> new_rho = copy(rho.cellValues());
     CellValue<Rd> new_u       = copy(u.cellValues());
     CellValue<double> new_E   = copy(E.cellValues());
 
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) {
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
         const auto& cell_nodes = cell_to_node_matrix[j];
 
         Rd momentum_fluxes   = zero;
@@ -463,7 +460,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
 
     parallel_for(
-      mesh->numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
+      mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
 
     return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
             std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
@@ -488,7 +485,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
       throw NormalError("acoustic solver expects P0 functions");
     }
 
-    return this->apply(dt, std::dynamic_pointer_cast<const MeshType>(i_mesh),
+    return this->apply(dt, dynamic_cast<const MeshType&>(*i_mesh),
                        *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
                        *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
                        *std::dynamic_pointer_cast<const DiscreteScalarFunction>(E));
@@ -502,7 +499,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
                  const std::shared_ptr<const IDiscreteFunction>& p,
                  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
     : AcousticSolver{solver_type,
-                     std::dynamic_pointer_cast<const Mesh<Connectivity<Dimension>>>(mesh),
+                     dynamic_cast<const MeshType&>(*mesh),
                      *std::dynamic_pointer_cast<const DiscreteScalarFunction>(rho),
                      *std::dynamic_pointer_cast<const DiscreteScalarFunction>(c),
                      *std::dynamic_pointer_cast<const DiscreteVectorFunction>(u),
diff --git a/src/scheme/AxisBoundaryConditionDescriptor.hpp b/src/scheme/AxisBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..a03a4bcdd5c2bfbcb23824d3da4e693d4bdbc345
--- /dev/null
+++ b/src/scheme/AxisBoundaryConditionDescriptor.hpp
@@ -0,0 +1,46 @@
+#ifndef AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class AxisBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "axis(" << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::axis;
+  }
+
+  AxisBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+    : m_boundary_descriptor(boundary_descriptor)
+  {
+    ;
+  }
+
+  AxisBoundaryConditionDescriptor(const AxisBoundaryConditionDescriptor&) = delete;
+  AxisBoundaryConditionDescriptor(AxisBoundaryConditionDescriptor&&)      = delete;
+
+  ~AxisBoundaryConditionDescriptor() = default;
+};
+
+#endif   // AXIS_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/FixedBoundaryConditionDescriptor.hpp b/src/scheme/FixedBoundaryConditionDescriptor.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..9167c24beb7a91d1810cff528d23a4f72e39a25a
--- /dev/null
+++ b/src/scheme/FixedBoundaryConditionDescriptor.hpp
@@ -0,0 +1,47 @@
+#ifndef FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class FixedBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "fixed(" << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::string_view m_name;
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+
+ public:
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::fixed;
+  }
+
+  FixedBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor)
+    : m_boundary_descriptor(boundary_descriptor)
+  {}
+
+  FixedBoundaryConditionDescriptor(const FixedBoundaryConditionDescriptor&) = delete;
+  FixedBoundaryConditionDescriptor(FixedBoundaryConditionDescriptor&&)      = delete;
+
+  ~FixedBoundaryConditionDescriptor() = default;
+};
+
+#endif   // FIXED_BOUNDARY_CONDITION_DESCRIPTOR_HPP
diff --git a/src/scheme/IBoundaryConditionDescriptor.hpp b/src/scheme/IBoundaryConditionDescriptor.hpp
index 30a25db66f1278e72eb2bea0ce7a017d8022da65..bf92cd450f000db208dfc87c340f2520b9671757 100644
--- a/src/scheme/IBoundaryConditionDescriptor.hpp
+++ b/src/scheme/IBoundaryConditionDescriptor.hpp
@@ -8,8 +8,10 @@ class IBoundaryConditionDescriptor
  public:
   enum class Type
   {
+    axis,
     dirichlet,
     fourier,
+    fixed,
     free,
     neumann,
     symmetry
diff --git a/src/utils/CMakeLists.txt b/src/utils/CMakeLists.txt
index 7639963d4fc73094f2e79f64e9f9c9383426ac43..19feb273c4a6e2ad09818f467caa9f0f79636b35 100644
--- a/src/utils/CMakeLists.txt
+++ b/src/utils/CMakeLists.txt
@@ -12,6 +12,7 @@ add_library(
   Partitioner.cpp
   PETScWrapper.cpp
   PugsUtils.cpp
+  RandomEngine.cpp
   RevisionInfo.cpp
   SignalManager.cpp
   SLEPcWrapper.cpp)
diff --git a/src/utils/RandomEngine.cpp b/src/utils/RandomEngine.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..508f77a39123492376a5141662f4b3daa41476f9
--- /dev/null
+++ b/src/utils/RandomEngine.cpp
@@ -0,0 +1,66 @@
+
+#include <utils/Messenger.hpp>
+#include <utils/RandomEngine.hpp>
+
+std::unique_ptr<RandomEngine> RandomEngine::m_instance;
+
+bool
+isSynchronized(const RandomEngine& random_engine)
+{
+  const uint64_t current_seed      = random_engine.getCurrentSeed();
+  const uint64_t parallel_max_seed = parallel::allReduceMax(current_seed);
+
+  return parallel::allReduceAnd(current_seed == parallel_max_seed);
+}
+
+uint64_t
+RandomEngine::getCurrentSeed() const
+{
+  std::ostringstream ostr;
+  ostr << m_random_engine;
+
+  std::istringstream istr(ostr.str());
+
+  uint64_t current_seed;
+  istr >> current_seed;
+
+  return current_seed;
+}
+
+void
+RandomEngine::create()
+{
+  m_instance = std::unique_ptr<RandomEngine>(new RandomEngine);
+}
+
+void
+RandomEngine::destroy()
+{
+  m_instance.reset();
+}
+
+RandomEngine::RandomEngine()
+{
+  uint64_t random_seed = std::random_device{}();
+  parallel::broadcast(random_seed, 0);
+
+  m_random_engine = std::default_random_engine(random_seed);
+
+  std::cout << " * setting " << rang::fgB::green << "random seed" << rang::style::reset << " to " << rang::fgB::yellow
+            << random_seed << rang::style::reset << '\n';
+}
+
+void
+RandomEngine::setRandomSeed(const uint64_t random_seed)
+{
+  m_random_engine = std::default_random_engine(random_seed);
+
+  std::cout << " * setting " << rang::fgB::green << "random seed" << rang::style::reset << " to " << rang::fgB::yellow
+            << random_seed << rang::style::reset << '\n';
+}
+
+void
+RandomEngine::resetRandomSeed()
+{
+  m_instance = std::unique_ptr<RandomEngine>(new RandomEngine);
+}
diff --git a/src/utils/RandomEngine.hpp b/src/utils/RandomEngine.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..77af2e4970e7c98f5ceee6ef065cfa284c8d8627
--- /dev/null
+++ b/src/utils/RandomEngine.hpp
@@ -0,0 +1,52 @@
+#ifndef RANDOM_ENGINE_HPP
+#define RANDOM_ENGINE_HPP
+
+#include <utils/PugsAssert.hpp>
+#include <utils/PugsMacros.hpp>
+
+#include <memory>
+#include <random>
+
+class RandomEngine
+{
+ private:
+  static std::unique_ptr<RandomEngine> m_instance;
+
+  std::default_random_engine m_random_engine;
+
+  RandomEngine();
+
+ public:
+  void setRandomSeed(const uint64_t random_seed);
+  void resetRandomSeed();
+
+  friend bool isSynchronized(const RandomEngine& random_engine);
+
+  static void create();
+
+  uint64_t getCurrentSeed() const;
+
+  PUGS_INLINE
+  std::default_random_engine&
+  engine()
+  {
+    return m_random_engine;
+  }
+
+  PUGS_INLINE
+  static RandomEngine&
+  instance()
+  {
+    Assert(m_instance.get() != nullptr, "undefined random engine instance");
+    return *m_instance;
+  }
+
+  static void destroy();
+
+  RandomEngine(const RandomEngine&) = delete;
+  RandomEngine(RandomEngine&&)      = delete;
+
+  ~RandomEngine() = default;
+};
+
+#endif   // RANDOM_ENGINE_HPP