From b68c7986ec34b21a15badc4ef07be0b7305b924c Mon Sep 17 00:00:00 2001 From: Connor Ethan Ouellette Date: Mon, 15 Jun 2026 14:14:16 -0600 Subject: [PATCH 1/3] Prevents circumscribing other nodes in the triangulation --- src/mesh/poly2tri_triangulator.C | 1743 ++++++++++++++---------------- 1 file changed, 824 insertions(+), 919 deletions(-) diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index e1280711f53..6bd39ec61fb 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/src/mesh/poly2tri_triangulator.C @@ -15,7 +15,6 @@ // License along with this library; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA - #include "libmesh/libmesh_config.h" #ifdef LIBMESH_HAVE_POLY2TRI @@ -53,19 +52,18 @@ struct P2TPointCompare } }; -p2t::Point to_p2t(const libMesh::Point & p) +p2t::Point +to_p2t(const libMesh::Point & p) { #if LIBMESH_DIM > 2 - libmesh_error_msg_if - (p(2) != 0, - "Poly2TriTriangulator only supports point sets in the XY plane"); + libmesh_error_msg_if(p(2) != 0, "Poly2TriTriangulator only supports point sets in the XY plane"); #endif return {double(p(0)), double(p(1))}; } -Real distance_from_circumcircle(const Elem & elem, - const Point & p) +Real +distance_from_circumcircle(const Elem & elem, const Point & p) { libmesh_assert_equal_to(elem.n_vertices(), 3); @@ -76,10 +74,8 @@ Real distance_from_circumcircle(const Elem & elem, return p_dist - radius; } - -bool in_circumcircle(const Elem & elem, - const Point & p, - const Real tol = 0) +bool +in_circumcircle(const Elem & elem, const Point & p, const Real tol = 0) { return (distance_from_circumcircle(elem, p) < tol); @@ -96,11 +92,8 @@ bool in_circumcircle(const Elem & elem, */ } - std::pair -can_delaunay_swap(const Elem & elem, - unsigned short side, - Real tol) +can_delaunay_swap(const Elem & elem, unsigned short side, Real tol) { const Elem * neigh = elem.neighbor_ptr(side); if (!neigh) @@ -110,42 +103,37 @@ can_delaunay_swap(const Elem & elem, // What neighbor node does elem not share? for (; nn < 3; ++nn) - { - const Node * neigh_node = neigh->node_ptr(nn); - if (neigh_node == elem.node_ptr(0) || - neigh_node == elem.node_ptr(1) || - neigh_node == elem.node_ptr(2)) - continue; - - // Might we need to do a diagonal swap here? Avoid - // undoing a borderline swap. - if (in_circumcircle(elem, *neigh_node, tol)) - break; - } + { + const Node * neigh_node = neigh->node_ptr(nn); + if (neigh_node == elem.node_ptr(0) || neigh_node == elem.node_ptr(1) || + neigh_node == elem.node_ptr(2)) + continue; + + // Might we need to do a diagonal swap here? Avoid + // undoing a borderline swap. + if (in_circumcircle(elem, *neigh_node, tol)) + break; + } if (nn == 3) return {false, 0}; - const unsigned short n = (side+2)%3; - const RealVectorValue right = - (elem.point((n+1)%3)-elem.point(n)).unit(); - const RealVectorValue mid = - (neigh->point(nn)-elem.point(n)).unit(); - const RealVectorValue left = - (elem.point((n+2)%3)-elem.point(n)).unit(); + const unsigned short n = (side + 2) % 3; + const RealVectorValue right = (elem.point((n + 1) % 3) - elem.point(n)).unit(); + const RealVectorValue mid = (neigh->point(nn) - elem.point(n)).unit(); + const RealVectorValue left = (elem.point((n + 2) % 3) - elem.point(n)).unit(); // If the "middle" vector isn't really in the middle, we can't do a // swap without involving other triangles (or we can't at all if // there's a domain boundary in the way) - if (mid*right < left*right || - left*mid < left*right) + if (mid * right < left * right || left * mid < left * right) return {false, 0}; return {true, nn}; } - -[[maybe_unused]] void libmesh_assert_locally_delaunay(const Elem & elem) +[[maybe_unused]] void +libmesh_assert_locally_delaunay(const Elem & elem) { libmesh_ignore(elem); @@ -158,9 +146,8 @@ can_delaunay_swap(const Elem & elem, } template -inline -void libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), - Container & new_elems) +inline void +libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), Container & new_elems) { libmesh_ignore(new_elems); #ifndef NDEBUG @@ -170,174 +157,168 @@ void libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), libmesh_assert_locally_delaunay(*elem); for (auto & [raw_elem, unique_elem] : new_elems) - { - libmesh_ignore(unique_elem); // avoid warnings on old gcc - libmesh_assert_locally_delaunay(*raw_elem); - } + { + libmesh_ignore(unique_elem); // avoid warnings on old gcc + libmesh_assert_locally_delaunay(*raw_elem); + } #endif } - // Restore a triangulation's Delaunay property, starting with a set of // all triangles that might initially not be locally Delaunay with // their neighbors. template -inline -void restore_delaunay(Container & check_delaunay_on, - BoundaryInfo & boundary_info) +inline void +restore_delaunay(Container & check_delaunay_on, BoundaryInfo & boundary_info) { LOG_SCOPE("restore_delaunay()", "Poly2TriTriangulator"); while (!check_delaunay_on.empty()) + { + Elem & elem = **check_delaunay_on.begin(); + check_delaunay_on.erase(&elem); + for (auto s : make_range(elem.n_sides())) { - Elem & elem = **check_delaunay_on.begin(); - check_delaunay_on.erase(&elem); - for (auto s : make_range(elem.n_sides())) - { - // Can we make a swap here? With what neighbor, with what - // far node? Use a negative tolerance to avoid swapping - // back and forth. - auto [can_swap, nn] = - can_delaunay_swap(elem, s, -TOLERANCE*TOLERANCE); - if (!can_swap) - continue; - - Elem * neigh = elem.neighbor_ptr(s); - - // If we made it here it's time to diagonal swap - const unsigned short n = (s+2)%3; - - const std::array nodes {elem.node_ptr(n), - elem.node_ptr((n+1)%3), neigh->node_ptr(nn), - elem.node_ptr((n+2)%3)}; - - // If we have to swap then either we or any of our neighbors - // might no longer be Delaunay - for (auto ds : make_range(3)) - { - if (elem.neighbor_ptr(ds)) - check_delaunay_on.insert(elem.neighbor_ptr(ds)); - if (neigh->neighbor_ptr(ds)) - check_delaunay_on.insert(neigh->neighbor_ptr(ds)); - } - - // An interior boundary between two newly triangulated - // triangles shouldn't have any bcids - libmesh_assert(!boundary_info.n_boundary_ids(neigh, (nn+1)%3)); - libmesh_assert(!boundary_info.n_boundary_ids(&elem, (n+1)%3)); - - // The two changing boundary sides might have bcids - std::vector bcids; - boundary_info.boundary_ids(&elem, (n+2)%3, bcids); - if (!bcids.empty()) - { - boundary_info.remove_side(&elem, (n+2)%3); - boundary_info.add_side(neigh, (nn+1)%3, bcids); - } + // Can we make a swap here? With what neighbor, with what + // far node? Use a negative tolerance to avoid swapping + // back and forth. + auto [can_swap, nn] = can_delaunay_swap(elem, s, -TOLERANCE * TOLERANCE); + if (!can_swap) + continue; - boundary_info.boundary_ids(neigh, (nn+2)%3, bcids); - if (!bcids.empty()) - { - boundary_info.remove_side(neigh, (nn+2)%3); - boundary_info.add_side(&elem, (n+1)%3, bcids); - } + Elem * neigh = elem.neighbor_ptr(s); - elem.set_node((n+2)%3, nodes[2]); - neigh->set_node((nn+2)%3, nodes[0]); + // If we made it here it's time to diagonal swap + const unsigned short n = (s + 2) % 3; - // No need for a temporary array to swap these around, if we - // do it in the right order. - // - // Watch me neigh->neigh - Elem * neighneigh = neigh->neighbor_ptr((nn+2)%3); - if (neighneigh) - { - unsigned int snn = neighneigh->which_neighbor_am_i(neigh); - neighneigh->set_neighbor(snn, &elem); - } + const std::array nodes{elem.node_ptr(n), + elem.node_ptr((n + 1) % 3), + neigh->node_ptr(nn), + elem.node_ptr((n + 2) % 3)}; - Elem * elemoldneigh = elem.neighbor_ptr((n+2)%3); - if (elemoldneigh) - { - unsigned int seon = elemoldneigh->which_neighbor_am_i(&elem); - elemoldneigh->set_neighbor(seon, neigh); - } + // If we have to swap then either we or any of our neighbors + // might no longer be Delaunay + for (auto ds : make_range(3)) + { + if (elem.neighbor_ptr(ds)) + check_delaunay_on.insert(elem.neighbor_ptr(ds)); + if (neigh->neighbor_ptr(ds)) + check_delaunay_on.insert(neigh->neighbor_ptr(ds)); + } + + // An interior boundary between two newly triangulated + // triangles shouldn't have any bcids + libmesh_assert(!boundary_info.n_boundary_ids(neigh, (nn + 1) % 3)); + libmesh_assert(!boundary_info.n_boundary_ids(&elem, (n + 1) % 3)); + + // The two changing boundary sides might have bcids + std::vector bcids; + boundary_info.boundary_ids(&elem, (n + 2) % 3, bcids); + if (!bcids.empty()) + { + boundary_info.remove_side(&elem, (n + 2) % 3); + boundary_info.add_side(neigh, (nn + 1) % 3, bcids); + } - elem.set_neighbor((n+1)%3, neigh->neighbor_ptr((nn+2)%3)); - neigh->set_neighbor((nn+1)%3, elem.neighbor_ptr((n+2)%3)); - elem.set_neighbor((n+2)%3, neigh); - neigh->set_neighbor((nn+2)%3, &elem); + boundary_info.boundary_ids(neigh, (nn + 2) % 3, bcids); + if (!bcids.empty()) + { + boundary_info.remove_side(neigh, (nn + 2) % 3); + boundary_info.add_side(&elem, (n + 1) % 3, bcids); + } + + elem.set_node((n + 2) % 3, nodes[2]); + neigh->set_node((nn + 2) % 3, nodes[0]); + + // No need for a temporary array to swap these around, if we + // do it in the right order. + // + // Watch me neigh->neigh + Elem * neighneigh = neigh->neighbor_ptr((nn + 2) % 3); + if (neighneigh) + { + unsigned int snn = neighneigh->which_neighbor_am_i(neigh); + neighneigh->set_neighbor(snn, &elem); + } - // Start over after this much change, don't just loop to the - // next neighbor - break; - } + Elem * elemoldneigh = elem.neighbor_ptr((n + 2) % 3); + if (elemoldneigh) + { + unsigned int seon = elemoldneigh->which_neighbor_am_i(&elem); + elemoldneigh->set_neighbor(seon, neigh); + } + + elem.set_neighbor((n + 1) % 3, neigh->neighbor_ptr((nn + 2) % 3)); + neigh->set_neighbor((nn + 1) % 3, elem.neighbor_ptr((n + 2) % 3)); + elem.set_neighbor((n + 2) % 3, neigh); + neigh->set_neighbor((nn + 2) % 3, &elem); + + // Start over after this much change, don't just loop to the + // next neighbor + break; } + } } - -unsigned int segment_intersection(const Elem & elem, - Point & source, - const Point & target, - unsigned int source_side) +unsigned int +segment_intersection(const Elem & elem, + Point & source, + const Point & target, + unsigned int source_side) { libmesh_assert_equal_to(elem.dim(), 2); const auto ns = elem.n_sides(); for (auto s : make_range(ns)) - { - // Don't go backwards just because some FP roundoff said to - if (s == source_side) - continue; + { + // Don't go backwards just because some FP roundoff said to + if (s == source_side) + continue; - const Point v0 = elem.point(s); - const Point v1 = elem.point((s+1)%ns); + const Point v0 = elem.point(s); + const Point v1 = elem.point((s + 1) % ns); - // Calculate intersection parameters (fractions of the distance - // along each segment) - const Real raydx = target(0)-source(0), - raydy = target(1)-source(1), - edgedx = v1(0)-v0(0), - edgedy = v1(1)-v0(1); - const Real denom = edgedx * raydy - edgedy * raydx; + // Calculate intersection parameters (fractions of the distance + // along each segment) + const Real raydx = target(0) - source(0), raydy = target(1) - source(1), edgedx = v1(0) - v0(0), + edgedy = v1(1) - v0(1); + const Real denom = edgedx * raydy - edgedy * raydx; - // divide-by-zero means the segments are parallel - if (denom == 0) - continue; + // divide-by-zero means the segments are parallel + if (denom == 0) + continue; - const Real one_over_denom = 1 / denom; + const Real one_over_denom = 1 / denom; - const Real targetsdx = v1(0)-target(0), - targetsdy = v1(1)-target(1); + const Real targetsdx = v1(0) - target(0), targetsdy = v1(1) - target(1); - const Real t_num = targetsdx * raydy - - targetsdy * raydx; - const Real t = t_num * one_over_denom; + const Real t_num = targetsdx * raydy - targetsdy * raydx; + const Real t = t_num * one_over_denom; - if (t < -TOLERANCE*TOLERANCE || t > 1 + TOLERANCE*TOLERANCE) - continue; + if (t < -TOLERANCE * TOLERANCE || t > 1 + TOLERANCE * TOLERANCE) + continue; - const Real u_num = targetsdx * edgedy - targetsdy * edgedx; - const Real u = u_num * one_over_denom; + const Real u_num = targetsdx * edgedy - targetsdy * edgedx; + const Real u = u_num * one_over_denom; - if (u < -TOLERANCE*TOLERANCE || u > 1 + TOLERANCE*TOLERANCE) - continue; + if (u < -TOLERANCE * TOLERANCE || u > 1 + TOLERANCE * TOLERANCE) + continue; -/* - // Partial workaround for an old poly2tri bug (issue #39): if we - // end up with boundary points that are nearly-collinear but - // infinitesimally concave, p2t::CDT::Triangulate throws a "null - // triangle" exception. So let's try to be infinitesimally - // convex instead. - const Real ray_fraction = (1-u) * (1+TOLERANCE*TOLERANCE); -*/ - const Real ray_fraction = (1-u); - - source(0) += raydx * ray_fraction; - source(1) += raydy * ray_fraction; - return s; - } + /* + // Partial workaround for an old poly2tri bug (issue #39): if we + // end up with boundary points that are nearly-collinear but + // infinitesimally concave, p2t::CDT::Triangulate throws a "null + // triangle" exception. So let's try to be infinitesimally + // convex instead. + const Real ray_fraction = (1-u) * (1+TOLERANCE*TOLERANCE); + */ + const Real ray_fraction = (1 - u); + + source(0) += raydx * ray_fraction; + source(1) += raydy * ray_fraction; + return s; + } return libMesh::invalid_uint; } @@ -351,20 +332,16 @@ namespace libMesh // // Constructor -Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh, - dof_id_type n_boundary_nodes) - : TriangulatorInterface(mesh), - _n_boundary_nodes(n_boundary_nodes), - _refine_bdy_allowed(true) +Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh, dof_id_type n_boundary_nodes) + : TriangulatorInterface(mesh), _n_boundary_nodes(n_boundary_nodes), _refine_bdy_allowed(true) { } - Poly2TriTriangulator::~Poly2TriTriangulator() = default; - // Primary function responsible for performing the triangulation -void Poly2TriTriangulator::triangulate() +void +Poly2TriTriangulator::triangulate() { LOG_SCOPE("triangulate()", "Poly2TriTriangulator"); @@ -388,9 +365,7 @@ void Poly2TriTriangulator::triangulate() // We won't support quads any time soon, or 1D/3D in this interface // ever. - if (_elem_type != TRI3 && - _elem_type != TRI6 && - _elem_type != TRI7) + if (_elem_type != TRI3 && _elem_type != TRI6 && _elem_type != TRI7) libmesh_not_implemented(); // If we have no explicit segments defined, we may get them from @@ -411,11 +386,10 @@ void Poly2TriTriangulator::triangulate() // This is currently done redundantly in parallel; make sure no // processor quits before the others. do - { - libmesh_parallel_only(_mesh.comm()); - this->triangulate_current_points(); - } - while (this->insert_refinement_points()); + { + libmesh_parallel_only(_mesh.comm()); + this->triangulate_current_points(); + } while (this->insert_refinement_points()); libmesh_parallel_only(_mesh.comm()); @@ -442,9 +416,8 @@ void Poly2TriTriangulator::triangulate() _mesh.prepare_for_use(); } - -void Poly2TriTriangulator::set_desired_area_function - (FunctionBase * desired) +void +Poly2TriTriangulator::set_desired_area_function(FunctionBase * desired) { if (desired) _desired_area_func = desired->clone(); @@ -452,17 +425,16 @@ void Poly2TriTriangulator::set_desired_area_function _desired_area_func.reset(); } - -FunctionBase * Poly2TriTriangulator::get_desired_area_function () +FunctionBase * +Poly2TriTriangulator::get_desired_area_function() { return _desired_area_func.get(); } - -bool Poly2TriTriangulator::is_refine_boundary_allowed - (const BoundaryInfo & boundary_info, - const Elem & elem, - unsigned int side) +bool +Poly2TriTriangulator::is_refine_boundary_allowed(const BoundaryInfo & boundary_info, + const Elem & elem, + unsigned int side) { // We should only be calling this on a boundary side libmesh_assert(!elem.neighbor_ptr(side)); @@ -480,14 +452,14 @@ bool Poly2TriTriangulator::is_refine_boundary_allowed // side libmesh_assert(this->_holes); - const boundary_id_type hole_num = bcids[0]-1; + const boundary_id_type hole_num = bcids[0] - 1; libmesh_assert_less(hole_num, this->_holes->size()); const Hole * hole = (*this->_holes)[hole_num]; return hole->refine_boundary_allowed(); } - -void Poly2TriTriangulator::triangulate_current_points() +void +Poly2TriTriangulator::triangulate_current_points() { LOG_SCOPE("triangulate_current_points()", "Poly2TriTriangulator"); @@ -506,15 +478,13 @@ void Poly2TriTriangulator::triangulate_current_points() std::vector> inner_hole_points(n_holes); dof_id_type nn = _mesh.max_node_id(); - libmesh_error_msg_if - (!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!"); + libmesh_error_msg_if(!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!"); // Unless we're using an explicit segments list, we assume node ids // are contiguous here. if (this->segments.empty()) - libmesh_error_msg_if - (_mesh.n_nodes() != nn, - "Poly2TriTriangulator needs contiguous node ids or explicit segments!"); + libmesh_error_msg_if(_mesh.n_nodes() != nn, + "Poly2TriTriangulator needs contiguous node ids or explicit segments!"); // And if we have more nodes than outer boundary points, the rest // may be interior "Steiner points". We use a set here so we can @@ -525,92 +495,85 @@ void Poly2TriTriangulator::triangulate_current_points() // If we were asked to use all mesh nodes as boundary nodes, now's // the time to see how many that is. if (_n_boundary_nodes == DofObject::invalid_id) - { - _n_boundary_nodes = _mesh.n_nodes(); - libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes), - std::distance(_mesh.nodes_begin(), - _mesh.nodes_end())); - - } + { + _n_boundary_nodes = _mesh.n_nodes(); + libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes), + std::distance(_mesh.nodes_begin(), _mesh.nodes_end())); + } else - libmesh_assert_less_equal(_n_boundary_nodes, - _mesh.n_nodes()); + libmesh_assert_less_equal(_n_boundary_nodes, _mesh.n_nodes()); // Prepare poly2tri points for our nodes, sorted into outer boundary // points and interior Steiner points. if (this->segments.empty()) + { + // If we have no segments even after taking elems into account, + // the nodal id ordering defines our outer polyline ordering + for (auto & node : _mesh.node_ptr_range()) { - // If we have no segments even after taking elems into account, - // the nodal id ordering defines our outer polyline ordering - for (auto & node : _mesh.node_ptr_range()) - { - const p2t::Point pt = to_p2t(*node); + const p2t::Point pt = to_p2t(*node); - // If we're out of boundary nodes, the rest are going to be - // Steiner points or hole points - if (node->id() < _n_boundary_nodes) - outer_boundary_points.push_back(pt); - else - steiner_points.insert(pt); + // If we're out of boundary nodes, the rest are going to be + // Steiner points or hole points + if (node->id() < _n_boundary_nodes) + outer_boundary_points.push_back(pt); + else + steiner_points.insert(pt); - // We're not going to support overlapping nodes on the boundary - if (point_node_map.count(pt)) - libmesh_not_implemented(); + // We're not going to support overlapping nodes on the boundary + if (point_node_map.count(pt)) + libmesh_not_implemented(); - point_node_map.emplace(pt, node); - } + point_node_map.emplace(pt, node); } + } // If we have explicit segments defined, that's our outer polyline // ordering: else - { - // Let's make sure our segments are in order - dof_id_type last_id = DofObject::invalid_id; + { + // Let's make sure our segments are in order + dof_id_type last_id = DofObject::invalid_id; - // Add nodes from every segment, in order, to the outer polyline - for (auto [segment_start, segment_end] : this->segments) - { - if (last_id != DofObject::invalid_id) - libmesh_error_msg_if(segment_start != last_id, - "Disconnected triangulator segments"); - last_id = segment_end; + // Add nodes from every segment, in order, to the outer polyline + for (auto [segment_start, segment_end] : this->segments) + { + if (last_id != DofObject::invalid_id) + libmesh_error_msg_if(segment_start != last_id, "Disconnected triangulator segments"); + last_id = segment_end; - Node * node = _mesh.query_node_ptr(segment_start); + Node * node = _mesh.query_node_ptr(segment_start); - libmesh_error_msg_if(!node, - "Triangulator segments reference nonexistent node id " << - segment_start); + libmesh_error_msg_if(!node, + "Triangulator segments reference nonexistent node id " << segment_start); - outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1))); - p2t::Point * pt = &outer_boundary_points.back(); + outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1))); + p2t::Point * pt = &outer_boundary_points.back(); - // We're not going to support overlapping nodes on the boundary - if (point_node_map.count(*pt)) - libmesh_not_implemented_msg - ("Triangulating overlapping boundary nodes is unsupported"); + // We're not going to support overlapping nodes on the boundary + if (point_node_map.count(*pt)) + libmesh_not_implemented_msg("Triangulating overlapping boundary nodes is unsupported"); - point_node_map.emplace(*pt, node); - } + point_node_map.emplace(*pt, node); + } - libmesh_error_msg_if(last_id != this->segments[0].first, - "Non-closed-loop triangulator segments"); + libmesh_error_msg_if(last_id != this->segments[0].first, + "Non-closed-loop triangulator segments"); - // If we have points that aren't in any segments, those will be - // Steiner points - for (auto & node : _mesh.node_ptr_range()) - { - const p2t::Point pt = to_p2t(*node); - if (const auto it = point_node_map.find(pt); - it == point_node_map.end()) - { - steiner_points.insert(pt); - point_node_map.emplace(pt, node); - } - else - libmesh_assert_equal_to(it->second, node); - } + // If we have points that aren't in any segments, those will be + // Steiner points + for (auto & node : _mesh.node_ptr_range()) + { + const p2t::Point pt = to_p2t(*node); + if (const auto it = point_node_map.find(pt); it == point_node_map.end()) + { + steiner_points.insert(pt); + point_node_map.emplace(pt, node); + } + else + libmesh_assert_equal_to(it->second, node); } + } // If we have any elements from a previous triangulation, we're // going to replace them with our own. If we have any elements that @@ -622,23 +585,19 @@ void Poly2TriTriangulator::triangulate_current_points() // triangle. We'll give the outer boundary BC 0, and give holes ids // starting from 1. We've already got the point_node_map to find // nodes, so we can just key on pairs of node ids to identify a side. - std::unordered_map, - boundary_id_type, libMesh::hash> side_boundary_id; + std::unordered_map, boundary_id_type, libMesh::hash> + side_boundary_id; const boundary_id_type outer_bcid = 0; const std::size_t n_outer = outer_boundary_points.size(); for (auto i : make_range(n_outer)) - { - const Node * node1 = - libmesh_map_find(point_node_map, outer_boundary_points[i]), - * node2 = - libmesh_map_find(point_node_map, outer_boundary_points[(i+1)%n_outer]); - - side_boundary_id.emplace(std::make_pair(node1->id(), - node2->id()), - outer_bcid); - } + { + const Node *node1 = libmesh_map_find(point_node_map, outer_boundary_points[i]), + *node2 = libmesh_map_find(point_node_map, outer_boundary_points[(i + 1) % n_outer]); + + side_boundary_id.emplace(std::make_pair(node1->id(), node2->id()), outer_bcid); + } // Create poly2tri triangulator with our mesh points std::vector outer_boundary_pointers(n_outer); @@ -647,7 +606,6 @@ void Poly2TriTriangulator::triangulate_current_points() outer_boundary_pointers.begin(), [](p2t::Point & p) { return &p; }); - // Make sure shims for holes last as long as the CDT does; the // poly2tri headers don't make clear whether or not they're hanging // on to references to these vectors, and it would be reasonable for @@ -658,66 +616,58 @@ void Poly2TriTriangulator::triangulate_current_points() // Add any holes for (auto h : make_range(n_holes)) + { + const Hole * initial_hole = (*_holes)[h]; + auto it = replaced_holes.find(initial_hole); + const Hole & our_hole = (it == replaced_holes.end()) ? *initial_hole : *it->second; + auto & poly2tri_hole = inner_hole_points[h]; + + for (auto i : make_range(our_hole.n_points())) { - const Hole * initial_hole = (*_holes)[h]; - auto it = replaced_holes.find(initial_hole); - const Hole & our_hole = - (it == replaced_holes.end()) ? - *initial_hole : *it->second; - auto & poly2tri_hole = inner_hole_points[h]; - - for (auto i : make_range(our_hole.n_points())) - { - Point p = our_hole.point(i); - poly2tri_hole.emplace_back(to_p2t(p)); + Point p = our_hole.point(i); + poly2tri_hole.emplace_back(to_p2t(p)); - const auto & pt = poly2tri_hole.back(); + const auto & pt = poly2tri_hole.back(); - // This won't be a steiner point. - steiner_points.erase(pt); + // This won't be a steiner point. + steiner_points.erase(pt); - // If we see a hole point already in the mesh, we'll share - // that node. This might be a problem if it's a boundary - // node, but it might just be the same hole point already - // added during a previous triangulation refinement step. - if (point_node_map.count(pt)) - { - libmesh_assert_equal_to - (point_node_map[pt], - _mesh.query_node_ptr(point_node_map[pt]->id())); - } - else - { - Node * node = _mesh.add_point(p, nn++); - point_node_map[pt] = node; - } - } + // If we see a hole point already in the mesh, we'll share + // that node. This might be a problem if it's a boundary + // node, but it might just be the same hole point already + // added during a previous triangulation refinement step. + if (point_node_map.count(pt)) + { + libmesh_assert_equal_to(point_node_map[pt], _mesh.query_node_ptr(point_node_map[pt]->id())); + } + else + { + Node * node = _mesh.add_point(p, nn++); + point_node_map[pt] = node; + } + } - const boundary_id_type inner_bcid = h+1; - const std::size_t n_inner = poly2tri_hole.size(); + const boundary_id_type inner_bcid = h + 1; + const std::size_t n_inner = poly2tri_hole.size(); - for (auto i : make_range(n_inner)) - { - const Node * node1 = - libmesh_map_find(point_node_map, poly2tri_hole[i]), - * node2 = - libmesh_map_find(point_node_map, poly2tri_hole[(i+1)%n_inner]); - - side_boundary_id.emplace(std::make_pair(node1->id(), - node2->id()), - inner_bcid); - } + for (auto i : make_range(n_inner)) + { + const Node *node1 = libmesh_map_find(point_node_map, poly2tri_hole[i]), + *node2 = libmesh_map_find(point_node_map, poly2tri_hole[(i + 1) % n_inner]); + + side_boundary_id.emplace(std::make_pair(node1->id(), node2->id()), inner_bcid); + } - auto & poly2tri_ptrs = inner_hole_pointers[h]; - poly2tri_ptrs.resize(n_inner); + auto & poly2tri_ptrs = inner_hole_pointers[h]; + poly2tri_ptrs.resize(n_inner); - std::transform(poly2tri_hole.begin(), - poly2tri_hole.end(), - poly2tri_ptrs.begin(), - [](p2t::Point & p) { return &p; }); + std::transform(poly2tri_hole.begin(), + poly2tri_hole.end(), + poly2tri_ptrs.begin(), + [](p2t::Point & p) { return &p; }); - cdt.AddHole(poly2tri_ptrs); - } + cdt.AddHole(poly2tri_ptrs); + } // Add any steiner points. We had them in a set, but post-C++11 // that won't give us non-const element access (even if we @@ -743,60 +693,52 @@ void Poly2TriTriangulator::triangulate_current_points() // Add the triangles to our Mesh data structure. for (auto ptri_ptr : triangles) + { + p2t::Triangle & ptri = *ptri_ptr; + + // We always use TRI3 here, since that's what we have nodes for; + // if we need a higher order we can convert at the end. + auto elem = Elem::build_with_id(TRI3, next_id++); + for (auto v : make_range(3)) { - p2t::Triangle & ptri = *ptri_ptr; + const p2t::Point & vertex = *ptri.GetPoint(v); - // We always use TRI3 here, since that's what we have nodes for; - // if we need a higher order we can convert at the end. - auto elem = Elem::build_with_id(TRI3, next_id++); - for (auto v : make_range(3)) - { - const p2t::Point & vertex = *ptri.GetPoint(v); + Node * node = libmesh_map_find(point_node_map, vertex); + libmesh_assert(node); + elem->set_node(v, node); + } - Node * node = libmesh_map_find(point_node_map, vertex); - libmesh_assert(node); - elem->set_node(v, node); - } + // We expect a consistent triangle orientation + libmesh_assert(!elem->is_flipped()); - // We expect a consistent triangle orientation - libmesh_assert(!elem->is_flipped()); + Elem * added_elem = _mesh.add_elem(std::move(elem)); - Elem * added_elem = _mesh.add_elem(std::move(elem)); + for (auto v : make_range(3)) + { + const Node &node1 = added_elem->node_ref(v), &node2 = added_elem->node_ref((v + 1) % 3); - for (auto v : make_range(3)) - { - const Node & node1 = added_elem->node_ref(v), - & node2 = added_elem->node_ref((v+1)%3); - - auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id())); - if (it == side_boundary_id.end()) - it = side_boundary_id.find(std::make_pair(node2.id(), node1.id())); - if (it != side_boundary_id.end()) - boundary_info.add_side(added_elem, v, it->second); - } + auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id())); + if (it == side_boundary_id.end()) + it = side_boundary_id.find(std::make_pair(node2.id(), node1.id())); + if (it != side_boundary_id.end()) + boundary_info.add_side(added_elem, v, it->second); } + } } - - -bool Poly2TriTriangulator::insert_refinement_points() +bool +Poly2TriTriangulator::insert_refinement_points() { LOG_SCOPE("insert_refinement_points()", "Poly2TriTriangulator"); if (this->minimum_angle() != 0) libmesh_not_implemented(); + BoundaryInfo & boundary_info = _mesh.get_boundary_info(); // We need neighbor pointers for ray casting and cavity finding UnstructuredMesh & mesh = dynamic_cast(this->_mesh); mesh.find_neighbors(); - if (this->desired_area() == 0 && - this->get_desired_area_function() == nullptr && - !this->has_auto_area_function()) - return false; - - BoundaryInfo & boundary_info = _mesh.get_boundary_info(); - // We won't immediately add these, lest we invalidate iterators on a // ReplicatedMesh. They'll still be in the mesh neighbor topology // for the purpose of doing Delaunay cavity stuff, so we need to @@ -814,8 +756,10 @@ bool Poly2TriTriangulator::insert_refinement_points() // not-yet-added elements so we can't iterate based on proper // element ids ... but we can set a temporary element id to use for // the purpose. - struct cmp { - bool operator()(Elem * a, Elem * b) const { + struct cmp + { + bool operator()(Elem * a, Elem * b) const + { libmesh_assert(a == b || a->id() != b->id()); return (a->id() < b->id()); } @@ -833,14 +777,17 @@ bool Poly2TriTriangulator::insert_refinement_points() // about the isomorphisms! If a triangle's nodes are permuted on // one processor vs another that's an issue. So sort our input // carefully. - std::set all_elems - { mesh.elements_begin(), mesh.elements_end(), comp }; + std::set all_elems{mesh.elements_begin(), mesh.elements_end(), comp}; restore_delaunay(all_elems, boundary_info); libmesh_assert_delaunay(mesh, new_elems); } + if (this->desired_area() == 0 && this->get_desired_area_function() == nullptr && + !this->has_auto_area_function()) + return false; + // Map of which points follow which in the boundary polylines. If // we have to add new boundary points, we'll use this to construct // an updated this->segments to retriangulate with. If we have to @@ -858,10 +805,10 @@ bool Poly2TriTriangulator::insert_refinement_points() #ifdef DEBUG std::unordered_set mesh_points; for (const Node * node : mesh.node_ptr_range()) - { - libmesh_assert(!mesh_points.count(*node)); - mesh_points.insert(*node); - } + { + libmesh_assert(!mesh_points.count(*node)); + mesh_points.insert(*node); + } #endif auto add_point = [&mesh, @@ -869,621 +816,586 @@ bool Poly2TriTriangulator::insert_refinement_points() &mesh_points, #endif &nn](const Point & p) - { + { #ifdef DEBUG - libmesh_assert(!mesh_points.count(p)); - mesh_points.insert(p); + libmesh_assert(!mesh_points.count(p)); + mesh_points.insert(p); #endif - return mesh.add_point(p, nn++); - }; + return mesh.add_point(p, nn++); + }; for (auto & elem : mesh.element_ptr_range()) - { - // element_ptr_range skips deleted elements ... right? - libmesh_assert(elem); - libmesh_assert(elem->valid_id()); + { + // element_ptr_range skips deleted elements ... right? + libmesh_assert(elem); + libmesh_assert(elem->valid_id()); - // We only handle triangles in our triangulation - libmesh_assert_equal_to(elem->level(), 0u); - libmesh_assert_equal_to(elem->type(), TRI3); + // We only handle triangles in our triangulation + libmesh_assert_equal_to(elem->level(), 0u); + libmesh_assert_equal_to(elem->type(), TRI3); - // If this triangle is as small as we desire, move along - if (!should_refine_elem(*elem)) - continue; + // If this triangle is as small as we desire, move along + if (!should_refine_elem(*elem)) + continue; - // Otherwise add a Steiner point. We'd like to add the - // circumcenter ... - Point new_pt = elem->quasicircumcenter(); + // Otherwise add a Steiner point. We'd like to add the + // circumcenter ... + Point new_pt = elem->quasicircumcenter(); - // And to give it a node; - Node * new_node = nullptr; + // And to give it a node; + Node * new_node = nullptr; - // But that might be outside our triangle, or even outside the - // boundary. We'll find a triangle that should contain our new - // point - Elem * cavity_elem = elem; // Start looking at elem anyway + // But that might be outside our triangle, or even outside the + // boundary. We'll find a triangle that should contain our new + // point + Elem * cavity_elem = elem; // Start looking at elem anyway - // We'll refine a boundary later if necessary. - auto boundary_refine = [this, &next_boundary_node, - &cavity_elem, &new_node] - (unsigned int side) + // We'll refine a boundary later if necessary. + auto boundary_refine = [this, &next_boundary_node, &cavity_elem, &new_node](unsigned int side) + { + libmesh_ignore(this); // Only used in dbg/devel + libmesh_assert(new_node); + libmesh_assert(new_node->valid_id()); + + Node *old_segment_start = cavity_elem->node_ptr(side), + *old_segment_end = cavity_elem->node_ptr((side + 1) % 3); + libmesh_assert(old_segment_start); + libmesh_assert(old_segment_start->valid_id()); + libmesh_assert(old_segment_end); + libmesh_assert(old_segment_end->valid_id()); + + if (auto it = next_boundary_node.find(*old_segment_start); it != next_boundary_node.end()) { - libmesh_ignore(this); // Only used in dbg/devel - libmesh_assert(new_node); - libmesh_assert(new_node->valid_id()); - - Node * old_segment_start = cavity_elem->node_ptr(side), - * old_segment_end = cavity_elem->node_ptr((side+1)%3); - libmesh_assert(old_segment_start); - libmesh_assert(old_segment_start->valid_id()); - libmesh_assert(old_segment_end); - libmesh_assert(old_segment_end->valid_id()); - - if (auto it = next_boundary_node.find(*old_segment_start); - it != next_boundary_node.end()) - { - libmesh_assert(it->second == old_segment_end); - it->second = new_node; - } - else - { - // This would be an O(N) sanity check if we already - // have a segments vector or any holes. :-P - libmesh_assert(!this->segments.empty() || - (_holes && !_holes->empty()) || - (old_segment_end->id() == - old_segment_start->id() + 1)); - next_boundary_node[*old_segment_start] = new_node; - } + libmesh_assert(it->second == old_segment_end); + it->second = new_node; + } + else + { + // This would be an O(N) sanity check if we already + // have a segments vector or any holes. :-P + libmesh_assert(!this->segments.empty() || (_holes && !_holes->empty()) || + (old_segment_end->id() == old_segment_start->id() + 1)); + next_boundary_node[*old_segment_start] = new_node; + } + + next_boundary_node[*new_node] = old_segment_end; + }; + + // Let's find a triangle containing our new point, or at least + // containing the end of a ray leading from our current triangle + // to the new point. + Point ray_start = elem->vertex_average(); - next_boundary_node[*new_node] = old_segment_end; - }; + // What side are we coming from, and what side are we going to? + unsigned int source_side = invalid_uint; + unsigned int side = invalid_uint; - // Let's find a triangle containing our new point, or at least - // containing the end of a ray leading from our current triangle - // to the new point. - Point ray_start = elem->vertex_average(); + while (!cavity_elem->contains_point(new_pt)) + { + side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side); - // What side are we coming from, and what side are we going to? - unsigned int source_side = invalid_uint; - unsigned int side = invalid_uint; + libmesh_assert_not_equal_to(side, invalid_uint); - while (!cavity_elem->contains_point(new_pt)) + Elem * neigh = cavity_elem->neighbor_ptr(side); + // If we're on a boundary, stop there. Refine the boundary + // if we're allowed, the boundary element otherwise. + if (!neigh) + { + if (this->is_refine_boundary_allowed(boundary_info, *cavity_elem, side)) { - side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side); + new_pt = ray_start; + new_node = add_point(new_pt); + boundary_refine(side); + } + else + { + // Should we just add the vertex average of the + // boundary element, to minimize the number of + // slivers created? + // + // new_pt = cavity_elem->vertex_average(); + // + // That works for a while, but it + // seems to be able to "run away" and leave us with + // crazy slivers on boundaries if we push interior + // refinement too far while disabling boundary + // refinement. + // + // Let's go back to refining our original problem + // element. + cavity_elem = elem; + new_pt = cavity_elem->vertex_average(); + new_node = add_point(new_pt); + + // This was going to be a side refinement but it's + // now an internal refinement + side = invalid_uint; + } - libmesh_assert_not_equal_to (side, invalid_uint); + break; + } - Elem * neigh = cavity_elem->neighbor_ptr(side); - // If we're on a boundary, stop there. Refine the boundary - // if we're allowed, the boundary element otherwise. - if (!neigh) - { - if (this->is_refine_boundary_allowed(boundary_info, - *cavity_elem, - side)) - { - new_pt = ray_start; - new_node = add_point(new_pt); - boundary_refine(side); - } - else - { - // Should we just add the vertex average of the - // boundary element, to minimize the number of - // slivers created? - // - // new_pt = cavity_elem->vertex_average(); - // - // That works for a while, but it - // seems to be able to "run away" and leave us with - // crazy slivers on boundaries if we push interior - // refinement too far while disabling boundary - // refinement. - // - // Let's go back to refining our original problem - // element. - cavity_elem = elem; - new_pt = cavity_elem->vertex_average(); - new_node = add_point(new_pt); - - // This was going to be a side refinement but it's - // now an internal refinement - side = invalid_uint; - } + source_side = neigh->which_neighbor_am_i(cavity_elem); + cavity_elem = neigh; + side = invalid_uint; + } - break; - } + // If we're ready to create a new node and we're not on a + // boundary ... should we be? We don't want to create any + // sliver elements or confuse poly2tri or anything. + if (side == invalid_uint && !new_node) + { + unsigned int worst_side = libMesh::invalid_uint; + Real worst_cos = 1; + for (auto s : make_range(3u)) + { + // We never snap to a non-domain-boundary + if (cavity_elem->neighbor_ptr(s)) + continue; - source_side = neigh->which_neighbor_am_i(cavity_elem); - cavity_elem = neigh; - side = invalid_uint; + Real ax = cavity_elem->point(s)(0) - new_pt(0), ay = cavity_elem->point(s)(1) - new_pt(1), + bx = cavity_elem->point((s + 1) % 3)(0) - new_pt(0), + by = cavity_elem->point((s + 1) % 3)(1) - new_pt(1); + const Real my_cos = + (ax * bx + ay * by) / std::sqrt(ax * ax + ay * ay) / std::sqrt(bx * bx + by * by); + + if (my_cos < worst_cos) + { + worst_side = s; + worst_cos = my_cos; } + } + + // If we'd create a sliver element on the side, let's just + // refine the side instead, if we're allowed. + if (worst_cos < -0.6) // -0.5 is the best we could enforce? + { + side = worst_side; - // If we're ready to create a new node and we're not on a - // boundary ... should we be? We don't want to create any - // sliver elements or confuse poly2tri or anything. - if (side == invalid_uint && !new_node) + if (this->is_refine_boundary_allowed(boundary_info, *cavity_elem, side)) { - unsigned int worst_side = libMesh::invalid_uint; - Real worst_cos = 1; - for (auto s : make_range(3u)) - { - // We never snap to a non-domain-boundary - if (cavity_elem->neighbor_ptr(s)) - continue; - - Real ax = cavity_elem->point(s)(0) - new_pt(0), - ay = cavity_elem->point(s)(1) - new_pt(1), - bx = cavity_elem->point((s+1)%3)(0) - new_pt(0), - by = cavity_elem->point((s+1)%3)(1) - new_pt(1); - const Real my_cos = (ax*bx+ay*by) / - std::sqrt(ax*ax+ay*ay) / - std::sqrt(bx*bx+by*by); - - if (my_cos < worst_cos) - { - worst_side = s; - worst_cos = my_cos; - } - } + // Let's just try bisecting for now + new_pt = (cavity_elem->point(side) + cavity_elem->point((side + 1) % 3)) / 2; + new_node = add_point(new_pt); + boundary_refine(side); + } + else // Do the best we can under these restrictions + { + new_pt = cavity_elem->vertex_average(); + new_node = add_point(new_pt); - // If we'd create a sliver element on the side, let's just - // refine the side instead, if we're allowed. - if (worst_cos < -0.6) // -0.5 is the best we could enforce? - { - side = worst_side; - - if (this->is_refine_boundary_allowed(boundary_info, - *cavity_elem, - side)) - { - // Let's just try bisecting for now - new_pt = (cavity_elem->point(side) + - cavity_elem->point((side+1)%3)) / 2; - new_node = add_point(new_pt); - boundary_refine(side); - } - else // Do the best we can under these restrictions - { - new_pt = cavity_elem->vertex_average(); - new_node = add_point(new_pt); - - // This was going to be a side refinement but it's - // now an internal refinement - side = invalid_uint; - } - } - else - new_node = add_point(new_pt); + // This was going to be a side refinement but it's + // now an internal refinement + side = invalid_uint; } + } else - libmesh_assert(new_node); + new_node = add_point(new_pt); + } + else + libmesh_assert(new_node); - // Find the Delaunay cavity around the new point. - std::set cavity(comp); + // Find the Delaunay cavity around the new point. + std::set cavity(comp); - std::set unchecked_cavity ({cavity_elem}, comp); - while (!unchecked_cavity.empty()) + std::set unchecked_cavity({cavity_elem}, comp); + while (!unchecked_cavity.empty()) + { + std::set checking_cavity(comp); + checking_cavity.swap(unchecked_cavity); + for (Elem * checking_elem : checking_cavity) + { + for (auto s : make_range(3u)) { - std::set checking_cavity(comp); - checking_cavity.swap(unchecked_cavity); - for (Elem * checking_elem : checking_cavity) - { - for (auto s : make_range(3u)) - { - Elem * neigh = checking_elem->neighbor_ptr(s); - if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh)) - continue; - - if (in_circumcircle(*neigh, new_pt, TOLERANCE*TOLERANCE)) - unchecked_cavity.insert(neigh); - } - } + Elem * neigh = checking_elem->neighbor_ptr(s); + if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh)) + continue; - libmesh_merge_move(cavity, checking_cavity); + if (in_circumcircle(*neigh, new_pt, TOLERANCE * TOLERANCE)) + unchecked_cavity.insert(neigh); } + } - // Retriangulate the Delaunay cavity. - // Each of our cavity triangle edges that are exterior to the - // cavity will be a source of one new triangle. + libmesh_merge_move(cavity, checking_cavity); + } + + // Retriangulate the Delaunay cavity. + // Each of our cavity triangle edges that are exterior to the + // cavity will be a source of one new triangle. - // Set of elements that might need Delaunay swaps - std::set check_delaunay_on(comp); + // Set of elements that might need Delaunay swaps + std::set check_delaunay_on(comp); - // Keep maps for doing neighbor pointer assignment. Not going - // to iterate through these so hashing pointers is fine. - std::unordered_map> - neighbors_CCW, neighbors_CW; + // Keep maps for doing neighbor pointer assignment. Not going + // to iterate through these so hashing pointers is fine. + std::unordered_map> neighbors_CCW, neighbors_CW; - for (Elem * old_elem : cavity) + for (Elem * old_elem : cavity) + { + for (auto s : make_range(3u)) + { + Elem * neigh = old_elem->neighbor_ptr(s); + if (!neigh || !cavity.count(neigh)) { - for (auto s : make_range(3u)) + Node *node_CW = old_elem->node_ptr(s), *node_CCW = old_elem->node_ptr((s + 1) % 3); + + auto set_neighbors = [&neighbors_CW, &neighbors_CCW, &node_CW, &node_CCW, &boundary_info]( + Elem * new_neigh, boundary_id_type bcid) + { + // Set clockwise neighbor and vice-versa if possible + if (const auto CW_it = neighbors_CW.find(node_CW); CW_it == neighbors_CW.end()) + { + libmesh_assert(!neighbors_CCW.count(node_CW)); + neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid); + } + else { - Elem * neigh = old_elem->neighbor_ptr(s); - if (!neigh || !cavity.count(neigh)) - { - Node * node_CW = old_elem->node_ptr(s), - * node_CCW = old_elem->node_ptr((s+1)%3); - - auto set_neighbors = - [&neighbors_CW, &neighbors_CCW, &node_CW, - &node_CCW, &boundary_info] - (Elem * new_neigh, boundary_id_type bcid) - { - // Set clockwise neighbor and vice-versa if possible - if (const auto CW_it = neighbors_CW.find(node_CW); - CW_it == neighbors_CW.end()) - { - libmesh_assert(!neighbors_CCW.count(node_CW)); - neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid); - } - else - { - Elem * neigh_CW = CW_it->second.first; - if (new_neigh) - { - new_neigh->set_neighbor(0, neigh_CW); - boundary_id_type bcid_CW = CW_it->second.second; - if (bcid_CW != BoundaryInfo::invalid_id) - boundary_info.add_side(new_neigh, 0, bcid_CW); - - } - if (neigh_CW) - { - neigh_CW->set_neighbor(2, new_neigh); - if (bcid != BoundaryInfo::invalid_id) - boundary_info.add_side(neigh_CW, 2, bcid); - } - neighbors_CW.erase(CW_it); - } - - // Set counter-CW neighbor and vice-versa if possible - if (const auto CCW_it = neighbors_CCW.find(node_CCW); - CCW_it == neighbors_CCW.end()) - { - libmesh_assert(!neighbors_CW.count(node_CCW)); - neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid); - } - else - { - Elem * neigh_CCW = CCW_it->second.first; - if (new_neigh) - { - boundary_id_type bcid_CCW = CCW_it->second.second; - new_neigh->set_neighbor(2, neigh_CCW); - if (bcid_CCW != BoundaryInfo::invalid_id) - boundary_info.add_side(new_neigh, 2, bcid_CCW); - } - if (neigh_CCW) - { - neigh_CCW->set_neighbor(0, new_neigh); - if (bcid != BoundaryInfo::invalid_id) - boundary_info.add_side(neigh_CCW, 0, bcid); - } - neighbors_CCW.erase(CCW_it); - } - }; - - // We aren't going to try to add a sliver element if we - // have a new boundary node here. We do need to - // keep track of other elements' neighbors, though. - if (old_elem == cavity_elem && - s == side) - { - std::vector bcids; - boundary_info.boundary_ids(old_elem, s, bcids); - libmesh_assert_equal_to(bcids.size(), 1); - set_neighbors(nullptr, bcids[0]); - continue; - } - - auto new_elem = Elem::build_with_id(TRI3, ne++); - new_elem->set_node(0, new_node); - new_elem->set_node(1, node_CW); - new_elem->set_node(2, node_CCW); - libmesh_assert(!new_elem->is_flipped()); - - // Set in-and-out-of-cavity neighbor pointers - new_elem->set_neighbor(1, neigh); - if (neigh) - { - const unsigned int neigh_s = - neigh->which_neighbor_am_i(old_elem); - neigh->set_neighbor(neigh_s, new_elem.get()); - } - else - { - std::vector bcids; - boundary_info.boundary_ids(old_elem, s, bcids); - boundary_info.add_side(new_elem.get(), 1, bcids); - } - - // Set in-cavity neighbors' neighbor pointers - set_neighbors(new_elem.get(), BoundaryInfo::invalid_id); - - // C++ allows function argument evaluation in any - // order, but we need get() to precede move - Elem * new_elem_ptr = new_elem.get(); - new_elems.emplace(new_elem_ptr, std::move(new_elem)); - - check_delaunay_on.insert(new_elem_ptr); - } + Elem * neigh_CW = CW_it->second.first; + if (new_neigh) + { + new_neigh->set_neighbor(0, neigh_CW); + boundary_id_type bcid_CW = CW_it->second.second; + if (bcid_CW != BoundaryInfo::invalid_id) + boundary_info.add_side(new_neigh, 0, bcid_CW); + } + if (neigh_CW) + { + neigh_CW->set_neighbor(2, new_neigh); + if (bcid != BoundaryInfo::invalid_id) + boundary_info.add_side(neigh_CW, 2, bcid); + } + neighbors_CW.erase(CW_it); } - boundary_info.remove(old_elem); - } + // Set counter-CW neighbor and vice-versa if possible + if (const auto CCW_it = neighbors_CCW.find(node_CCW); CCW_it == neighbors_CCW.end()) + { + libmesh_assert(!neighbors_CW.count(node_CCW)); + neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid); + } + else + { + Elem * neigh_CCW = CCW_it->second.first; + if (new_neigh) + { + boundary_id_type bcid_CCW = CCW_it->second.second; + new_neigh->set_neighbor(2, neigh_CCW); + if (bcid_CCW != BoundaryInfo::invalid_id) + boundary_info.add_side(new_neigh, 2, bcid_CCW); + } + if (neigh_CCW) + { + neigh_CCW->set_neighbor(0, new_neigh); + if (bcid != BoundaryInfo::invalid_id) + boundary_info.add_side(neigh_CCW, 0, bcid); + } + neighbors_CCW.erase(CCW_it); + } + }; - // Now that we're done using our cavity elems (including with a - // cavity.find() that used a comparator that dereferences the - // elements!) it's safe to delete them. - for (Elem * old_elem : cavity) - { - if (const auto it = new_elems.find(old_elem); - it == new_elems.end()) - mesh.delete_elem(old_elem); + // We aren't going to try to add a sliver element if we + // have a new boundary node here. We do need to + // keep track of other elements' neighbors, though. + if (old_elem == cavity_elem && s == side) + { + std::vector bcids; + boundary_info.boundary_ids(old_elem, s, bcids); + libmesh_assert_equal_to(bcids.size(), 1); + set_neighbors(nullptr, bcids[0]); + continue; + } + + auto new_elem = Elem::build_with_id(TRI3, ne++); + new_elem->set_node(0, new_node); + new_elem->set_node(1, node_CW); + new_elem->set_node(2, node_CCW); + libmesh_assert(!new_elem->is_flipped()); + + // Set in-and-out-of-cavity neighbor pointers + new_elem->set_neighbor(1, neigh); + if (neigh) + { + const unsigned int neigh_s = neigh->which_neighbor_am_i(old_elem); + neigh->set_neighbor(neigh_s, new_elem.get()); + } else - new_elems.erase(it); + { + std::vector bcids; + boundary_info.boundary_ids(old_elem, s, bcids); + boundary_info.add_side(new_elem.get(), 1, bcids); + } + + // Set in-cavity neighbors' neighbor pointers + set_neighbors(new_elem.get(), BoundaryInfo::invalid_id); + + // C++ allows function argument evaluation in any + // order, but we need get() to precede move + Elem * new_elem_ptr = new_elem.get(); + new_elems.emplace(new_elem_ptr, std::move(new_elem)); + + check_delaunay_on.insert(new_elem_ptr); } + } - // Everybody found their match? - libmesh_assert(neighbors_CW.empty()); - libmesh_assert(neighbors_CCW.empty()); + boundary_info.remove(old_elem); + } - // Because we're preserving boundaries here, our naive cavity - // triangulation might not be a Delaunay triangulation. Let's - // check and if necessary fix that; we depend on it when doing - // future point insertions. - restore_delaunay(check_delaunay_on, boundary_info); + // Now that we're done using our cavity elems (including with a + // cavity.find() that used a comparator that dereferences the + // elements!) it's safe to delete them. + for (Elem * old_elem : cavity) + { + if (const auto it = new_elems.find(old_elem); it == new_elems.end()) + mesh.delete_elem(old_elem); + else + new_elems.erase(it); + } + + // Everybody found their match? + libmesh_assert(neighbors_CW.empty()); + libmesh_assert(neighbors_CCW.empty()); - // This is too expensive to do on every cavity in devel mode + // Because we're preserving boundaries here, our naive cavity + // triangulation might not be a Delaunay triangulation. Let's + // check and if necessary fix that; we depend on it when doing + // future point insertions. + restore_delaunay(check_delaunay_on, boundary_info); + + // This is too expensive to do on every cavity in devel mode #ifdef DEBUG - libmesh_assert_delaunay(mesh, new_elems); + libmesh_assert_delaunay(mesh, new_elems); #endif - } + } // If we added any new boundary nodes, we're going to need to keep // track of the changes they made to the outer polyline and/or to // any holes. if (!next_boundary_node.empty()) + { + auto checked_emplace = [this](dof_id_type new_first, dof_id_type new_second) { - auto checked_emplace = [this](dof_id_type new_first, - dof_id_type new_second) - { #ifdef DEBUG - for (auto [first, second] : this->segments) - { - libmesh_assert_not_equal_to(first, new_first); - libmesh_assert_not_equal_to(second, new_second); - } - if (!this->segments.empty()) - libmesh_assert_equal_to(this->segments.back().second, new_first); + for (auto [first, second] : this->segments) + { + libmesh_assert_not_equal_to(first, new_first); + libmesh_assert_not_equal_to(second, new_second); + } + if (!this->segments.empty()) + libmesh_assert_equal_to(this->segments.back().second, new_first); #endif - libmesh_assert_not_equal_to(new_first, new_second); + libmesh_assert_not_equal_to(new_first, new_second); - this->segments.emplace_back (new_first, new_second); - }; + this->segments.emplace_back(new_first, new_second); + }; - // Keep track of the outer polyline - if (this->segments.empty()) - { - dof_id_type last_id = DofObject::invalid_id; + // Keep track of the outer polyline + if (this->segments.empty()) + { + dof_id_type last_id = DofObject::invalid_id; - // Custom loop because we increment node_it 1+ times inside - for (auto node_it = _mesh.nodes_begin(), - node_end = _mesh.nodes_end(); - node_it != node_end;) - { - Node & node = **node_it; - ++node_it; + // Custom loop because we increment node_it 1+ times inside + for (auto node_it = _mesh.nodes_begin(), node_end = _mesh.nodes_end(); node_it != node_end;) + { + Node & node = **node_it; + ++node_it; - const dof_id_type node_id = node.id(); + const dof_id_type node_id = node.id(); - // Don't add Steiner points - if (node_id >= _n_boundary_nodes) - break; + // Don't add Steiner points + if (node_id >= _n_boundary_nodes) + break; - // Connect up the previous node, if we didn't already - // connect it after some newly inserted nodes - if (!this->segments.empty()) - last_id = this->segments.back().second; + // Connect up the previous node, if we didn't already + // connect it after some newly inserted nodes + if (!this->segments.empty()) + last_id = this->segments.back().second; - if (last_id != DofObject::invalid_id && - last_id != node_id) - checked_emplace(last_id, node_id); + if (last_id != DofObject::invalid_id && last_id != node_id) + checked_emplace(last_id, node_id); - last_id = node_id; + last_id = node_id; - // Connect to any newly inserted nodes - Node * this_node = &node; - auto it = next_boundary_node.find(*this_node); - while (it != next_boundary_node.end()) - { - libmesh_assert(this_node->valid_id()); - Node * next_node = it->second; - libmesh_assert(next_node->valid_id()); + // Connect to any newly inserted nodes + Node * this_node = &node; + auto it = next_boundary_node.find(*this_node); + while (it != next_boundary_node.end()) + { + libmesh_assert(this_node->valid_id()); + Node * next_node = it->second; + libmesh_assert(next_node->valid_id()); - if (node_it != node_end && - next_node == *node_it) - ++node_it; + if (node_it != node_end && next_node == *node_it) + ++node_it; - checked_emplace(this_node->id(), next_node->id()); + checked_emplace(this_node->id(), next_node->id()); - this_node = next_node; - if (this_node->id() == this->segments.front().first) - break; + this_node = next_node; + if (this_node->id() == this->segments.front().first) + break; - it = next_boundary_node.find(*this_node); - } - } + it = next_boundary_node.find(*this_node); + } + } - // We expect a closed loop here - if (this->segments.back().second != this->segments.front().first) - checked_emplace(this->segments.back().second, - this->segments.front().first); + // We expect a closed loop here + if (this->segments.back().second != this->segments.front().first) + checked_emplace(this->segments.back().second, this->segments.front().first); + } + else + { + std::vector> old_segments; + old_segments.swap(this->segments); + + auto old_it = old_segments.begin(); + + const Node * node = _mesh.node_ptr(old_it->first); + const Node * const first_node = node; + + do + { + const dof_id_type node_id = node->id(); + if (const auto it = next_boundary_node.find(*node); it == next_boundary_node.end()) + { + while (node_id != old_it->first) + { + ++old_it; + libmesh_assert(old_it != old_segments.end()); + } + node = mesh.node_ptr(old_it->second); } - else + else { - std::vector> old_segments; - old_segments.swap(this->segments); + node = it->second; + } - auto old_it = old_segments.begin(); + checked_emplace(node_id, node->id()); + } while (node != first_node); + } - const Node * node = _mesh.node_ptr(old_it->first); - const Node * const first_node = node; + // Keep track of any holes + if (this->_holes) + { + // Do we have any holes that need to be newly replaced? + for (const Hole * hole : *this->_holes) + { + if (this->replaced_holes.count(hole)) + continue; - do - { - const dof_id_type node_id = node->id(); - if (const auto it = next_boundary_node.find(*node); - it == next_boundary_node.end()) - { - while (node_id != old_it->first) - { - ++old_it; - libmesh_assert(old_it != old_segments.end()); - } - node = mesh.node_ptr(old_it->second); - } - else - { - node = it->second; - } - - checked_emplace(node_id, node->id()); - } - while (node != first_node); - } + bool hole_point_insertion = false; + for (auto p : make_range(hole->n_points())) + if (next_boundary_node.count(hole->point(p))) + { + hole_point_insertion = true; + break; + } + if (hole_point_insertion) + this->replaced_holes.emplace(hole, std::make_unique(*hole)); + } - // Keep track of any holes - if (this->_holes) - { - // Do we have any holes that need to be newly replaced? - for (const Hole * hole : *this->_holes) - { - if (this->replaced_holes.count(hole)) - continue; - - bool hole_point_insertion = false; - for (auto p : make_range(hole->n_points())) - if (next_boundary_node.count(hole->point(p))) - { - hole_point_insertion = true; - break; - } - if (hole_point_insertion) - this->replaced_holes.emplace - (hole, std::make_unique(*hole)); - } + // If we have any holes that are being replaced, make sure + // their replacements are up to date. + for (const Hole * hole : *this->_holes) + { + auto hole_it = replaced_holes.find(hole); + if (hole_it == replaced_holes.end()) + continue; - // If we have any holes that are being replaced, make sure - // their replacements are up to date. - for (const Hole * hole : *this->_holes) - { - auto hole_it = replaced_holes.find(hole); - if (hole_it == replaced_holes.end()) - continue; - - ArbitraryHole & arb = *hole_it->second; - - // We only need to update a replacement that's just had - // new points inserted - bool point_inserted = false; - for (const Point & point : arb.get_points()) - if (next_boundary_node.count(point)) - { - point_inserted = true; - break; - } - - if (!point_inserted) - continue; - - // Find all points in the replacement hole - std::vector new_points; - - // Our outer polyline is expected to have points in - // counter-clockwise order, so it proceeds "to the left" - // from the point of view of rays inside the domain - // pointing outward, and our next_boundary_node ordering - // was filled accordingly. - // - // Our inner holes are expected to have points in - // counter-clockwise order, but for holes "to the left - // as viewed from the hole interior is the *opposite* of - // "to the left as viewed from the domain interior". We - // need to build the updated hole ordering "backwards". - - // We should never see duplicate points when we add one - // to a hole; if we do then we did something wrong. - auto push_back_new_point = [&new_points](const Point & p) { - // O(1) assert in devel - libmesh_assert(new_points.empty() || - new_points.back() != p); + ArbitraryHole & arb = *hole_it->second; + + // We only need to update a replacement that's just had + // new points inserted + bool point_inserted = false; + for (const Point & point : arb.get_points()) + if (next_boundary_node.count(point)) + { + point_inserted = true; + break; + } + + if (!point_inserted) + continue; + + // Find all points in the replacement hole + std::vector new_points; + + // Our outer polyline is expected to have points in + // counter-clockwise order, so it proceeds "to the left" + // from the point of view of rays inside the domain + // pointing outward, and our next_boundary_node ordering + // was filled accordingly. + // + // Our inner holes are expected to have points in + // counter-clockwise order, but for holes "to the left + // as viewed from the hole interior is the *opposite* of + // "to the left as viewed from the domain interior". We + // need to build the updated hole ordering "backwards". + + // We should never see duplicate points when we add one + // to a hole; if we do then we did something wrong. + auto push_back_new_point = [&new_points](const Point & p) + { + // O(1) assert in devel + libmesh_assert(new_points.empty() || new_points.back() != p); #ifdef DEBUG - // O(N) asserts in dbg - for (auto old_p : new_points) - libmesh_assert_not_equal_to(old_p, p); + // O(N) asserts in dbg + for (auto old_p : new_points) + libmesh_assert_not_equal_to(old_p, p); #endif - new_points.push_back(p); - }; - - for (auto point_it = arb.get_points().rbegin(), - point_end = arb.get_points().rend(); - point_it != point_end;) - { - Point point = *point_it; - ++point_it; - - if (new_points.empty() || - (point != new_points.back() && - point != new_points.front())) - push_back_new_point(point); - - auto it = next_boundary_node.find(point); - while (it != next_boundary_node.end()) - { - point = *it->second; - if (point == new_points.front()) - break; - if (point_it != point_end && - point == *point_it) - ++point_it; - push_back_new_point(point); - it = next_boundary_node.find(point); - } - } - - std::reverse(new_points.begin(), new_points.end()); - - arb.set_points(std::move(new_points)); - } + new_points.push_back(p); + }; + + for (auto point_it = arb.get_points().rbegin(), point_end = arb.get_points().rend(); + point_it != point_end;) + { + Point point = *point_it; + ++point_it; + + if (new_points.empty() || (point != new_points.back() && point != new_points.front())) + push_back_new_point(point); + + auto it = next_boundary_node.find(point); + while (it != next_boundary_node.end()) + { + point = *it->second; + if (point == new_points.front()) + break; + if (point_it != point_end && point == *point_it) + ++point_it; + push_back_new_point(point); + it = next_boundary_node.find(point); + } } + + std::reverse(new_points.begin(), new_points.end()); + + arb.set_points(std::move(new_points)); + } } + } // Okay, *now* we can add the new elements. for (auto & [raw_elem, unique_elem] : new_elems) - { - libmesh_assert_equal_to(raw_elem, unique_elem.get()); - libmesh_assert(!raw_elem->is_flipped()); - libmesh_ignore(raw_elem); // Old gcc warns "unused variable" - mesh.add_elem(std::move(unique_elem)); - } + { + libmesh_assert_equal_to(raw_elem, unique_elem.get()); + libmesh_assert(!raw_elem->is_flipped()); + libmesh_ignore(raw_elem); // Old gcc warns "unused variable" + mesh.add_elem(std::move(unique_elem)); + } // Did we add anything? return !new_elems.empty(); } - -bool Poly2TriTriangulator::should_refine_elem(Elem & elem) +bool +Poly2TriTriangulator::should_refine_elem(Elem & elem) { const Real min_area_target = this->desired_area(); - FunctionBase *area_func = this->has_auto_area_function() ? this->get_auto_area_function() : this->get_desired_area_function(); + FunctionBase * area_func = this->has_auto_area_function() + ? this->get_auto_area_function() + : this->get_desired_area_function(); // If this isn't a question, why are we here? - libmesh_assert(min_area_target > 0 || - area_func != nullptr || - this->has_auto_area_function()); + libmesh_assert(min_area_target > 0 || area_func != nullptr || this->has_auto_area_function()); const Real area = elem.volume(); @@ -1491,40 +1403,33 @@ bool Poly2TriTriangulator::should_refine_elem(Elem & elem) // decision quickly if (!area_func && !this->has_auto_area_function()) return (area > min_area_target); - else if(area_func && this->has_auto_area_function()) - libmesh_warning("WARNING: both desired are function and automatic area function are set. Using automatic area function."); + else if (area_func && this->has_auto_area_function()) + libmesh_warning("WARNING: both desired are function and automatic area function are set. " + "Using automatic area function."); // If we do? // // See if we're meeting the local area target at all the elem // vertices first for (auto v : make_range(elem.n_vertices())) - { - // If we have an auto area function, we'll use it and override other area options - const Real local_area_target = (*area_func)(elem.point(v)); - libmesh_error_msg_if - (local_area_target <= 0, - "Non-positive desired element areas are unachievable"); - if (area > local_area_target) - return true; - } + { + // If we have an auto area function, we'll use it and override other area options + const Real local_area_target = (*area_func)(elem.point(v)); + libmesh_error_msg_if(local_area_target <= 0, + "Non-positive desired element areas are unachievable"); + if (area > local_area_target) + return true; + } // If our vertices are happy, it's still possible that our interior // isn't. Are we allowed not to bother checking it? if (!min_area_target) return false; - libmesh_not_implemented_msg - ("Combining a minimum desired_area with an area function isn't yet supported."); + libmesh_not_implemented_msg( + "Combining a minimum desired_area with an area function isn't yet supported."); } - } // namespace libMesh - - - - - - #endif // LIBMESH_HAVE_POLY2TRI From aeba92cbf937dc6aebf9259010f2904eeed04278 Mon Sep 17 00:00:00 2001 From: Connor Ethan Ouellette Date: Mon, 15 Jun 2026 14:44:27 -0600 Subject: [PATCH 2/3] Restored original file to prevent clang reformat spam --- src/mesh/poly2tri_triangulator.C | 1743 ++++++++++++++++-------------- 1 file changed, 919 insertions(+), 824 deletions(-) diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index 6bd39ec61fb..e1280711f53 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/src/mesh/poly2tri_triangulator.C @@ -15,6 +15,7 @@ // License along with this library; if not, write to the Free Software // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + #include "libmesh/libmesh_config.h" #ifdef LIBMESH_HAVE_POLY2TRI @@ -52,18 +53,19 @@ struct P2TPointCompare } }; -p2t::Point -to_p2t(const libMesh::Point & p) +p2t::Point to_p2t(const libMesh::Point & p) { #if LIBMESH_DIM > 2 - libmesh_error_msg_if(p(2) != 0, "Poly2TriTriangulator only supports point sets in the XY plane"); + libmesh_error_msg_if + (p(2) != 0, + "Poly2TriTriangulator only supports point sets in the XY plane"); #endif return {double(p(0)), double(p(1))}; } -Real -distance_from_circumcircle(const Elem & elem, const Point & p) +Real distance_from_circumcircle(const Elem & elem, + const Point & p) { libmesh_assert_equal_to(elem.n_vertices(), 3); @@ -74,8 +76,10 @@ distance_from_circumcircle(const Elem & elem, const Point & p) return p_dist - radius; } -bool -in_circumcircle(const Elem & elem, const Point & p, const Real tol = 0) + +bool in_circumcircle(const Elem & elem, + const Point & p, + const Real tol = 0) { return (distance_from_circumcircle(elem, p) < tol); @@ -92,8 +96,11 @@ in_circumcircle(const Elem & elem, const Point & p, const Real tol = 0) */ } + std::pair -can_delaunay_swap(const Elem & elem, unsigned short side, Real tol) +can_delaunay_swap(const Elem & elem, + unsigned short side, + Real tol) { const Elem * neigh = elem.neighbor_ptr(side); if (!neigh) @@ -103,37 +110,42 @@ can_delaunay_swap(const Elem & elem, unsigned short side, Real tol) // What neighbor node does elem not share? for (; nn < 3; ++nn) - { - const Node * neigh_node = neigh->node_ptr(nn); - if (neigh_node == elem.node_ptr(0) || neigh_node == elem.node_ptr(1) || - neigh_node == elem.node_ptr(2)) - continue; - - // Might we need to do a diagonal swap here? Avoid - // undoing a borderline swap. - if (in_circumcircle(elem, *neigh_node, tol)) - break; - } + { + const Node * neigh_node = neigh->node_ptr(nn); + if (neigh_node == elem.node_ptr(0) || + neigh_node == elem.node_ptr(1) || + neigh_node == elem.node_ptr(2)) + continue; + + // Might we need to do a diagonal swap here? Avoid + // undoing a borderline swap. + if (in_circumcircle(elem, *neigh_node, tol)) + break; + } if (nn == 3) return {false, 0}; - const unsigned short n = (side + 2) % 3; - const RealVectorValue right = (elem.point((n + 1) % 3) - elem.point(n)).unit(); - const RealVectorValue mid = (neigh->point(nn) - elem.point(n)).unit(); - const RealVectorValue left = (elem.point((n + 2) % 3) - elem.point(n)).unit(); + const unsigned short n = (side+2)%3; + const RealVectorValue right = + (elem.point((n+1)%3)-elem.point(n)).unit(); + const RealVectorValue mid = + (neigh->point(nn)-elem.point(n)).unit(); + const RealVectorValue left = + (elem.point((n+2)%3)-elem.point(n)).unit(); // If the "middle" vector isn't really in the middle, we can't do a // swap without involving other triangles (or we can't at all if // there's a domain boundary in the way) - if (mid * right < left * right || left * mid < left * right) + if (mid*right < left*right || + left*mid < left*right) return {false, 0}; return {true, nn}; } -[[maybe_unused]] void -libmesh_assert_locally_delaunay(const Elem & elem) + +[[maybe_unused]] void libmesh_assert_locally_delaunay(const Elem & elem) { libmesh_ignore(elem); @@ -146,8 +158,9 @@ libmesh_assert_locally_delaunay(const Elem & elem) } template -inline void -libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), Container & new_elems) +inline +void libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), + Container & new_elems) { libmesh_ignore(new_elems); #ifndef NDEBUG @@ -157,168 +170,174 @@ libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh), Container & new_elems) libmesh_assert_locally_delaunay(*elem); for (auto & [raw_elem, unique_elem] : new_elems) - { - libmesh_ignore(unique_elem); // avoid warnings on old gcc - libmesh_assert_locally_delaunay(*raw_elem); - } + { + libmesh_ignore(unique_elem); // avoid warnings on old gcc + libmesh_assert_locally_delaunay(*raw_elem); + } #endif } + // Restore a triangulation's Delaunay property, starting with a set of // all triangles that might initially not be locally Delaunay with // their neighbors. template -inline void -restore_delaunay(Container & check_delaunay_on, BoundaryInfo & boundary_info) +inline +void restore_delaunay(Container & check_delaunay_on, + BoundaryInfo & boundary_info) { LOG_SCOPE("restore_delaunay()", "Poly2TriTriangulator"); while (!check_delaunay_on.empty()) - { - Elem & elem = **check_delaunay_on.begin(); - check_delaunay_on.erase(&elem); - for (auto s : make_range(elem.n_sides())) { - // Can we make a swap here? With what neighbor, with what - // far node? Use a negative tolerance to avoid swapping - // back and forth. - auto [can_swap, nn] = can_delaunay_swap(elem, s, -TOLERANCE * TOLERANCE); - if (!can_swap) - continue; + Elem & elem = **check_delaunay_on.begin(); + check_delaunay_on.erase(&elem); + for (auto s : make_range(elem.n_sides())) + { + // Can we make a swap here? With what neighbor, with what + // far node? Use a negative tolerance to avoid swapping + // back and forth. + auto [can_swap, nn] = + can_delaunay_swap(elem, s, -TOLERANCE*TOLERANCE); + if (!can_swap) + continue; - Elem * neigh = elem.neighbor_ptr(s); + Elem * neigh = elem.neighbor_ptr(s); - // If we made it here it's time to diagonal swap - const unsigned short n = (s + 2) % 3; + // If we made it here it's time to diagonal swap + const unsigned short n = (s+2)%3; - const std::array nodes{elem.node_ptr(n), - elem.node_ptr((n + 1) % 3), - neigh->node_ptr(nn), - elem.node_ptr((n + 2) % 3)}; + const std::array nodes {elem.node_ptr(n), + elem.node_ptr((n+1)%3), neigh->node_ptr(nn), + elem.node_ptr((n+2)%3)}; - // If we have to swap then either we or any of our neighbors - // might no longer be Delaunay - for (auto ds : make_range(3)) - { - if (elem.neighbor_ptr(ds)) - check_delaunay_on.insert(elem.neighbor_ptr(ds)); - if (neigh->neighbor_ptr(ds)) - check_delaunay_on.insert(neigh->neighbor_ptr(ds)); - } - - // An interior boundary between two newly triangulated - // triangles shouldn't have any bcids - libmesh_assert(!boundary_info.n_boundary_ids(neigh, (nn + 1) % 3)); - libmesh_assert(!boundary_info.n_boundary_ids(&elem, (n + 1) % 3)); - - // The two changing boundary sides might have bcids - std::vector bcids; - boundary_info.boundary_ids(&elem, (n + 2) % 3, bcids); - if (!bcids.empty()) - { - boundary_info.remove_side(&elem, (n + 2) % 3); - boundary_info.add_side(neigh, (nn + 1) % 3, bcids); - } + // If we have to swap then either we or any of our neighbors + // might no longer be Delaunay + for (auto ds : make_range(3)) + { + if (elem.neighbor_ptr(ds)) + check_delaunay_on.insert(elem.neighbor_ptr(ds)); + if (neigh->neighbor_ptr(ds)) + check_delaunay_on.insert(neigh->neighbor_ptr(ds)); + } - boundary_info.boundary_ids(neigh, (nn + 2) % 3, bcids); - if (!bcids.empty()) - { - boundary_info.remove_side(neigh, (nn + 2) % 3); - boundary_info.add_side(&elem, (n + 1) % 3, bcids); - } - - elem.set_node((n + 2) % 3, nodes[2]); - neigh->set_node((nn + 2) % 3, nodes[0]); - - // No need for a temporary array to swap these around, if we - // do it in the right order. - // - // Watch me neigh->neigh - Elem * neighneigh = neigh->neighbor_ptr((nn + 2) % 3); - if (neighneigh) - { - unsigned int snn = neighneigh->which_neighbor_am_i(neigh); - neighneigh->set_neighbor(snn, &elem); - } + // An interior boundary between two newly triangulated + // triangles shouldn't have any bcids + libmesh_assert(!boundary_info.n_boundary_ids(neigh, (nn+1)%3)); + libmesh_assert(!boundary_info.n_boundary_ids(&elem, (n+1)%3)); - Elem * elemoldneigh = elem.neighbor_ptr((n + 2) % 3); - if (elemoldneigh) - { - unsigned int seon = elemoldneigh->which_neighbor_am_i(&elem); - elemoldneigh->set_neighbor(seon, neigh); - } - - elem.set_neighbor((n + 1) % 3, neigh->neighbor_ptr((nn + 2) % 3)); - neigh->set_neighbor((nn + 1) % 3, elem.neighbor_ptr((n + 2) % 3)); - elem.set_neighbor((n + 2) % 3, neigh); - neigh->set_neighbor((nn + 2) % 3, &elem); - - // Start over after this much change, don't just loop to the - // next neighbor - break; + // The two changing boundary sides might have bcids + std::vector bcids; + boundary_info.boundary_ids(&elem, (n+2)%3, bcids); + if (!bcids.empty()) + { + boundary_info.remove_side(&elem, (n+2)%3); + boundary_info.add_side(neigh, (nn+1)%3, bcids); + } + + boundary_info.boundary_ids(neigh, (nn+2)%3, bcids); + if (!bcids.empty()) + { + boundary_info.remove_side(neigh, (nn+2)%3); + boundary_info.add_side(&elem, (n+1)%3, bcids); + } + + elem.set_node((n+2)%3, nodes[2]); + neigh->set_node((nn+2)%3, nodes[0]); + + // No need for a temporary array to swap these around, if we + // do it in the right order. + // + // Watch me neigh->neigh + Elem * neighneigh = neigh->neighbor_ptr((nn+2)%3); + if (neighneigh) + { + unsigned int snn = neighneigh->which_neighbor_am_i(neigh); + neighneigh->set_neighbor(snn, &elem); + } + + Elem * elemoldneigh = elem.neighbor_ptr((n+2)%3); + if (elemoldneigh) + { + unsigned int seon = elemoldneigh->which_neighbor_am_i(&elem); + elemoldneigh->set_neighbor(seon, neigh); + } + + elem.set_neighbor((n+1)%3, neigh->neighbor_ptr((nn+2)%3)); + neigh->set_neighbor((nn+1)%3, elem.neighbor_ptr((n+2)%3)); + elem.set_neighbor((n+2)%3, neigh); + neigh->set_neighbor((nn+2)%3, &elem); + + // Start over after this much change, don't just loop to the + // next neighbor + break; + } } - } } -unsigned int -segment_intersection(const Elem & elem, - Point & source, - const Point & target, - unsigned int source_side) + +unsigned int segment_intersection(const Elem & elem, + Point & source, + const Point & target, + unsigned int source_side) { libmesh_assert_equal_to(elem.dim(), 2); const auto ns = elem.n_sides(); for (auto s : make_range(ns)) - { - // Don't go backwards just because some FP roundoff said to - if (s == source_side) - continue; - - const Point v0 = elem.point(s); - const Point v1 = elem.point((s + 1) % ns); + { + // Don't go backwards just because some FP roundoff said to + if (s == source_side) + continue; - // Calculate intersection parameters (fractions of the distance - // along each segment) - const Real raydx = target(0) - source(0), raydy = target(1) - source(1), edgedx = v1(0) - v0(0), - edgedy = v1(1) - v0(1); - const Real denom = edgedx * raydy - edgedy * raydx; + const Point v0 = elem.point(s); + const Point v1 = elem.point((s+1)%ns); - // divide-by-zero means the segments are parallel - if (denom == 0) - continue; + // Calculate intersection parameters (fractions of the distance + // along each segment) + const Real raydx = target(0)-source(0), + raydy = target(1)-source(1), + edgedx = v1(0)-v0(0), + edgedy = v1(1)-v0(1); + const Real denom = edgedx * raydy - edgedy * raydx; - const Real one_over_denom = 1 / denom; + // divide-by-zero means the segments are parallel + if (denom == 0) + continue; - const Real targetsdx = v1(0) - target(0), targetsdy = v1(1) - target(1); + const Real one_over_denom = 1 / denom; - const Real t_num = targetsdx * raydy - targetsdy * raydx; - const Real t = t_num * one_over_denom; + const Real targetsdx = v1(0)-target(0), + targetsdy = v1(1)-target(1); - if (t < -TOLERANCE * TOLERANCE || t > 1 + TOLERANCE * TOLERANCE) - continue; + const Real t_num = targetsdx * raydy - + targetsdy * raydx; + const Real t = t_num * one_over_denom; - const Real u_num = targetsdx * edgedy - targetsdy * edgedx; - const Real u = u_num * one_over_denom; + if (t < -TOLERANCE*TOLERANCE || t > 1 + TOLERANCE*TOLERANCE) + continue; - if (u < -TOLERANCE * TOLERANCE || u > 1 + TOLERANCE * TOLERANCE) - continue; + const Real u_num = targetsdx * edgedy - targetsdy * edgedx; + const Real u = u_num * one_over_denom; - /* - // Partial workaround for an old poly2tri bug (issue #39): if we - // end up with boundary points that are nearly-collinear but - // infinitesimally concave, p2t::CDT::Triangulate throws a "null - // triangle" exception. So let's try to be infinitesimally - // convex instead. - const Real ray_fraction = (1-u) * (1+TOLERANCE*TOLERANCE); - */ - const Real ray_fraction = (1 - u); + if (u < -TOLERANCE*TOLERANCE || u > 1 + TOLERANCE*TOLERANCE) + continue; - source(0) += raydx * ray_fraction; - source(1) += raydy * ray_fraction; - return s; - } +/* + // Partial workaround for an old poly2tri bug (issue #39): if we + // end up with boundary points that are nearly-collinear but + // infinitesimally concave, p2t::CDT::Triangulate throws a "null + // triangle" exception. So let's try to be infinitesimally + // convex instead. + const Real ray_fraction = (1-u) * (1+TOLERANCE*TOLERANCE); +*/ + const Real ray_fraction = (1-u); + + source(0) += raydx * ray_fraction; + source(1) += raydy * ray_fraction; + return s; + } return libMesh::invalid_uint; } @@ -332,16 +351,20 @@ namespace libMesh // // Constructor -Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh, dof_id_type n_boundary_nodes) - : TriangulatorInterface(mesh), _n_boundary_nodes(n_boundary_nodes), _refine_bdy_allowed(true) +Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh, + dof_id_type n_boundary_nodes) + : TriangulatorInterface(mesh), + _n_boundary_nodes(n_boundary_nodes), + _refine_bdy_allowed(true) { } + Poly2TriTriangulator::~Poly2TriTriangulator() = default; + // Primary function responsible for performing the triangulation -void -Poly2TriTriangulator::triangulate() +void Poly2TriTriangulator::triangulate() { LOG_SCOPE("triangulate()", "Poly2TriTriangulator"); @@ -365,7 +388,9 @@ Poly2TriTriangulator::triangulate() // We won't support quads any time soon, or 1D/3D in this interface // ever. - if (_elem_type != TRI3 && _elem_type != TRI6 && _elem_type != TRI7) + if (_elem_type != TRI3 && + _elem_type != TRI6 && + _elem_type != TRI7) libmesh_not_implemented(); // If we have no explicit segments defined, we may get them from @@ -386,10 +411,11 @@ Poly2TriTriangulator::triangulate() // This is currently done redundantly in parallel; make sure no // processor quits before the others. do - { - libmesh_parallel_only(_mesh.comm()); - this->triangulate_current_points(); - } while (this->insert_refinement_points()); + { + libmesh_parallel_only(_mesh.comm()); + this->triangulate_current_points(); + } + while (this->insert_refinement_points()); libmesh_parallel_only(_mesh.comm()); @@ -416,8 +442,9 @@ Poly2TriTriangulator::triangulate() _mesh.prepare_for_use(); } -void -Poly2TriTriangulator::set_desired_area_function(FunctionBase * desired) + +void Poly2TriTriangulator::set_desired_area_function + (FunctionBase * desired) { if (desired) _desired_area_func = desired->clone(); @@ -425,16 +452,17 @@ Poly2TriTriangulator::set_desired_area_function(FunctionBase * desired) _desired_area_func.reset(); } -FunctionBase * -Poly2TriTriangulator::get_desired_area_function() + +FunctionBase * Poly2TriTriangulator::get_desired_area_function () { return _desired_area_func.get(); } -bool -Poly2TriTriangulator::is_refine_boundary_allowed(const BoundaryInfo & boundary_info, - const Elem & elem, - unsigned int side) + +bool Poly2TriTriangulator::is_refine_boundary_allowed + (const BoundaryInfo & boundary_info, + const Elem & elem, + unsigned int side) { // We should only be calling this on a boundary side libmesh_assert(!elem.neighbor_ptr(side)); @@ -452,14 +480,14 @@ Poly2TriTriangulator::is_refine_boundary_allowed(const BoundaryInfo & boundary_i // side libmesh_assert(this->_holes); - const boundary_id_type hole_num = bcids[0] - 1; + const boundary_id_type hole_num = bcids[0]-1; libmesh_assert_less(hole_num, this->_holes->size()); const Hole * hole = (*this->_holes)[hole_num]; return hole->refine_boundary_allowed(); } -void -Poly2TriTriangulator::triangulate_current_points() + +void Poly2TriTriangulator::triangulate_current_points() { LOG_SCOPE("triangulate_current_points()", "Poly2TriTriangulator"); @@ -478,13 +506,15 @@ Poly2TriTriangulator::triangulate_current_points() std::vector> inner_hole_points(n_holes); dof_id_type nn = _mesh.max_node_id(); - libmesh_error_msg_if(!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!"); + libmesh_error_msg_if + (!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!"); // Unless we're using an explicit segments list, we assume node ids // are contiguous here. if (this->segments.empty()) - libmesh_error_msg_if(_mesh.n_nodes() != nn, - "Poly2TriTriangulator needs contiguous node ids or explicit segments!"); + libmesh_error_msg_if + (_mesh.n_nodes() != nn, + "Poly2TriTriangulator needs contiguous node ids or explicit segments!"); // And if we have more nodes than outer boundary points, the rest // may be interior "Steiner points". We use a set here so we can @@ -495,85 +525,92 @@ Poly2TriTriangulator::triangulate_current_points() // If we were asked to use all mesh nodes as boundary nodes, now's // the time to see how many that is. if (_n_boundary_nodes == DofObject::invalid_id) - { - _n_boundary_nodes = _mesh.n_nodes(); - libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes), - std::distance(_mesh.nodes_begin(), _mesh.nodes_end())); - } + { + _n_boundary_nodes = _mesh.n_nodes(); + libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes), + std::distance(_mesh.nodes_begin(), + _mesh.nodes_end())); + + } else - libmesh_assert_less_equal(_n_boundary_nodes, _mesh.n_nodes()); + libmesh_assert_less_equal(_n_boundary_nodes, + _mesh.n_nodes()); // Prepare poly2tri points for our nodes, sorted into outer boundary // points and interior Steiner points. if (this->segments.empty()) - { - // If we have no segments even after taking elems into account, - // the nodal id ordering defines our outer polyline ordering - for (auto & node : _mesh.node_ptr_range()) { - const p2t::Point pt = to_p2t(*node); + // If we have no segments even after taking elems into account, + // the nodal id ordering defines our outer polyline ordering + for (auto & node : _mesh.node_ptr_range()) + { + const p2t::Point pt = to_p2t(*node); - // If we're out of boundary nodes, the rest are going to be - // Steiner points or hole points - if (node->id() < _n_boundary_nodes) - outer_boundary_points.push_back(pt); - else - steiner_points.insert(pt); + // If we're out of boundary nodes, the rest are going to be + // Steiner points or hole points + if (node->id() < _n_boundary_nodes) + outer_boundary_points.push_back(pt); + else + steiner_points.insert(pt); - // We're not going to support overlapping nodes on the boundary - if (point_node_map.count(pt)) - libmesh_not_implemented(); + // We're not going to support overlapping nodes on the boundary + if (point_node_map.count(pt)) + libmesh_not_implemented(); - point_node_map.emplace(pt, node); + point_node_map.emplace(pt, node); + } } - } // If we have explicit segments defined, that's our outer polyline // ordering: else - { - // Let's make sure our segments are in order - dof_id_type last_id = DofObject::invalid_id; - - // Add nodes from every segment, in order, to the outer polyline - for (auto [segment_start, segment_end] : this->segments) { - if (last_id != DofObject::invalid_id) - libmesh_error_msg_if(segment_start != last_id, "Disconnected triangulator segments"); - last_id = segment_end; + // Let's make sure our segments are in order + dof_id_type last_id = DofObject::invalid_id; - Node * node = _mesh.query_node_ptr(segment_start); + // Add nodes from every segment, in order, to the outer polyline + for (auto [segment_start, segment_end] : this->segments) + { + if (last_id != DofObject::invalid_id) + libmesh_error_msg_if(segment_start != last_id, + "Disconnected triangulator segments"); + last_id = segment_end; - libmesh_error_msg_if(!node, - "Triangulator segments reference nonexistent node id " << segment_start); + Node * node = _mesh.query_node_ptr(segment_start); - outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1))); - p2t::Point * pt = &outer_boundary_points.back(); + libmesh_error_msg_if(!node, + "Triangulator segments reference nonexistent node id " << + segment_start); - // We're not going to support overlapping nodes on the boundary - if (point_node_map.count(*pt)) - libmesh_not_implemented_msg("Triangulating overlapping boundary nodes is unsupported"); + outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1))); + p2t::Point * pt = &outer_boundary_points.back(); - point_node_map.emplace(*pt, node); - } + // We're not going to support overlapping nodes on the boundary + if (point_node_map.count(*pt)) + libmesh_not_implemented_msg + ("Triangulating overlapping boundary nodes is unsupported"); + + point_node_map.emplace(*pt, node); + } - libmesh_error_msg_if(last_id != this->segments[0].first, - "Non-closed-loop triangulator segments"); + libmesh_error_msg_if(last_id != this->segments[0].first, + "Non-closed-loop triangulator segments"); - // If we have points that aren't in any segments, those will be - // Steiner points - for (auto & node : _mesh.node_ptr_range()) - { - const p2t::Point pt = to_p2t(*node); - if (const auto it = point_node_map.find(pt); it == point_node_map.end()) - { - steiner_points.insert(pt); - point_node_map.emplace(pt, node); - } - else - libmesh_assert_equal_to(it->second, node); + // If we have points that aren't in any segments, those will be + // Steiner points + for (auto & node : _mesh.node_ptr_range()) + { + const p2t::Point pt = to_p2t(*node); + if (const auto it = point_node_map.find(pt); + it == point_node_map.end()) + { + steiner_points.insert(pt); + point_node_map.emplace(pt, node); + } + else + libmesh_assert_equal_to(it->second, node); + } } - } // If we have any elements from a previous triangulation, we're // going to replace them with our own. If we have any elements that @@ -585,19 +622,23 @@ Poly2TriTriangulator::triangulate_current_points() // triangle. We'll give the outer boundary BC 0, and give holes ids // starting from 1. We've already got the point_node_map to find // nodes, so we can just key on pairs of node ids to identify a side. - std::unordered_map, boundary_id_type, libMesh::hash> - side_boundary_id; + std::unordered_map, + boundary_id_type, libMesh::hash> side_boundary_id; const boundary_id_type outer_bcid = 0; const std::size_t n_outer = outer_boundary_points.size(); for (auto i : make_range(n_outer)) - { - const Node *node1 = libmesh_map_find(point_node_map, outer_boundary_points[i]), - *node2 = libmesh_map_find(point_node_map, outer_boundary_points[(i + 1) % n_outer]); - - side_boundary_id.emplace(std::make_pair(node1->id(), node2->id()), outer_bcid); - } + { + const Node * node1 = + libmesh_map_find(point_node_map, outer_boundary_points[i]), + * node2 = + libmesh_map_find(point_node_map, outer_boundary_points[(i+1)%n_outer]); + + side_boundary_id.emplace(std::make_pair(node1->id(), + node2->id()), + outer_bcid); + } // Create poly2tri triangulator with our mesh points std::vector outer_boundary_pointers(n_outer); @@ -606,6 +647,7 @@ Poly2TriTriangulator::triangulate_current_points() outer_boundary_pointers.begin(), [](p2t::Point & p) { return &p; }); + // Make sure shims for holes last as long as the CDT does; the // poly2tri headers don't make clear whether or not they're hanging // on to references to these vectors, and it would be reasonable for @@ -616,58 +658,66 @@ Poly2TriTriangulator::triangulate_current_points() // Add any holes for (auto h : make_range(n_holes)) - { - const Hole * initial_hole = (*_holes)[h]; - auto it = replaced_holes.find(initial_hole); - const Hole & our_hole = (it == replaced_holes.end()) ? *initial_hole : *it->second; - auto & poly2tri_hole = inner_hole_points[h]; - - for (auto i : make_range(our_hole.n_points())) { - Point p = our_hole.point(i); - poly2tri_hole.emplace_back(to_p2t(p)); - - const auto & pt = poly2tri_hole.back(); + const Hole * initial_hole = (*_holes)[h]; + auto it = replaced_holes.find(initial_hole); + const Hole & our_hole = + (it == replaced_holes.end()) ? + *initial_hole : *it->second; + auto & poly2tri_hole = inner_hole_points[h]; + + for (auto i : make_range(our_hole.n_points())) + { + Point p = our_hole.point(i); + poly2tri_hole.emplace_back(to_p2t(p)); - // This won't be a steiner point. - steiner_points.erase(pt); + const auto & pt = poly2tri_hole.back(); - // If we see a hole point already in the mesh, we'll share - // that node. This might be a problem if it's a boundary - // node, but it might just be the same hole point already - // added during a previous triangulation refinement step. - if (point_node_map.count(pt)) - { - libmesh_assert_equal_to(point_node_map[pt], _mesh.query_node_ptr(point_node_map[pt]->id())); - } - else - { - Node * node = _mesh.add_point(p, nn++); - point_node_map[pt] = node; - } - } + // This won't be a steiner point. + steiner_points.erase(pt); - const boundary_id_type inner_bcid = h + 1; - const std::size_t n_inner = poly2tri_hole.size(); + // If we see a hole point already in the mesh, we'll share + // that node. This might be a problem if it's a boundary + // node, but it might just be the same hole point already + // added during a previous triangulation refinement step. + if (point_node_map.count(pt)) + { + libmesh_assert_equal_to + (point_node_map[pt], + _mesh.query_node_ptr(point_node_map[pt]->id())); + } + else + { + Node * node = _mesh.add_point(p, nn++); + point_node_map[pt] = node; + } + } - for (auto i : make_range(n_inner)) - { - const Node *node1 = libmesh_map_find(point_node_map, poly2tri_hole[i]), - *node2 = libmesh_map_find(point_node_map, poly2tri_hole[(i + 1) % n_inner]); + const boundary_id_type inner_bcid = h+1; + const std::size_t n_inner = poly2tri_hole.size(); - side_boundary_id.emplace(std::make_pair(node1->id(), node2->id()), inner_bcid); - } + for (auto i : make_range(n_inner)) + { + const Node * node1 = + libmesh_map_find(point_node_map, poly2tri_hole[i]), + * node2 = + libmesh_map_find(point_node_map, poly2tri_hole[(i+1)%n_inner]); + + side_boundary_id.emplace(std::make_pair(node1->id(), + node2->id()), + inner_bcid); + } - auto & poly2tri_ptrs = inner_hole_pointers[h]; - poly2tri_ptrs.resize(n_inner); + auto & poly2tri_ptrs = inner_hole_pointers[h]; + poly2tri_ptrs.resize(n_inner); - std::transform(poly2tri_hole.begin(), - poly2tri_hole.end(), - poly2tri_ptrs.begin(), - [](p2t::Point & p) { return &p; }); + std::transform(poly2tri_hole.begin(), + poly2tri_hole.end(), + poly2tri_ptrs.begin(), + [](p2t::Point & p) { return &p; }); - cdt.AddHole(poly2tri_ptrs); - } + cdt.AddHole(poly2tri_ptrs); + } // Add any steiner points. We had them in a set, but post-C++11 // that won't give us non-const element access (even if we @@ -693,52 +743,60 @@ Poly2TriTriangulator::triangulate_current_points() // Add the triangles to our Mesh data structure. for (auto ptri_ptr : triangles) - { - p2t::Triangle & ptri = *ptri_ptr; - - // We always use TRI3 here, since that's what we have nodes for; - // if we need a higher order we can convert at the end. - auto elem = Elem::build_with_id(TRI3, next_id++); - for (auto v : make_range(3)) { - const p2t::Point & vertex = *ptri.GetPoint(v); + p2t::Triangle & ptri = *ptri_ptr; - Node * node = libmesh_map_find(point_node_map, vertex); - libmesh_assert(node); - elem->set_node(v, node); - } + // We always use TRI3 here, since that's what we have nodes for; + // if we need a higher order we can convert at the end. + auto elem = Elem::build_with_id(TRI3, next_id++); + for (auto v : make_range(3)) + { + const p2t::Point & vertex = *ptri.GetPoint(v); - // We expect a consistent triangle orientation - libmesh_assert(!elem->is_flipped()); + Node * node = libmesh_map_find(point_node_map, vertex); + libmesh_assert(node); + elem->set_node(v, node); + } - Elem * added_elem = _mesh.add_elem(std::move(elem)); + // We expect a consistent triangle orientation + libmesh_assert(!elem->is_flipped()); - for (auto v : make_range(3)) - { - const Node &node1 = added_elem->node_ref(v), &node2 = added_elem->node_ref((v + 1) % 3); + Elem * added_elem = _mesh.add_elem(std::move(elem)); - auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id())); - if (it == side_boundary_id.end()) - it = side_boundary_id.find(std::make_pair(node2.id(), node1.id())); - if (it != side_boundary_id.end()) - boundary_info.add_side(added_elem, v, it->second); + for (auto v : make_range(3)) + { + const Node & node1 = added_elem->node_ref(v), + & node2 = added_elem->node_ref((v+1)%3); + + auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id())); + if (it == side_boundary_id.end()) + it = side_boundary_id.find(std::make_pair(node2.id(), node1.id())); + if (it != side_boundary_id.end()) + boundary_info.add_side(added_elem, v, it->second); + } } - } } -bool -Poly2TriTriangulator::insert_refinement_points() + + +bool Poly2TriTriangulator::insert_refinement_points() { LOG_SCOPE("insert_refinement_points()", "Poly2TriTriangulator"); if (this->minimum_angle() != 0) libmesh_not_implemented(); - BoundaryInfo & boundary_info = _mesh.get_boundary_info(); // We need neighbor pointers for ray casting and cavity finding UnstructuredMesh & mesh = dynamic_cast(this->_mesh); mesh.find_neighbors(); + if (this->desired_area() == 0 && + this->get_desired_area_function() == nullptr && + !this->has_auto_area_function()) + return false; + + BoundaryInfo & boundary_info = _mesh.get_boundary_info(); + // We won't immediately add these, lest we invalidate iterators on a // ReplicatedMesh. They'll still be in the mesh neighbor topology // for the purpose of doing Delaunay cavity stuff, so we need to @@ -756,10 +814,8 @@ Poly2TriTriangulator::insert_refinement_points() // not-yet-added elements so we can't iterate based on proper // element ids ... but we can set a temporary element id to use for // the purpose. - struct cmp - { - bool operator()(Elem * a, Elem * b) const - { + struct cmp { + bool operator()(Elem * a, Elem * b) const { libmesh_assert(a == b || a->id() != b->id()); return (a->id() < b->id()); } @@ -777,17 +833,14 @@ Poly2TriTriangulator::insert_refinement_points() // about the isomorphisms! If a triangle's nodes are permuted on // one processor vs another that's an issue. So sort our input // carefully. - std::set all_elems{mesh.elements_begin(), mesh.elements_end(), comp}; + std::set all_elems + { mesh.elements_begin(), mesh.elements_end(), comp }; restore_delaunay(all_elems, boundary_info); libmesh_assert_delaunay(mesh, new_elems); } - if (this->desired_area() == 0 && this->get_desired_area_function() == nullptr && - !this->has_auto_area_function()) - return false; - // Map of which points follow which in the boundary polylines. If // we have to add new boundary points, we'll use this to construct // an updated this->segments to retriangulate with. If we have to @@ -805,10 +858,10 @@ Poly2TriTriangulator::insert_refinement_points() #ifdef DEBUG std::unordered_set mesh_points; for (const Node * node : mesh.node_ptr_range()) - { - libmesh_assert(!mesh_points.count(*node)); - mesh_points.insert(*node); - } + { + libmesh_assert(!mesh_points.count(*node)); + mesh_points.insert(*node); + } #endif auto add_point = [&mesh, @@ -816,586 +869,621 @@ Poly2TriTriangulator::insert_refinement_points() &mesh_points, #endif &nn](const Point & p) - { + { #ifdef DEBUG - libmesh_assert(!mesh_points.count(p)); - mesh_points.insert(p); + libmesh_assert(!mesh_points.count(p)); + mesh_points.insert(p); #endif - return mesh.add_point(p, nn++); - }; + return mesh.add_point(p, nn++); + }; for (auto & elem : mesh.element_ptr_range()) - { - // element_ptr_range skips deleted elements ... right? - libmesh_assert(elem); - libmesh_assert(elem->valid_id()); + { + // element_ptr_range skips deleted elements ... right? + libmesh_assert(elem); + libmesh_assert(elem->valid_id()); - // We only handle triangles in our triangulation - libmesh_assert_equal_to(elem->level(), 0u); - libmesh_assert_equal_to(elem->type(), TRI3); + // We only handle triangles in our triangulation + libmesh_assert_equal_to(elem->level(), 0u); + libmesh_assert_equal_to(elem->type(), TRI3); - // If this triangle is as small as we desire, move along - if (!should_refine_elem(*elem)) - continue; + // If this triangle is as small as we desire, move along + if (!should_refine_elem(*elem)) + continue; - // Otherwise add a Steiner point. We'd like to add the - // circumcenter ... - Point new_pt = elem->quasicircumcenter(); + // Otherwise add a Steiner point. We'd like to add the + // circumcenter ... + Point new_pt = elem->quasicircumcenter(); - // And to give it a node; - Node * new_node = nullptr; + // And to give it a node; + Node * new_node = nullptr; - // But that might be outside our triangle, or even outside the - // boundary. We'll find a triangle that should contain our new - // point - Elem * cavity_elem = elem; // Start looking at elem anyway + // But that might be outside our triangle, or even outside the + // boundary. We'll find a triangle that should contain our new + // point + Elem * cavity_elem = elem; // Start looking at elem anyway - // We'll refine a boundary later if necessary. - auto boundary_refine = [this, &next_boundary_node, &cavity_elem, &new_node](unsigned int side) - { - libmesh_ignore(this); // Only used in dbg/devel - libmesh_assert(new_node); - libmesh_assert(new_node->valid_id()); - - Node *old_segment_start = cavity_elem->node_ptr(side), - *old_segment_end = cavity_elem->node_ptr((side + 1) % 3); - libmesh_assert(old_segment_start); - libmesh_assert(old_segment_start->valid_id()); - libmesh_assert(old_segment_end); - libmesh_assert(old_segment_end->valid_id()); - - if (auto it = next_boundary_node.find(*old_segment_start); it != next_boundary_node.end()) - { - libmesh_assert(it->second == old_segment_end); - it->second = new_node; - } - else + // We'll refine a boundary later if necessary. + auto boundary_refine = [this, &next_boundary_node, + &cavity_elem, &new_node] + (unsigned int side) { - // This would be an O(N) sanity check if we already - // have a segments vector or any holes. :-P - libmesh_assert(!this->segments.empty() || (_holes && !_holes->empty()) || - (old_segment_end->id() == old_segment_start->id() + 1)); - next_boundary_node[*old_segment_start] = new_node; - } - - next_boundary_node[*new_node] = old_segment_end; - }; - - // Let's find a triangle containing our new point, or at least - // containing the end of a ray leading from our current triangle - // to the new point. - Point ray_start = elem->vertex_average(); + libmesh_ignore(this); // Only used in dbg/devel + libmesh_assert(new_node); + libmesh_assert(new_node->valid_id()); + + Node * old_segment_start = cavity_elem->node_ptr(side), + * old_segment_end = cavity_elem->node_ptr((side+1)%3); + libmesh_assert(old_segment_start); + libmesh_assert(old_segment_start->valid_id()); + libmesh_assert(old_segment_end); + libmesh_assert(old_segment_end->valid_id()); + + if (auto it = next_boundary_node.find(*old_segment_start); + it != next_boundary_node.end()) + { + libmesh_assert(it->second == old_segment_end); + it->second = new_node; + } + else + { + // This would be an O(N) sanity check if we already + // have a segments vector or any holes. :-P + libmesh_assert(!this->segments.empty() || + (_holes && !_holes->empty()) || + (old_segment_end->id() == + old_segment_start->id() + 1)); + next_boundary_node[*old_segment_start] = new_node; + } - // What side are we coming from, and what side are we going to? - unsigned int source_side = invalid_uint; - unsigned int side = invalid_uint; + next_boundary_node[*new_node] = old_segment_end; + }; - while (!cavity_elem->contains_point(new_pt)) - { - side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side); + // Let's find a triangle containing our new point, or at least + // containing the end of a ray leading from our current triangle + // to the new point. + Point ray_start = elem->vertex_average(); - libmesh_assert_not_equal_to(side, invalid_uint); + // What side are we coming from, and what side are we going to? + unsigned int source_side = invalid_uint; + unsigned int side = invalid_uint; - Elem * neigh = cavity_elem->neighbor_ptr(side); - // If we're on a boundary, stop there. Refine the boundary - // if we're allowed, the boundary element otherwise. - if (!neigh) - { - if (this->is_refine_boundary_allowed(boundary_info, *cavity_elem, side)) + while (!cavity_elem->contains_point(new_pt)) { - new_pt = ray_start; - new_node = add_point(new_pt); - boundary_refine(side); - } - else - { - // Should we just add the vertex average of the - // boundary element, to minimize the number of - // slivers created? - // - // new_pt = cavity_elem->vertex_average(); - // - // That works for a while, but it - // seems to be able to "run away" and leave us with - // crazy slivers on boundaries if we push interior - // refinement too far while disabling boundary - // refinement. - // - // Let's go back to refining our original problem - // element. - cavity_elem = elem; - new_pt = cavity_elem->vertex_average(); - new_node = add_point(new_pt); - - // This was going to be a side refinement but it's - // now an internal refinement - side = invalid_uint; - } - - break; - } + side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side); - source_side = neigh->which_neighbor_am_i(cavity_elem); - cavity_elem = neigh; - side = invalid_uint; - } + libmesh_assert_not_equal_to (side, invalid_uint); - // If we're ready to create a new node and we're not on a - // boundary ... should we be? We don't want to create any - // sliver elements or confuse poly2tri or anything. - if (side == invalid_uint && !new_node) - { - unsigned int worst_side = libMesh::invalid_uint; - Real worst_cos = 1; - for (auto s : make_range(3u)) - { - // We never snap to a non-domain-boundary - if (cavity_elem->neighbor_ptr(s)) - continue; + Elem * neigh = cavity_elem->neighbor_ptr(side); + // If we're on a boundary, stop there. Refine the boundary + // if we're allowed, the boundary element otherwise. + if (!neigh) + { + if (this->is_refine_boundary_allowed(boundary_info, + *cavity_elem, + side)) + { + new_pt = ray_start; + new_node = add_point(new_pt); + boundary_refine(side); + } + else + { + // Should we just add the vertex average of the + // boundary element, to minimize the number of + // slivers created? + // + // new_pt = cavity_elem->vertex_average(); + // + // That works for a while, but it + // seems to be able to "run away" and leave us with + // crazy slivers on boundaries if we push interior + // refinement too far while disabling boundary + // refinement. + // + // Let's go back to refining our original problem + // element. + cavity_elem = elem; + new_pt = cavity_elem->vertex_average(); + new_node = add_point(new_pt); + + // This was going to be a side refinement but it's + // now an internal refinement + side = invalid_uint; + } - Real ax = cavity_elem->point(s)(0) - new_pt(0), ay = cavity_elem->point(s)(1) - new_pt(1), - bx = cavity_elem->point((s + 1) % 3)(0) - new_pt(0), - by = cavity_elem->point((s + 1) % 3)(1) - new_pt(1); - const Real my_cos = - (ax * bx + ay * by) / std::sqrt(ax * ax + ay * ay) / std::sqrt(bx * bx + by * by); + break; + } - if (my_cos < worst_cos) - { - worst_side = s; - worst_cos = my_cos; + source_side = neigh->which_neighbor_am_i(cavity_elem); + cavity_elem = neigh; + side = invalid_uint; } - } - - // If we'd create a sliver element on the side, let's just - // refine the side instead, if we're allowed. - if (worst_cos < -0.6) // -0.5 is the best we could enforce? - { - side = worst_side; - if (this->is_refine_boundary_allowed(boundary_info, *cavity_elem, side)) + // If we're ready to create a new node and we're not on a + // boundary ... should we be? We don't want to create any + // sliver elements or confuse poly2tri or anything. + if (side == invalid_uint && !new_node) { - // Let's just try bisecting for now - new_pt = (cavity_elem->point(side) + cavity_elem->point((side + 1) % 3)) / 2; - new_node = add_point(new_pt); - boundary_refine(side); - } - else // Do the best we can under these restrictions - { - new_pt = cavity_elem->vertex_average(); - new_node = add_point(new_pt); + unsigned int worst_side = libMesh::invalid_uint; + Real worst_cos = 1; + for (auto s : make_range(3u)) + { + // We never snap to a non-domain-boundary + if (cavity_elem->neighbor_ptr(s)) + continue; + + Real ax = cavity_elem->point(s)(0) - new_pt(0), + ay = cavity_elem->point(s)(1) - new_pt(1), + bx = cavity_elem->point((s+1)%3)(0) - new_pt(0), + by = cavity_elem->point((s+1)%3)(1) - new_pt(1); + const Real my_cos = (ax*bx+ay*by) / + std::sqrt(ax*ax+ay*ay) / + std::sqrt(bx*bx+by*by); + + if (my_cos < worst_cos) + { + worst_side = s; + worst_cos = my_cos; + } + } - // This was going to be a side refinement but it's - // now an internal refinement - side = invalid_uint; + // If we'd create a sliver element on the side, let's just + // refine the side instead, if we're allowed. + if (worst_cos < -0.6) // -0.5 is the best we could enforce? + { + side = worst_side; + + if (this->is_refine_boundary_allowed(boundary_info, + *cavity_elem, + side)) + { + // Let's just try bisecting for now + new_pt = (cavity_elem->point(side) + + cavity_elem->point((side+1)%3)) / 2; + new_node = add_point(new_pt); + boundary_refine(side); + } + else // Do the best we can under these restrictions + { + new_pt = cavity_elem->vertex_average(); + new_node = add_point(new_pt); + + // This was going to be a side refinement but it's + // now an internal refinement + side = invalid_uint; + } + } + else + new_node = add_point(new_pt); } - } else - new_node = add_point(new_pt); - } - else - libmesh_assert(new_node); + libmesh_assert(new_node); - // Find the Delaunay cavity around the new point. - std::set cavity(comp); + // Find the Delaunay cavity around the new point. + std::set cavity(comp); - std::set unchecked_cavity({cavity_elem}, comp); - while (!unchecked_cavity.empty()) - { - std::set checking_cavity(comp); - checking_cavity.swap(unchecked_cavity); - for (Elem * checking_elem : checking_cavity) - { - for (auto s : make_range(3u)) + std::set unchecked_cavity ({cavity_elem}, comp); + while (!unchecked_cavity.empty()) { - Elem * neigh = checking_elem->neighbor_ptr(s); - if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh)) - continue; + std::set checking_cavity(comp); + checking_cavity.swap(unchecked_cavity); + for (Elem * checking_elem : checking_cavity) + { + for (auto s : make_range(3u)) + { + Elem * neigh = checking_elem->neighbor_ptr(s); + if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh)) + continue; + + if (in_circumcircle(*neigh, new_pt, TOLERANCE*TOLERANCE)) + unchecked_cavity.insert(neigh); + } + } - if (in_circumcircle(*neigh, new_pt, TOLERANCE * TOLERANCE)) - unchecked_cavity.insert(neigh); + libmesh_merge_move(cavity, checking_cavity); } - } - libmesh_merge_move(cavity, checking_cavity); - } - - // Retriangulate the Delaunay cavity. - // Each of our cavity triangle edges that are exterior to the - // cavity will be a source of one new triangle. + // Retriangulate the Delaunay cavity. + // Each of our cavity triangle edges that are exterior to the + // cavity will be a source of one new triangle. - // Set of elements that might need Delaunay swaps - std::set check_delaunay_on(comp); + // Set of elements that might need Delaunay swaps + std::set check_delaunay_on(comp); - // Keep maps for doing neighbor pointer assignment. Not going - // to iterate through these so hashing pointers is fine. - std::unordered_map> neighbors_CCW, neighbors_CW; + // Keep maps for doing neighbor pointer assignment. Not going + // to iterate through these so hashing pointers is fine. + std::unordered_map> + neighbors_CCW, neighbors_CW; - for (Elem * old_elem : cavity) - { - for (auto s : make_range(3u)) - { - Elem * neigh = old_elem->neighbor_ptr(s); - if (!neigh || !cavity.count(neigh)) + for (Elem * old_elem : cavity) { - Node *node_CW = old_elem->node_ptr(s), *node_CCW = old_elem->node_ptr((s + 1) % 3); - - auto set_neighbors = [&neighbors_CW, &neighbors_CCW, &node_CW, &node_CCW, &boundary_info]( - Elem * new_neigh, boundary_id_type bcid) - { - // Set clockwise neighbor and vice-versa if possible - if (const auto CW_it = neighbors_CW.find(node_CW); CW_it == neighbors_CW.end()) - { - libmesh_assert(!neighbors_CCW.count(node_CW)); - neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid); - } - else + for (auto s : make_range(3u)) { - Elem * neigh_CW = CW_it->second.first; - if (new_neigh) - { - new_neigh->set_neighbor(0, neigh_CW); - boundary_id_type bcid_CW = CW_it->second.second; - if (bcid_CW != BoundaryInfo::invalid_id) - boundary_info.add_side(new_neigh, 0, bcid_CW); - } - if (neigh_CW) - { - neigh_CW->set_neighbor(2, new_neigh); - if (bcid != BoundaryInfo::invalid_id) - boundary_info.add_side(neigh_CW, 2, bcid); - } - neighbors_CW.erase(CW_it); + Elem * neigh = old_elem->neighbor_ptr(s); + if (!neigh || !cavity.count(neigh)) + { + Node * node_CW = old_elem->node_ptr(s), + * node_CCW = old_elem->node_ptr((s+1)%3); + + auto set_neighbors = + [&neighbors_CW, &neighbors_CCW, &node_CW, + &node_CCW, &boundary_info] + (Elem * new_neigh, boundary_id_type bcid) + { + // Set clockwise neighbor and vice-versa if possible + if (const auto CW_it = neighbors_CW.find(node_CW); + CW_it == neighbors_CW.end()) + { + libmesh_assert(!neighbors_CCW.count(node_CW)); + neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid); + } + else + { + Elem * neigh_CW = CW_it->second.first; + if (new_neigh) + { + new_neigh->set_neighbor(0, neigh_CW); + boundary_id_type bcid_CW = CW_it->second.second; + if (bcid_CW != BoundaryInfo::invalid_id) + boundary_info.add_side(new_neigh, 0, bcid_CW); + + } + if (neigh_CW) + { + neigh_CW->set_neighbor(2, new_neigh); + if (bcid != BoundaryInfo::invalid_id) + boundary_info.add_side(neigh_CW, 2, bcid); + } + neighbors_CW.erase(CW_it); + } + + // Set counter-CW neighbor and vice-versa if possible + if (const auto CCW_it = neighbors_CCW.find(node_CCW); + CCW_it == neighbors_CCW.end()) + { + libmesh_assert(!neighbors_CW.count(node_CCW)); + neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid); + } + else + { + Elem * neigh_CCW = CCW_it->second.first; + if (new_neigh) + { + boundary_id_type bcid_CCW = CCW_it->second.second; + new_neigh->set_neighbor(2, neigh_CCW); + if (bcid_CCW != BoundaryInfo::invalid_id) + boundary_info.add_side(new_neigh, 2, bcid_CCW); + } + if (neigh_CCW) + { + neigh_CCW->set_neighbor(0, new_neigh); + if (bcid != BoundaryInfo::invalid_id) + boundary_info.add_side(neigh_CCW, 0, bcid); + } + neighbors_CCW.erase(CCW_it); + } + }; + + // We aren't going to try to add a sliver element if we + // have a new boundary node here. We do need to + // keep track of other elements' neighbors, though. + if (old_elem == cavity_elem && + s == side) + { + std::vector bcids; + boundary_info.boundary_ids(old_elem, s, bcids); + libmesh_assert_equal_to(bcids.size(), 1); + set_neighbors(nullptr, bcids[0]); + continue; + } + + auto new_elem = Elem::build_with_id(TRI3, ne++); + new_elem->set_node(0, new_node); + new_elem->set_node(1, node_CW); + new_elem->set_node(2, node_CCW); + libmesh_assert(!new_elem->is_flipped()); + + // Set in-and-out-of-cavity neighbor pointers + new_elem->set_neighbor(1, neigh); + if (neigh) + { + const unsigned int neigh_s = + neigh->which_neighbor_am_i(old_elem); + neigh->set_neighbor(neigh_s, new_elem.get()); + } + else + { + std::vector bcids; + boundary_info.boundary_ids(old_elem, s, bcids); + boundary_info.add_side(new_elem.get(), 1, bcids); + } + + // Set in-cavity neighbors' neighbor pointers + set_neighbors(new_elem.get(), BoundaryInfo::invalid_id); + + // C++ allows function argument evaluation in any + // order, but we need get() to precede move + Elem * new_elem_ptr = new_elem.get(); + new_elems.emplace(new_elem_ptr, std::move(new_elem)); + + check_delaunay_on.insert(new_elem_ptr); + } } - // Set counter-CW neighbor and vice-versa if possible - if (const auto CCW_it = neighbors_CCW.find(node_CCW); CCW_it == neighbors_CCW.end()) - { - libmesh_assert(!neighbors_CW.count(node_CCW)); - neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid); - } - else - { - Elem * neigh_CCW = CCW_it->second.first; - if (new_neigh) - { - boundary_id_type bcid_CCW = CCW_it->second.second; - new_neigh->set_neighbor(2, neigh_CCW); - if (bcid_CCW != BoundaryInfo::invalid_id) - boundary_info.add_side(new_neigh, 2, bcid_CCW); - } - if (neigh_CCW) - { - neigh_CCW->set_neighbor(0, new_neigh); - if (bcid != BoundaryInfo::invalid_id) - boundary_info.add_side(neigh_CCW, 0, bcid); - } - neighbors_CCW.erase(CCW_it); - } - }; - - // We aren't going to try to add a sliver element if we - // have a new boundary node here. We do need to - // keep track of other elements' neighbors, though. - if (old_elem == cavity_elem && s == side) - { - std::vector bcids; - boundary_info.boundary_ids(old_elem, s, bcids); - libmesh_assert_equal_to(bcids.size(), 1); - set_neighbors(nullptr, bcids[0]); - continue; - } - - auto new_elem = Elem::build_with_id(TRI3, ne++); - new_elem->set_node(0, new_node); - new_elem->set_node(1, node_CW); - new_elem->set_node(2, node_CCW); - libmesh_assert(!new_elem->is_flipped()); + boundary_info.remove(old_elem); + } - // Set in-and-out-of-cavity neighbor pointers - new_elem->set_neighbor(1, neigh); - if (neigh) - { - const unsigned int neigh_s = neigh->which_neighbor_am_i(old_elem); - neigh->set_neighbor(neigh_s, new_elem.get()); - } + // Now that we're done using our cavity elems (including with a + // cavity.find() that used a comparator that dereferences the + // elements!) it's safe to delete them. + for (Elem * old_elem : cavity) + { + if (const auto it = new_elems.find(old_elem); + it == new_elems.end()) + mesh.delete_elem(old_elem); else - { - std::vector bcids; - boundary_info.boundary_ids(old_elem, s, bcids); - boundary_info.add_side(new_elem.get(), 1, bcids); - } - - // Set in-cavity neighbors' neighbor pointers - set_neighbors(new_elem.get(), BoundaryInfo::invalid_id); - - // C++ allows function argument evaluation in any - // order, but we need get() to precede move - Elem * new_elem_ptr = new_elem.get(); - new_elems.emplace(new_elem_ptr, std::move(new_elem)); - - check_delaunay_on.insert(new_elem_ptr); + new_elems.erase(it); } - } - boundary_info.remove(old_elem); - } + // Everybody found their match? + libmesh_assert(neighbors_CW.empty()); + libmesh_assert(neighbors_CCW.empty()); - // Now that we're done using our cavity elems (including with a - // cavity.find() that used a comparator that dereferences the - // elements!) it's safe to delete them. - for (Elem * old_elem : cavity) - { - if (const auto it = new_elems.find(old_elem); it == new_elems.end()) - mesh.delete_elem(old_elem); - else - new_elems.erase(it); - } - - // Everybody found their match? - libmesh_assert(neighbors_CW.empty()); - libmesh_assert(neighbors_CCW.empty()); + // Because we're preserving boundaries here, our naive cavity + // triangulation might not be a Delaunay triangulation. Let's + // check and if necessary fix that; we depend on it when doing + // future point insertions. + restore_delaunay(check_delaunay_on, boundary_info); - // Because we're preserving boundaries here, our naive cavity - // triangulation might not be a Delaunay triangulation. Let's - // check and if necessary fix that; we depend on it when doing - // future point insertions. - restore_delaunay(check_delaunay_on, boundary_info); - - // This is too expensive to do on every cavity in devel mode + // This is too expensive to do on every cavity in devel mode #ifdef DEBUG - libmesh_assert_delaunay(mesh, new_elems); + libmesh_assert_delaunay(mesh, new_elems); #endif - } + } // If we added any new boundary nodes, we're going to need to keep // track of the changes they made to the outer polyline and/or to // any holes. if (!next_boundary_node.empty()) - { - auto checked_emplace = [this](dof_id_type new_first, dof_id_type new_second) { -#ifdef DEBUG - for (auto [first, second] : this->segments) + auto checked_emplace = [this](dof_id_type new_first, + dof_id_type new_second) { - libmesh_assert_not_equal_to(first, new_first); - libmesh_assert_not_equal_to(second, new_second); - } - if (!this->segments.empty()) - libmesh_assert_equal_to(this->segments.back().second, new_first); +#ifdef DEBUG + for (auto [first, second] : this->segments) + { + libmesh_assert_not_equal_to(first, new_first); + libmesh_assert_not_equal_to(second, new_second); + } + if (!this->segments.empty()) + libmesh_assert_equal_to(this->segments.back().second, new_first); #endif - libmesh_assert_not_equal_to(new_first, new_second); + libmesh_assert_not_equal_to(new_first, new_second); - this->segments.emplace_back(new_first, new_second); - }; + this->segments.emplace_back (new_first, new_second); + }; - // Keep track of the outer polyline - if (this->segments.empty()) - { - dof_id_type last_id = DofObject::invalid_id; - - // Custom loop because we increment node_it 1+ times inside - for (auto node_it = _mesh.nodes_begin(), node_end = _mesh.nodes_end(); node_it != node_end;) - { - Node & node = **node_it; - ++node_it; - - const dof_id_type node_id = node.id(); - - // Don't add Steiner points - if (node_id >= _n_boundary_nodes) - break; + // Keep track of the outer polyline + if (this->segments.empty()) + { + dof_id_type last_id = DofObject::invalid_id; - // Connect up the previous node, if we didn't already - // connect it after some newly inserted nodes - if (!this->segments.empty()) - last_id = this->segments.back().second; + // Custom loop because we increment node_it 1+ times inside + for (auto node_it = _mesh.nodes_begin(), + node_end = _mesh.nodes_end(); + node_it != node_end;) + { + Node & node = **node_it; + ++node_it; - if (last_id != DofObject::invalid_id && last_id != node_id) - checked_emplace(last_id, node_id); + const dof_id_type node_id = node.id(); - last_id = node_id; + // Don't add Steiner points + if (node_id >= _n_boundary_nodes) + break; - // Connect to any newly inserted nodes - Node * this_node = &node; - auto it = next_boundary_node.find(*this_node); - while (it != next_boundary_node.end()) - { - libmesh_assert(this_node->valid_id()); - Node * next_node = it->second; - libmesh_assert(next_node->valid_id()); + // Connect up the previous node, if we didn't already + // connect it after some newly inserted nodes + if (!this->segments.empty()) + last_id = this->segments.back().second; - if (node_it != node_end && next_node == *node_it) - ++node_it; + if (last_id != DofObject::invalid_id && + last_id != node_id) + checked_emplace(last_id, node_id); - checked_emplace(this_node->id(), next_node->id()); + last_id = node_id; - this_node = next_node; - if (this_node->id() == this->segments.front().first) - break; + // Connect to any newly inserted nodes + Node * this_node = &node; + auto it = next_boundary_node.find(*this_node); + while (it != next_boundary_node.end()) + { + libmesh_assert(this_node->valid_id()); + Node * next_node = it->second; + libmesh_assert(next_node->valid_id()); - it = next_boundary_node.find(*this_node); - } - } + if (node_it != node_end && + next_node == *node_it) + ++node_it; - // We expect a closed loop here - if (this->segments.back().second != this->segments.front().first) - checked_emplace(this->segments.back().second, this->segments.front().first); - } - else - { - std::vector> old_segments; - old_segments.swap(this->segments); + checked_emplace(this_node->id(), next_node->id()); - auto old_it = old_segments.begin(); + this_node = next_node; + if (this_node->id() == this->segments.front().first) + break; - const Node * node = _mesh.node_ptr(old_it->first); - const Node * const first_node = node; + it = next_boundary_node.find(*this_node); + } + } - do - { - const dof_id_type node_id = node->id(); - if (const auto it = next_boundary_node.find(*node); it == next_boundary_node.end()) - { - while (node_id != old_it->first) - { - ++old_it; - libmesh_assert(old_it != old_segments.end()); - } - node = mesh.node_ptr(old_it->second); + // We expect a closed loop here + if (this->segments.back().second != this->segments.front().first) + checked_emplace(this->segments.back().second, + this->segments.front().first); } - else + else { - node = it->second; - } - - checked_emplace(node_id, node->id()); - } while (node != first_node); - } - - // Keep track of any holes - if (this->_holes) - { - // Do we have any holes that need to be newly replaced? - for (const Hole * hole : *this->_holes) - { - if (this->replaced_holes.count(hole)) - continue; - - bool hole_point_insertion = false; - for (auto p : make_range(hole->n_points())) - if (next_boundary_node.count(hole->point(p))) - { - hole_point_insertion = true; - break; - } - if (hole_point_insertion) - this->replaced_holes.emplace(hole, std::make_unique(*hole)); - } + std::vector> old_segments; + old_segments.swap(this->segments); - // If we have any holes that are being replaced, make sure - // their replacements are up to date. - for (const Hole * hole : *this->_holes) - { - auto hole_it = replaced_holes.find(hole); - if (hole_it == replaced_holes.end()) - continue; + auto old_it = old_segments.begin(); - ArbitraryHole & arb = *hole_it->second; + const Node * node = _mesh.node_ptr(old_it->first); + const Node * const first_node = node; - // We only need to update a replacement that's just had - // new points inserted - bool point_inserted = false; - for (const Point & point : arb.get_points()) - if (next_boundary_node.count(point)) - { - point_inserted = true; - break; - } + do + { + const dof_id_type node_id = node->id(); + if (const auto it = next_boundary_node.find(*node); + it == next_boundary_node.end()) + { + while (node_id != old_it->first) + { + ++old_it; + libmesh_assert(old_it != old_segments.end()); + } + node = mesh.node_ptr(old_it->second); + } + else + { + node = it->second; + } + + checked_emplace(node_id, node->id()); + } + while (node != first_node); + } - if (!point_inserted) - continue; - - // Find all points in the replacement hole - std::vector new_points; - - // Our outer polyline is expected to have points in - // counter-clockwise order, so it proceeds "to the left" - // from the point of view of rays inside the domain - // pointing outward, and our next_boundary_node ordering - // was filled accordingly. - // - // Our inner holes are expected to have points in - // counter-clockwise order, but for holes "to the left - // as viewed from the hole interior is the *opposite* of - // "to the left as viewed from the domain interior". We - // need to build the updated hole ordering "backwards". - - // We should never see duplicate points when we add one - // to a hole; if we do then we did something wrong. - auto push_back_new_point = [&new_points](const Point & p) + // Keep track of any holes + if (this->_holes) { - // O(1) assert in devel - libmesh_assert(new_points.empty() || new_points.back() != p); + // Do we have any holes that need to be newly replaced? + for (const Hole * hole : *this->_holes) + { + if (this->replaced_holes.count(hole)) + continue; + + bool hole_point_insertion = false; + for (auto p : make_range(hole->n_points())) + if (next_boundary_node.count(hole->point(p))) + { + hole_point_insertion = true; + break; + } + if (hole_point_insertion) + this->replaced_holes.emplace + (hole, std::make_unique(*hole)); + } + + // If we have any holes that are being replaced, make sure + // their replacements are up to date. + for (const Hole * hole : *this->_holes) + { + auto hole_it = replaced_holes.find(hole); + if (hole_it == replaced_holes.end()) + continue; + + ArbitraryHole & arb = *hole_it->second; + + // We only need to update a replacement that's just had + // new points inserted + bool point_inserted = false; + for (const Point & point : arb.get_points()) + if (next_boundary_node.count(point)) + { + point_inserted = true; + break; + } + + if (!point_inserted) + continue; + + // Find all points in the replacement hole + std::vector new_points; + + // Our outer polyline is expected to have points in + // counter-clockwise order, so it proceeds "to the left" + // from the point of view of rays inside the domain + // pointing outward, and our next_boundary_node ordering + // was filled accordingly. + // + // Our inner holes are expected to have points in + // counter-clockwise order, but for holes "to the left + // as viewed from the hole interior is the *opposite* of + // "to the left as viewed from the domain interior". We + // need to build the updated hole ordering "backwards". + + // We should never see duplicate points when we add one + // to a hole; if we do then we did something wrong. + auto push_back_new_point = [&new_points](const Point & p) { + // O(1) assert in devel + libmesh_assert(new_points.empty() || + new_points.back() != p); #ifdef DEBUG - // O(N) asserts in dbg - for (auto old_p : new_points) - libmesh_assert_not_equal_to(old_p, p); + // O(N) asserts in dbg + for (auto old_p : new_points) + libmesh_assert_not_equal_to(old_p, p); #endif - new_points.push_back(p); - }; - - for (auto point_it = arb.get_points().rbegin(), point_end = arb.get_points().rend(); - point_it != point_end;) - { - Point point = *point_it; - ++point_it; - - if (new_points.empty() || (point != new_points.back() && point != new_points.front())) - push_back_new_point(point); - - auto it = next_boundary_node.find(point); - while (it != next_boundary_node.end()) - { - point = *it->second; - if (point == new_points.front()) - break; - if (point_it != point_end && point == *point_it) - ++point_it; - push_back_new_point(point); - it = next_boundary_node.find(point); - } + new_points.push_back(p); + }; + + for (auto point_it = arb.get_points().rbegin(), + point_end = arb.get_points().rend(); + point_it != point_end;) + { + Point point = *point_it; + ++point_it; + + if (new_points.empty() || + (point != new_points.back() && + point != new_points.front())) + push_back_new_point(point); + + auto it = next_boundary_node.find(point); + while (it != next_boundary_node.end()) + { + point = *it->second; + if (point == new_points.front()) + break; + if (point_it != point_end && + point == *point_it) + ++point_it; + push_back_new_point(point); + it = next_boundary_node.find(point); + } + } + + std::reverse(new_points.begin(), new_points.end()); + + arb.set_points(std::move(new_points)); + } } - - std::reverse(new_points.begin(), new_points.end()); - - arb.set_points(std::move(new_points)); - } } - } // Okay, *now* we can add the new elements. for (auto & [raw_elem, unique_elem] : new_elems) - { - libmesh_assert_equal_to(raw_elem, unique_elem.get()); - libmesh_assert(!raw_elem->is_flipped()); - libmesh_ignore(raw_elem); // Old gcc warns "unused variable" - mesh.add_elem(std::move(unique_elem)); - } + { + libmesh_assert_equal_to(raw_elem, unique_elem.get()); + libmesh_assert(!raw_elem->is_flipped()); + libmesh_ignore(raw_elem); // Old gcc warns "unused variable" + mesh.add_elem(std::move(unique_elem)); + } // Did we add anything? return !new_elems.empty(); } -bool -Poly2TriTriangulator::should_refine_elem(Elem & elem) + +bool Poly2TriTriangulator::should_refine_elem(Elem & elem) { const Real min_area_target = this->desired_area(); - FunctionBase * area_func = this->has_auto_area_function() - ? this->get_auto_area_function() - : this->get_desired_area_function(); + FunctionBase *area_func = this->has_auto_area_function() ? this->get_auto_area_function() : this->get_desired_area_function(); // If this isn't a question, why are we here? - libmesh_assert(min_area_target > 0 || area_func != nullptr || this->has_auto_area_function()); + libmesh_assert(min_area_target > 0 || + area_func != nullptr || + this->has_auto_area_function()); const Real area = elem.volume(); @@ -1403,33 +1491,40 @@ Poly2TriTriangulator::should_refine_elem(Elem & elem) // decision quickly if (!area_func && !this->has_auto_area_function()) return (area > min_area_target); - else if (area_func && this->has_auto_area_function()) - libmesh_warning("WARNING: both desired are function and automatic area function are set. " - "Using automatic area function."); + else if(area_func && this->has_auto_area_function()) + libmesh_warning("WARNING: both desired are function and automatic area function are set. Using automatic area function."); // If we do? // // See if we're meeting the local area target at all the elem // vertices first for (auto v : make_range(elem.n_vertices())) - { - // If we have an auto area function, we'll use it and override other area options - const Real local_area_target = (*area_func)(elem.point(v)); - libmesh_error_msg_if(local_area_target <= 0, - "Non-positive desired element areas are unachievable"); - if (area > local_area_target) - return true; - } + { + // If we have an auto area function, we'll use it and override other area options + const Real local_area_target = (*area_func)(elem.point(v)); + libmesh_error_msg_if + (local_area_target <= 0, + "Non-positive desired element areas are unachievable"); + if (area > local_area_target) + return true; + } // If our vertices are happy, it's still possible that our interior // isn't. Are we allowed not to bother checking it? if (!min_area_target) return false; - libmesh_not_implemented_msg( - "Combining a minimum desired_area with an area function isn't yet supported."); + libmesh_not_implemented_msg + ("Combining a minimum desired_area with an area function isn't yet supported."); } + } // namespace libMesh + + + + + + #endif // LIBMESH_HAVE_POLY2TRI From 078a19fc2e57e6f8db50b825843d8e1618139bca Mon Sep 17 00:00:00 2001 From: Connor Ethan Ouellette Date: Mon, 15 Jun 2026 14:50:53 -0600 Subject: [PATCH 3/3] Restored edit, without clang reformatting spam --- src/mesh/poly2tri_triangulator.C | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mesh/poly2tri_triangulator.C b/src/mesh/poly2tri_triangulator.C index e1280711f53..1734a3df09f 100644 --- a/src/mesh/poly2tri_triangulator.C +++ b/src/mesh/poly2tri_triangulator.C @@ -790,11 +790,6 @@ bool Poly2TriTriangulator::insert_refinement_points() UnstructuredMesh & mesh = dynamic_cast(this->_mesh); mesh.find_neighbors(); - if (this->desired_area() == 0 && - this->get_desired_area_function() == nullptr && - !this->has_auto_area_function()) - return false; - BoundaryInfo & boundary_info = _mesh.get_boundary_info(); // We won't immediately add these, lest we invalidate iterators on a @@ -841,6 +836,11 @@ bool Poly2TriTriangulator::insert_refinement_points() libmesh_assert_delaunay(mesh, new_elems); } + if (this->desired_area() == 0 && + this->get_desired_area_function() == nullptr && + !this->has_auto_area_function()) + return false; + // Map of which points follow which in the boundary polylines. If // we have to add new boundary points, we'll use this to construct // an updated this->segments to retriangulate with. If we have to