From 4c68f366052f7bdc0a61ec90caddd5c935b0f682 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com> Date: Fri, 7 Jul 2023 19:21:20 +0200 Subject: [PATCH] [ci-skip] Begin boundary condition treatment --- src/mesh/MeshSmootherEscobar.cpp | 402 +++++++++++++++---------------- 1 file changed, 189 insertions(+), 213 deletions(-) diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp index f9a18d442..736a0e037 100644 --- a/src/mesh/MeshSmootherEscobar.cpp +++ b/src/mesh/MeshSmootherEscobar.cpp @@ -40,9 +40,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar using BoundaryConditionList = std::vector<BoundaryCondition>; BoundaryConditionList m_boundary_condition_list; + NodeValue<const bool> m_is_fixed_node; + BoundaryConditionList _getBCList(const MeshType& mesh, - const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) + const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const { BoundaryConditionList bc_list; @@ -77,8 +79,83 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar return bc_list; } + NodeValue<const bool> + _getFixedNodes() const + { + NodeValue<bool> is_fixed{m_given_mesh.connectivity()}; + is_fixed.fill(false); + + 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, FixedBoundaryCondition>) or + ((Dimension == 2) and std::is_same_v<BCType, AxisBoundaryCondition>)) { + const Array<const NodeId>& node_list = bc.nodeList(); + parallel_for( + node_list.size(), PUGS_LAMBDA(size_t i_node) { is_fixed[node_list[i_node]] = true; }); + } + }, + boundary_condition); + } + +#warning treat the axis line case in dimension 3 before this + + synchronize(is_fixed); + + NodeValue<int> bc_number{m_given_mesh.connectivity()}; + bc_number.fill(-1); + + { + int bc_id = 0; + + 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 Array<const NodeId>& node_list = bc.nodeList(); + parallel_for( + node_list.size(), PUGS_LAMBDA(size_t i_node) { + const NodeId node_id = node_list[i_node]; + if (not is_fixed[node_id]) { + if (bc_number[node_id] < 0) { + bc_number[node_id] = bc_id; + } else { + bc_number[node_id] = -2; // error tag + } + } + }); + } + }, + boundary_condition); + bc_id++; + } + + if (min(bc_number) < -1) { + throw NormalError("Smoothing boundary conditions overlap. One must precise fixed points and smoothing axis"); + } + } + + bool missing_bc = false; + auto is_boundary = m_given_mesh.connectivity().isBoundaryNode(); + parallel_for(m_given_mesh.numberOfNodes(), [&](NodeId node_id) { + if (is_boundary[node_id]) { + if (not is_fixed[node_id] and (bc_number[node_id] == -1)) { + missing_bc = true; + } + } + }); + + if (parallel::allReduceOr(missing_bc)) { + throw NormalError("Some boundary conditions are missing (boundary nodes without BC)"); + } + + return is_fixed; + } + void - _applyBC(NodeValue<Rd>& shift) const + _applyBC(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const { for (auto&& boundary_condition : m_boundary_condition_list) { std::visit( @@ -96,22 +173,12 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar 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]; + new_xr[node_id] = P * new_xr[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]; - }); + throw NotImplementedError("Escobar: axis boundary conditions"); } else { throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1"); } @@ -121,7 +188,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar parallel_for( node_list.size(), PUGS_LAMBDA(const size_t i_node) { const NodeId node_id = node_list[i_node]; - shift[node_id] = zero; + new_xr[node_id] = old_xr[node_id]; }); } else { @@ -132,94 +199,32 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar } } - 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; - }); - - NodeValue<Rd> shift_r{connectivity}; - - parallel_for( - m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { - const auto& node_cell_list = node_to_cell_matrix[node_id]; - Rd mean_position(zero); - size_t number_of_neighbours = 0; - - 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]; - mean_position += given_xr[cell_node_id]; - number_of_neighbours++; - } - } - } - mean_position = 1. / number_of_neighbours * mean_position; - shift_r[node_id] = mean_position - given_xr[node_id]; - }); - - this->_applyBC(shift_r); - - synchronize(shift_r); - - return shift_r; - } - public: std::shared_ptr<const ItemValueVariant> getQuality() const { - if constexpr (Dimension == 2) { - const ConnectivityType& connectivity = m_given_mesh.connectivity(); - NodeValue<const Rd> xr = m_given_mesh.xr(); + return std::make_shared<ItemValueVariant>(m_given_mesh.xr()); + } + void + _getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const + { + const ConnectivityType& connectivity = m_given_mesh.connectivity(); + + auto is_boundary_node = connectivity.isBoundaryNode(); + auto node_is_owned = connectivity.nodeIsOwned(); + + if constexpr (Dimension == 2) { auto cell_to_node_matrix = connectivity.cellToNodeMatrix(); auto node_to_cell_matrix = connectivity.nodeToCellMatrix(); auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells(); - auto is_boundary_node = connectivity.isBoundaryNode(); NodeValue<double> quality{connectivity}; constexpr double eps = 1E-15; quality.fill(2); - auto f_inner = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double { + auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double { auto cell_list = node_to_cell_matrix[node_id]; auto node_number_in_cell = node_number_in_their_cells[node_id]; @@ -241,8 +246,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]]; const size_t cell_nb_nodes = cell_node_list.size(); - const TinyVector a = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]]; - const TinyVector b = xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]]; + const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]]; + const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]]; const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], // a[1] - x[1], b[1] - x[1]}; @@ -262,9 +267,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); }; - // TinyVector<Dimension> f_gradient = zero; - // TinyMatrix<Dimension> f_hessian = zero; - double final_f = 0; for (size_t i_iter = 0; i_iter < 100; ++i_iter) { @@ -274,8 +276,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]]; const size_t cell_nb_nodes = cell_node_list.size(); - const TinyVector a = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]]; - const TinyVector b = xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]]; + const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]]; + const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]]; const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0], // a[1] - x[1], b[1] - x[1]}; @@ -310,22 +312,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), // trace(Sigma * S_gradient[1])}; - const std::array<TinyMatrix<Dimension>, Dimension> // - Sigma_gradient_old{sigma_gradient[0] * inverse(S) - inverse(S) * S_gradient[0] * Sigma, - sigma_gradient[1] * inverse(S) - inverse(S) * S_gradient[1] * Sigma}; - - const std::array<TinyMatrix<Dimension>, Dimension> // - Sigma_gradient_new{TinyMatrix<Dimension>{0, 1. / std::sin(alpha - 1. / std::tan(alpha)), // - 0, -1}, - TinyMatrix<Dimension>{-1. / std::sin(alpha) + 1. / std::tan(alpha), 0, // - 1, 0}}; - const auto Sigma_gradient = Sigma_gradient_new; - std::cout << "Sigma_gradient_old[0] = " << Sigma_gradient_old[0] << '\n'; - std::cout << "Sigma_gradient_new[0] = " << Sigma_gradient_new[0] << '\n'; - std::cout << "Sigma_gradient_old[1] = " << Sigma_gradient_old[1] << '\n'; - std::cout << "Sigma_gradient_new[1] = " << Sigma_gradient_new[1] << '\n'; - - // TinyVector<Dimension> h_gradient = h / (h - sigma_list[i_cell]) * sigma_gradient; + const std::array<TinyMatrix<Dimension>, Dimension> // + Sigma_gradient{TinyMatrix<Dimension>{0, 1. / std::sin(alpha - 1. / std::tan(alpha)), // + 0, -1}, + TinyMatrix<Dimension>{-1. / std::sin(alpha) + 1. / std::tan(alpha), 0, // + 1, 0}}; TinyVector<Dimension> g{trace(transpose(S) * S_gradient[0]) / S_norm2 // + trace(transpose(Sigma) * Sigma_gradient[0]) / Sigma_norm2 // @@ -360,19 +351,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar f_hessian += f_jr_hessian; } - std::cout << "f = " << f << '\n'; - std::cout << "grad(f) = " << f_gradient << '\n'; - std::cout << "hess(f) = " << f_hessian << " | hess(f)^T -hess(f) = " << transpose(f_hessian) - f_hessian - << '\n'; - - std::cout << "inv(H) = " << inverse(f_hessian) << '\n'; - std::cout << "inv(H)*grad(f) = " << inverse(f_hessian) * f_gradient << '\n'; - - std::cout << rang::fgB::yellow << "x = " << x << " -> " << x - inverse(f_hessian) * f_gradient - << rang::fg::reset << '\n'; - - std::cout << rang::fgB::green << i_iter << ": l2Norm(f_gradient) = " << l2Norm(f_gradient) << rang::fg::reset - << '\n'; if (l2Norm(f_gradient) < 1E-6) { break; } @@ -386,132 +364,128 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar parallel_for( connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) { - // auto cell_list = node_to_cell_matrix[node_id]; - // auto node_number_in_cell = node_number_in_their_cells[node_id]; - - if (is_boundary_node[node_id]) { - quality[node_id] = 1; - } else { - TinyVector x = xr[node_id]; - quality[node_id] = f_inner(node_id, x); - - std::exit(0); - - // TinyMatrix<Dimension> B = identity; + if (node_is_owned[node_id] and not is_boundary_node[node_id]) { + quality[node_id] = smooth(node_id, new_xr[node_id]); } }); - - return std::make_shared<ItemValueVariant>(quality); } else { throw NotImplementedError("Dimension != 2"); } + + this->_applyBC(old_xr, new_xr); } std::shared_ptr<const IMesh> getSmoothedMesh() const { NodeValue<const Rd> given_xr = m_given_mesh.xr(); + NodeValue<Rd> new_xr = copy(given_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]; }); + this->_getNewPositions(given_xr, new_xr); - return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr); + return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), new_xr); } std::shared_ptr<const IMesh> getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const { - NodeValue<const Rd> given_xr = m_given_mesh.xr(); + // 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<const bool> is_displaced = + // InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr); - NodeValue<Rd> xr = this->_getDisplacement(); + // 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]; }); + // 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); + // return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr); + return nullptr; } std::shared_ptr<const IMesh> getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const { - NodeValue<const Rd> given_xr = m_given_mesh.xr(); - - auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); - - NodeValue<bool> is_displaced{m_given_mesh.connectivity()}; - is_displaced.fill(false); - - for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) { - MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]); - const auto cell_list = cell_zone.cellList(); - CellValue<bool> is_zone_cell{m_given_mesh.connectivity()}; - is_zone_cell.fill(false); - parallel_for( - cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; }); - parallel_for( - m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { - auto node_cell_list = node_to_cell_matrix[node_id]; - bool displace = true; - for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) { - const CellId cell_id = node_cell_list[i_node_cell]; - displace &= is_zone_cell[cell_id]; - } - if (displace) { - is_displaced[node_id] = true; - } - }); - } - synchronize(is_displaced); - 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); + // NodeValue<const Rd> given_xr = m_given_mesh.xr(); + + // auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); + + // NodeValue<bool> is_displaced{m_given_mesh.connectivity()}; + // is_displaced.fill(false); + + // for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) { + // MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]); + // const auto cell_list = cell_zone.cellList(); + // CellValue<bool> is_zone_cell{m_given_mesh.connectivity()}; + // is_zone_cell.fill(false); + // parallel_for( + // cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; }); + // parallel_for( + // m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { + // auto node_cell_list = node_to_cell_matrix[node_id]; + // bool displace = true; + // for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) { + // const CellId cell_id = node_cell_list[i_node_cell]; + // displace &= is_zone_cell[cell_id]; + // } + // if (displace) { + // is_displaced[node_id] = true; + // } + // }); + // } + // synchronize(is_displaced); + // 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); + return nullptr; } std::shared_ptr<const IMesh> getSmoothedMesh( const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const { - NodeValue<const Rd> given_xr = m_given_mesh.xr(); + // NodeValue<const Rd> given_xr = m_given_mesh.xr(); - auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); + // auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); - NodeValue<bool> is_displaced{m_given_mesh.connectivity()}; - is_displaced.fill(false); + // NodeValue<bool> is_displaced{m_given_mesh.connectivity()}; + // is_displaced.fill(false); - for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) { - auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>(); + // for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) { + // auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const + // double>>(); - parallel_for( - m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { - auto node_cell_list = node_to_cell_matrix[node_id]; - bool displace = true; - for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) { - const CellId cell_id = node_cell_list[i_node_cell]; - displace &= (is_zone_cell[cell_id] != 0); - } - if (displace) { - is_displaced[node_id] = true; - } - }); - } - synchronize(is_displaced); - NodeValue<Rd> xr = this->_getDisplacement(); + // parallel_for( + // m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { + // auto node_cell_list = node_to_cell_matrix[node_id]; + // bool displace = true; + // for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) { + // const CellId cell_id = node_cell_list[i_node_cell]; + // displace &= (is_zone_cell[cell_id] != 0); + // } + // if (displace) { + // is_displaced[node_id] = true; + // } + // }); + // } + // synchronize(is_displaced); + // 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]; + // }); - 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); - return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr); + return nullptr; } MeshSmootherEscobar(const MeshSmootherEscobar&) = delete; @@ -519,7 +493,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar MeshSmootherEscobar(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)) + : m_given_mesh{given_mesh}, + m_boundary_condition_list{this->_getBCList(given_mesh, bc_descriptor_list)}, + m_is_fixed_node{this->_getFixedNodes()} {} ~MeshSmootherEscobar() = default; -- GitLab