Add the MeshRandomizer functionality

It provides random perturbations of mesh in 1d, 2d and 3d.
It currently supports fixed and symmetry boundaries.

The obtained mesh is independent of the number of processors

src/language/modules/SchemeModule.cpp

@@ -10,6 +10,7 @@
 #include <mesh/Mesh.hpp>
 #include <mesh/MeshData.hpp>
 #include <mesh/MeshDataManager.hpp>
+#include <mesh/MeshRandomizer.hpp>
 #include <scheme/AcousticSolver.hpp>
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionDescriptorP0.hpp>
@@ -92,6 +93,44 @@ 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("boundaryName",
                             std::make_shared<
                               BuiltinFunctionEmbedder<std::shared_ptr<const IBoundaryDescriptor>(const std::string&)>>(

src/mesh/MeshRandomizer.hpp

+#ifndef MESH_RANDOMIZER_HPP
+#define MESH_RANDOMIZER_HPP
+
+#include <algebra/TinyMatrix.hpp>
+#include <algebra/TinyVector.hpp>
+#include <mesh/Connectivity.hpp>
+#include <mesh/ItemValueUtils.hpp>
+#include <mesh/Mesh.hpp>
+#include <mesh/MeshNodeBoundary.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 SymmetryBoundaryCondition;
+  class FixedBoundaryCondition;
+
+  using BoundaryCondition = std::variant<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 FaceType = [] {
+      if constexpr (Dimension > 1) {
+        return ItemType::face;
+      } else {
+        return ItemType::node;
+      }
+    }();
+
+    for (const auto& bc_descriptor : bc_descriptor_list) {
+      bool is_valid_boundary_condition = true;
+
+      switch (bc_descriptor->type()) {
+      case IBoundaryConditionDescriptor::Type::symmetry: {
+        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()) {
+            SymmetryBoundaryCondition{MeshFlatNodeBoundary<Dimension>(mesh, ref_face_list)};
+            bc_list.push_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()) {
+            MeshNodeBoundary<Dimension> mesh_node_boundary{mesh, ref_face_list};
+
+            const auto& node_list = mesh_node_boundary.nodeList();
+
+            bc_list.push_back(FixedBoundaryCondition{node_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
+  {
+    std::cout << "applying " << m_boundary_condition_list.size() << " boundary conditions\n";
+    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, 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);
+    }
+  }
+
+ public:
+  std::shared_ptr<const MeshType>
+  getRandomizedMesh() const
+  {
+    const ConnectivityType& connectivity = m_given_mesh.connectivity();
+    NodeValue<const Rd> given_xr         = m_given_mesh.xr();
+
+    auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+    EdgeValue<double> edge_length{connectivity};
+    parallel_for(
+      connectivity.numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
+        const auto& edge_nodes = edge_to_node_matrix[edge_id];
+        const NodeId node_0    = edge_nodes[0];
+        const NodeId node_1    = edge_nodes[1];
+        const Rd delta_x       = given_xr[node_0] - given_xr[node_1];
+
+        edge_length[edge_id] = std::sqrt(dot(delta_x, delta_x));
+      });
+
+    auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
+    NodeValue<double> max_delta_xr{connectivity};
+    parallel_for(
+      connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        const auto& node_edges = node_to_edge_matrix[node_id];
+        double max_delta       = edge_length[node_edges[0]];
+
+        for (size_t i_edge = 1; i_edge < node_edges.size(); ++i_edge) {
+          const EdgeId edge_id = node_edges[i_edge];
+          max_delta            = std::min(max_delta, edge_length[edge_id]);
+        }
+        max_delta_xr[node_id] = max_delta;
+      });
+
+    synchronize(max_delta_xr);
+
+    std::uniform_real_distribution<> dis(-0.3, 0.3);
+
+    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();
+
+    NodeValue<Rd> xr{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];
+      for (; i_node_number < node_number; ++i_node_number) {
+        for (size_t j = 0; j < Dimension; ++j) {
+          dis(random_engine.engine());
+        }
+      }
+
+      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());
+      }
+
+      xr[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(xr);
+
+    parallel_for(
+      connectivity.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);
+  }
+
+  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>::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(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
+    : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
+  {
+    ;
+  }
+
+  ~SymmetryBoundaryCondition() = default;
+};
+
+template <size_t Dimension>
+class MeshRandomizer<Dimension>::FixedBoundaryCondition
+{
+ private:
+  const Array<const NodeId> m_node_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_node_list;
+  }
+
+  FixedBoundaryCondition(const Array<const NodeId>& node_list) : m_node_list{node_list} {}
+
+  ~FixedBoundaryCondition() = default;
+};
+
+#endif   // MESH_RANDOMIZER_HPP