diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp index cc36d64587eb52acaa0e7c25cb11f0c3851c036f..eca67723dd72c733107cdd5b0caeef04ae24e117 100644 --- a/src/mesh/MeshSmootherEscobar.cpp +++ b/src/mesh/MeshSmootherEscobar.cpp @@ -37,6 +37,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar class FixedBoundaryCondition; class SymmetryBoundaryCondition; + static constexpr double m_min_quality = 0.2; + using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>; using BoundaryConditionList = std::vector<BoundaryCondition>; @@ -428,7 +430,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar return non_convex_cell_tag; } - NodeValue<const bool> + [[deprecated]] NodeValue<const bool> _getSmoothedNode(const CellValue<const int>& non_convex_cell_tag, const NodeValue<const bool>& m_is_fixed_node) const { auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); @@ -451,6 +453,20 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar return is_smoothed; } + NodeValue<const bool> + _getSmoothedNode(const NodeValue<const double>& node_quality, const NodeValue<const bool>& m_is_fixed_node) const + { + auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix(); + NodeValue<bool> is_smoothed{m_given_mesh.connectivity()}; + is_smoothed.fill(false); + parallel_for( + m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { + is_smoothed[node_id] = ((not m_is_fixed_node[node_id]) and (node_quality[node_id] < m_min_quality)); + }); + + return is_smoothed; + } + void get_f_grad_hess(const NodeValue<const TinyVector<Dimension>>& old_xr, double& sigma_min, @@ -786,14 +802,14 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) { const CellId cell_id = node_cells[i_cell]; - const size_t i_A = node_number[i_cell]; - const size_t i_B = (i_A + 1) % node_cells.size(); - const size_t i_C = (i_A + node_cells.size() - 1) % node_cells.size(); - const size_t i_D = (i_A + 2) % node_cells.size(); - const size_t i_E = (i_A + node_cells.size() - 2) % node_cells.size(); - const auto cell_nodes = cell_to_node_matrix[cell_id]; + const size_t i_A = node_number[i_cell]; + const size_t i_B = (i_A + 1) % cell_nodes.size(); + const size_t i_C = (i_A + cell_nodes.size() - 1) % cell_nodes.size(); + const size_t i_D = (i_A + 2) % cell_nodes.size(); + const size_t i_E = (i_A + cell_nodes.size() - 2) % cell_nodes.size(); + const Rd A = xr[cell_nodes[i_A]]; const Rd B = xr[cell_nodes[i_B]]; const Rd C = xr[cell_nodes[i_C]]; @@ -802,15 +818,23 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar const Rd AB = B - A; const Rd AC = C - A; - const Rd AD = D - A; - const Rd AE = E - A; + + const Rd CA = A - C; + const Rd CE = E - C; + + const Rd BD = D - B; + const Rd BA = A - B; const double l_AB = l2Norm(AB); const double l_AC = l2Norm(AC); - const double l_AD = l2Norm(AD); - const double l_AE = l2Norm(AE); - if ((l_AB == 0) or (l_AC == 0) or (l_AD == 0) or (l_AE == 0)) { + const double l_CA = l2Norm(CA); + const double l_CE = l2Norm(CE); + + const double l_BD = l2Norm(BD); + const double l_BA = l2Norm(BA); + + if ((l_AB == 0) or (l_AC == 0) or (l_CA == 0) or (l_CE == 0) or (l_BD == 0) or (l_BA == 0)) { min_sin = -1; max_abs_sin = 1; break; @@ -818,39 +842,41 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar const Rd u_AB = 1. / l_AB * AB; const Rd u_AC = 1. / l_AC * AC; - const Rd u_AD = 1. / l_AD * AD; - const Rd u_AE = 1. / l_AE * AE; + const Rd u_CA = 1. / l_CA * CA; + const Rd u_CE = 1. / l_CE * CE; + const Rd u_BD = 1. / l_BD * BD; + const Rd u_BA = 1. / l_BA * BA; const double sin_CAB = det(TinyMatrix{u_AB, u_AC}); - const double sin_DAB = det(TinyMatrix{u_AB, u_AD}); - const double sin_CAE = det(TinyMatrix{u_AE, u_AC}); + const double sin_ECA = det(TinyMatrix{u_CA, u_CE}); + const double sin_ABD = det(TinyMatrix{u_BD, u_BA}); const double abs_sin_CAB = std::abs(sin_CAB); - const double abs_sin_DAB = std::abs(sin_DAB); - const double abs_sin_CAE = std::abs(sin_CAE); + const double abs_sin_ECA = std::abs(sin_ECA); + const double abs_sin_ABD = std::abs(sin_ABD); if (sin_CAB < min_sin) { min_sin = sin_CAB; } - if (sin_DAB < min_sin) { - min_sin = sin_DAB; + if (sin_ECA < min_sin) { + min_sin = sin_ECA; } - if (sin_CAE < min_sin) { - min_sin = sin_CAE; + if (sin_ABD < min_sin) { + min_sin = sin_ABD; } if (abs_sin_CAB > max_abs_sin) { max_abs_sin = abs_sin_CAB; } - if (abs_sin_DAB > max_abs_sin) { - max_abs_sin = abs_sin_DAB; + if (abs_sin_ECA > max_abs_sin) { + max_abs_sin = abs_sin_ECA; } - if (abs_sin_CAE > max_abs_sin) { - max_abs_sin = abs_sin_CAE; + if (abs_sin_ABD > max_abs_sin) { + max_abs_sin = abs_sin_ABD; } } @@ -1046,7 +1072,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar const double dfi_dn = dot(f_gradient_i, n); - // std::cout << "=== i=" << i << " alpha_i=" << alpha_i << " alpha_i_1=" << alpha_i_1 << " f0=" << f0 + // std::cout << "=== i=" << i << " alpha_i=" << alpha_i << " alpha_i_1=" << alpha_i_1 << " f0=" << + // f0 // << " f_i=" << f_i << " f_i_1=" << f_i_1 << '\n'; if ((f_i > f0 + c1 * alpha_i * df0_dn) or ((f_i >= f_i_1) and (i > 1))) { // std::cout << "case 1: zoom[" << alpha_i << "," << alpha_i_1 << "] f_i=" << f_i << " f0=" << f0 @@ -1087,7 +1114,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar // std::cout << "x -> " << x << '\n'; // std::cout << rang::fgB::red << "f(x)=" << f << " f(x+dx)=" << new_f << " delta = " << delta_x // << rang::fg::reset << '\n'; - // std::cout << rang::fgB::magenta << "l2Norm(delta_x) = " << l2Norm(delta_x) << rang::fg::reset << '\n'; + // std::cout << rang::fgB::magenta << "l2Norm(delta_x) = " << l2Norm(delta_x) << rang::fg::reset << + // '\n'; if (l2Norm(delta_x) < 1E-8 * max_edge_lenght) { break; @@ -1165,7 +1193,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar for (CellId cell_id = 0; cell_id < m_given_mesh.numberOfCells(); ++cell_id) { for (size_t i_node = 0; i_node < cell_to_node_matrix[cell_id].size(); ++i_node) { const NodeId node_id = cell_to_node_matrix[cell_id][i_node]; - if (node_quality[node_id] <= 0.15) { + if (node_quality[node_id] <= m_min_quality) { cell_tag[cell_id] = 1; break; } @@ -1180,7 +1208,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar std::cout << " remaining non convex cells: " << nb_non_convex_cells << '\n'; if (nb_non_convex_cells > 0) { - m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node); + m_is_smoothed_node = _getSmoothedNode(node_quality, m_is_fixed_node); this->_getNewPositions(xr, new_xr); xr = new_xr; }