diff --git a/src/mesh/MeshSmoother.cpp b/src/mesh/MeshSmoother.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4e5ed50eb70cf534cc3afbfb1bc8a0d1dea87f63 --- /dev/null +++ b/src/mesh/MeshSmoother.cpp @@ -0,0 +1,423 @@ +#include <mesh/MeshSmoother.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 MeshSmootherHandler::MeshSmoother +{ + 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 smoother"; + 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>) { + if constexpr (Dimension > 1) { + 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 { + throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1"); + } + + } 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 smoothation"); + + 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("smoothed mesh coordinates are not synchronized"); + } +#endif // NDEBUG + + Assert(isSynchronized(random_engine), "seed is not synchronized after mesh smoothation"); + + return shift_r; + } + + public: + std::shared_ptr<const IMesh> + getSmoothedMesh() 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> + getSmoothedMesh(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); + } + + MeshSmoother(const MeshSmoother&) = delete; + MeshSmoother(MeshSmoother&&) = delete; + + MeshSmoother(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)) + {} + + ~MeshSmoother() = default; +}; + +template <size_t Dimension> +class MeshSmootherHandler::MeshSmoother<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 <> +class MeshSmootherHandler::MeshSmoother<1>::AxisBoundaryCondition +{ + public: + AxisBoundaryCondition() = default; + ~AxisBoundaryCondition() = default; +}; + +template <size_t Dimension> +class MeshSmootherHandler::MeshSmoother<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 MeshSmootherHandler::MeshSmoother<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> +MeshSmootherHandler::getSmoothedMesh( + 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>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(); + } + case 2: { + constexpr size_t Dimension = 2; + using MeshType = Mesh<Connectivity<Dimension>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(); + } + case 3: { + constexpr size_t Dimension = 3; + using MeshType = Mesh<Connectivity<Dimension>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} + +std::shared_ptr<const IMesh> +MeshSmootherHandler::getSmoothedMesh( + 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>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(function_symbol_id); + } + case 2: { + constexpr size_t Dimension = 2; + using MeshType = Mesh<Connectivity<Dimension>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(function_symbol_id); + } + case 3: { + constexpr size_t Dimension = 3; + using MeshType = Mesh<Connectivity<Dimension>>; + MeshSmoother smoother(dynamic_cast<const MeshType&>(mesh), bc_descriptor_list); + return smoother.getSmoothedMesh(function_symbol_id); + } + default: { + throw UnexpectedError("invalid mesh dimension"); + } + } +} diff --git a/src/mesh/MeshSmoother.hpp b/src/mesh/MeshSmoother.hpp new file mode 100644 index 0000000000000000000000000000000000000000..32ba5c9e64b85ace4290906edf7da05713bd92b8 --- /dev/null +++ b/src/mesh/MeshSmoother.hpp @@ -0,0 +1,33 @@ +#ifndef MESH_SMOOTHER_HPP +#define MESH_SMOOTHER_HPP + +#include <mesh/IMesh.hpp> +#include <scheme/IBoundaryConditionDescriptor.hpp> + +#include <memory> +#include <vector> + +class FunctionSymbolId; + +class MeshSmootherHandler +{ + private: + template <size_t Dimension> + class MeshSmoother; + + public: + std::shared_ptr<const IMesh> getSmoothedMesh( + const IMesh& mesh, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const; + + std::shared_ptr<const IMesh> getSmoothedMesh( + const IMesh& mesh, + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list, + const FunctionSymbolId& function_symbol_id) const; + + MeshSmootherHandler() = default; + MeshSmootherHandler(MeshSmootherHandler&&) = default; + ~MeshSmootherHandler() = default; +}; + +#endif // MESH_RANDOMIZER_HPP