diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp index 81812defdaf2d71101a7c35fb36b466b3a20b33a..d989bff45dffa0c0f361607c1488f75f4dc9f0f4 100644 --- a/src/language/modules/SchemeModule.cpp +++ b/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); } )); diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt index 97b03eda3487e6506c2d18f77a0ac5fd13d9b780..eb8e73a01d9374f07b764d18b0e1dc9a5ece1bf1 100644 --- a/src/mesh/CMakeLists.txt +++ b/src/mesh/CMakeLists.txt @@ -21,6 +21,7 @@ add_library( MeshFlatNodeBoundary.cpp MeshLineNodeBoundary.cpp MeshNodeBoundary.cpp + MeshRandomizer.cpp SynchronizerManager.cpp) # Additional dependencies diff --git a/src/mesh/MeshRandomizer.cpp b/src/mesh/MeshRandomizer.cpp new file mode 100644 index 0000000000000000000000000000000000000000..34087ff69d78f1c2d2bfc83589350c376d5fea70 --- /dev/null +++ b/src/mesh/MeshRandomizer.cpp @@ -0,0 +1,411 @@ +#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"); + } + } +} diff --git a/src/mesh/MeshRandomizer.hpp b/src/mesh/MeshRandomizer.hpp index fb29175a8b6254931e055d773bcf345b940de054..f09aa78634919818fe6fe8f581f2b12976dea3b9 100644 --- a/src/mesh/MeshRandomizer.hpp +++ b/src/mesh/MeshRandomizer.hpp @@ -1,355 +1,33 @@ #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); +class FunctionSymbolId; -#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 +class MeshRandomizerHandler { - public: - using Rd = TinyVector<Dimension, double>; - private: - const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary; + template <size_t Dimension> + class MeshRandomizer; 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; + 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