diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index b4022d660475cf93eb9f9e1cb291622de22dba22..9ba22459f1884535c176f11325ed4e447a8fbf71 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/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&)>>(
diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..7c77b7328c28cbf3d64756ecd59b1d9d28c787ff
--- /dev/null
+++ b/src/mesh/MeshRandomizer.hpp
@@ -0,0 +1,291 @@
+#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