Use an handler to encapsulate MeshRandomizer (simplify compilation)

src/language/modules/SchemeModule.cpp

@@ -21,6 +21,7 @@
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
@@ -102,32 +103,8 @@ SchemeModule::SchemeModule()
                               [](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");
-                                }
-                                }
+                                MeshRandomizerHandler handler;
+                                return handler.getRandomizedMesh(*p_mesh, bc_descriptor_list);
                               }
 
                               ));
@@ -142,32 +119,8 @@ SchemeModule::SchemeModule()
                                  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");
-                                }
-                                }
+                                MeshRandomizerHandler handler;
+                                return handler.getRandomizedMesh(*p_mesh, bc_descriptor_list, function_symbol_id);
                               }
 
                               ));

src/mesh/CMakeLists.txt

@@ -21,6 +21,7 @@ add_library(
   MeshFlatNodeBoundary.cpp
   MeshLineNodeBoundary.cpp
   MeshNodeBoundary.cpp
+  MeshRandomizer.cpp
   SynchronizerManager.cpp)
 
 # Additional dependencies

src/mesh/MeshRandomizer.cpp

+#include <mesh/MeshRandomizer.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/MeshFlatNodeBoundary.hpp>
+#include <mesh/MeshLineNodeBoundary.hpp>
+#include <mesh/MeshNodeBoundary.hpp>
+#include <scheme/AxisBoundaryConditionDescriptor.hpp>
+#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/RandomEngine.hpp>
+
+#include <variant>
+
+template <size_t Dimension>
+class MeshRandomizerHandler::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;
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::axis: {
+        if constexpr (Dimension == 1) {
+          bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        } else {
+          bc_list.emplace_back(
+            AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        }
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        bc_list.emplace_back(
+          SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        break;
+      }
+      case IBoundaryConditionDescriptor::Type::fixed: {
+        bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
+        break;
+      }
+      default: {
+        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<> distribution(-0.45, 0.45);
+
+    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) {
+          distribution(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 * distribution(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) {
+        distribution(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 IMesh>
+  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 IMesh>
+  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 MeshRandomizerHandler::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 MeshRandomizerHandler::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 MeshRandomizerHandler::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;
+};
+
+std::shared_ptr<const IMesh>
+MeshRandomizerHandler::getRandomizedMesh(
+  const IMesh& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
+{
+  switch (mesh.dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh();
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh();
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh();
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}
+
+std::shared_ptr<const IMesh>
+MeshRandomizerHandler::getRandomizedMesh(
+  const IMesh& mesh,
+  const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+  const FunctionSymbolId& function_symbol_id) const
+{
+  switch (mesh.dimension()) {
+  case 1: {
+    constexpr size_t Dimension = 1;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh(function_symbol_id);
+  }
+  case 2: {
+    constexpr size_t Dimension = 2;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh(function_symbol_id);
+  }
+  case 3: {
+    constexpr size_t Dimension = 3;
+    using MeshType             = Mesh<Connectivity<Dimension>>;
+    MeshRandomizer randomizer(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list);
+    return randomizer.getRandomizedMesh(function_symbol_id);
+  }
+  default: {
+    throw UnexpectedError("invalid mesh dimension");
+  }
+  }
+}

src/mesh/MeshRandomizer.hpp

 #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/MeshFlatNodeBoundary.hpp>
-#include <mesh/MeshLineNodeBoundary.hpp>
-#include <mesh/MeshNodeBoundary.hpp>
-#include <scheme/AxisBoundaryConditionDescriptor.hpp>
-#include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <mesh/IMesh.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
-#include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
-#include <utils/RandomEngine.hpp>
 
-#include <variant>
+#include <memory>
 #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;
-
-    for (const auto& bc_descriptor : bc_descriptor_list) {
-      switch (bc_descriptor->type()) {
-      case IBoundaryConditionDescriptor::Type::axis: {
-        if constexpr (Dimension == 1) {
-          bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
-        } else {
-          bc_list.emplace_back(
-            AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
-        }
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::symmetry: {
-        bc_list.emplace_back(
-          SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
-        break;
-      }
-      case IBoundaryConditionDescriptor::Type::fixed: {
-        bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
-        break;
-      }
-      default: {
-        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<> distribution(-0.45, 0.45);
-
-    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) {
-          distribution(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 * distribution(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) {
-        distribution(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;
-  }
+class FunctionSymbolId;
 
- public:
-  std::shared_ptr<const MeshType>
-  getRandomizedMesh() const
+class MeshRandomizerHandler
 {
-    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;
+  class MeshRandomizer;
 
  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;
+  std::shared_ptr<const IMesh> getRandomizedMesh(
+    const IMesh& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const;
+
+  std::shared_ptr<const IMesh> getRandomizedMesh(
+    const IMesh& mesh,
+    const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+    const FunctionSymbolId& function_symbol_id) const;
+
+  MeshRandomizerHandler()                        = default;
+  MeshRandomizerHandler(MeshRandomizerHandler&&) = default;
+  ~MeshRandomizerHandler()                       = default;
 };
 
 #endif   // MESH_RANDOMIZER_HPP