Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 : #include "libmesh/libmesh_config.h"
20 :
21 : #ifdef LIBMESH_HAVE_POLY2TRI
22 :
23 : // We'd like reproduceability here even when different FP rounding can
24 : // lead to different triangle splitting decisions
25 : #include "libmesh/enforce_ieee754.h"
26 :
27 : // libmesh includes
28 : #include "libmesh/poly2tri_triangulator.h"
29 :
30 : #include "libmesh/boundary_info.h"
31 : #include "libmesh/elem.h"
32 : #include "libmesh/enum_elem_type.h"
33 : #include "libmesh/function_base.h"
34 : #include "libmesh/hashing.h"
35 : #include "libmesh/libmesh_logging.h"
36 : #include "libmesh/mesh_serializer.h"
37 : #include "libmesh/mesh_smoother_laplace.h"
38 : #include "libmesh/mesh_triangle_holes.h"
39 : #include "libmesh/unstructured_mesh.h"
40 : #include "libmesh/utility.h"
41 :
42 : // poly2tri includes
43 : #include "libmesh/ignore_warnings.h" // utf-8 comments should be fine...
44 : #include "poly2tri/poly2tri.h"
45 : #include "libmesh/restore_warnings.h"
46 :
47 : // Anonymous namespace - poly2tri doesn't define operator<(Point,Point)
48 : namespace
49 : {
50 : using namespace libMesh;
51 :
52 : struct P2TPointCompare
53 : {
54 1122820 : bool operator()(const p2t::Point & a, const p2t::Point & b) const
55 : {
56 39680804 : return a.x < b.x || (a.x == b.x && a.y < b.y);
57 : }
58 : };
59 :
60 548333 : p2t::Point to_p2t(const libMesh::Point & p)
61 : {
62 : #if LIBMESH_DIM > 2
63 548333 : libmesh_error_msg_if
64 : (p(2) != 0,
65 : "Poly2TriTriangulator only supports point sets in the XY plane");
66 : #endif
67 :
68 548333 : return {double(p(0)), double(p(1))};
69 : }
70 :
71 5608600 : Real distance_from_circumcircle(const Elem & elem,
72 : const Point & p)
73 : {
74 1970920 : libmesh_assert_equal_to(elem.n_vertices(), 3);
75 :
76 5608600 : const Point circumcenter = elem.quasicircumcenter();
77 5608600 : const Real radius = (elem.point(0) - circumcenter).norm();
78 5503160 : const Real p_dist = (p - circumcenter).norm();
79 :
80 5608600 : return p_dist - radius;
81 : }
82 :
83 :
84 1970920 : bool in_circumcircle(const Elem & elem,
85 : const Point & p,
86 : const Real tol = 0)
87 : {
88 2506701 : return (distance_from_circumcircle(elem, p) < tol);
89 :
90 : /*
91 : libmesh_assert_equal_to(elem.n_vertices(), 3);
92 :
93 : const Point pv0 = elem.point(0) - p;
94 : const Point pv1 = elem.point(1) - p;
95 : const Point pv2 = elem.point(2) - p;
96 :
97 : return ((pv0.norm_sq() * (pv1(0)*pv2(1)-pv2(0)*pv1(1))) -
98 : (pv1.norm_sq() * (pv0(0)*pv2(1)-pv2(0)*pv0(1))) +
99 : (pv2.norm_sq() * (pv0(0)*pv1(1)-pv1(0)*pv0(1)))) > 0;
100 : */
101 : }
102 :
103 :
104 : std::pair<bool, unsigned short>
105 5535979 : can_delaunay_swap(const Elem & elem,
106 : unsigned short side,
107 : Real tol)
108 : {
109 5535979 : const Elem * neigh = elem.neighbor_ptr(side);
110 5535979 : if (!neigh)
111 383412 : return {false, 0};
112 :
113 1958074 : unsigned short nn = 0;
114 :
115 : // What neighbor node does elem not share?
116 20609984 : for (; nn < 3; ++nn)
117 : {
118 6151996 : const Node * neigh_node = neigh->node_ptr(nn);
119 13364278 : if (neigh_node == elem.node_ptr(0) ||
120 24031671 : neigh_node == elem.node_ptr(1) ||
121 3255868 : neigh_node == elem.node_ptr(2))
122 10119808 : continue;
123 :
124 : // Might we need to do a diagonal swap here? Avoid
125 : // undoing a borderline swap.
126 5152567 : if (in_circumcircle(elem, *neigh_node, tol))
127 4 : break;
128 : }
129 :
130 5152567 : if (nn == 3)
131 5152425 : return {false, 0};
132 :
133 142 : const unsigned short n = (side+2)%3;
134 : const RealVectorValue right =
135 150 : (elem.point((n+1)%3)-elem.point(n)).unit();
136 : const RealVectorValue mid =
137 150 : (neigh->point(nn)-elem.point(n)).unit();
138 : const RealVectorValue left =
139 150 : (elem.point((n+2)%3)-elem.point(n)).unit();
140 :
141 : // If the "middle" vector isn't really in the middle, we can't do a
142 : // swap without involving other triangles (or we can't at all if
143 : // there's a domain boundary in the way)
144 146 : if (mid*right < left*right ||
145 4 : left*mid < left*right)
146 0 : return {false, 0};
147 :
148 142 : return {true, nn};
149 : }
150 :
151 :
152 660928 : [[maybe_unused]] void libmesh_assert_locally_delaunay(const Elem & elem)
153 : {
154 660928 : libmesh_ignore(elem);
155 :
156 : #ifndef NDEBUG
157 : // -TOLERANCE, because we're fine with something a little inside the
158 : // circumcircle
159 2643712 : for (auto s : make_range(elem.n_sides()))
160 1982784 : libmesh_assert(!can_delaunay_swap(elem, s, -TOLERANCE).first);
161 : #endif
162 660928 : }
163 :
164 : template <typename Container>
165 : inline
166 2356 : void libmesh_assert_delaunay(MeshBase & libmesh_dbg_var(mesh),
167 : Container & new_elems)
168 : {
169 2356 : libmesh_ignore(new_elems);
170 : #ifndef NDEBUG
171 4712 : LOG_SCOPE("libmesh_assert_delaunay()", "Poly2TriTriangulator");
172 :
173 436916 : for (auto & elem : mesh.element_ptr_range())
174 434560 : libmesh_assert_locally_delaunay(*elem);
175 :
176 228724 : for (auto & [raw_elem, unique_elem] : new_elems)
177 : {
178 226368 : libmesh_ignore(unique_elem); // avoid warnings on old gcc
179 226368 : libmesh_assert_locally_delaunay(*raw_elem);
180 : }
181 : #endif
182 2356 : }
183 :
184 :
185 : // Restore a triangulation's Delaunay property, starting with a set of
186 : // all triangles that might initially not be locally Delaunay with
187 : // their neighbors.
188 : template <typename Container>
189 : inline
190 83638 : void restore_delaunay(Container & check_delaunay_on,
191 : BoundaryInfo & boundary_info)
192 : {
193 4712 : LOG_SCOPE("restore_delaunay()", "Poly2TriTriangulator");
194 :
195 1346986 : while (!check_delaunay_on.empty())
196 : {
197 1184422 : Elem & elem = **check_delaunay_on.begin();
198 1184422 : check_delaunay_on.erase(&elem);
199 4737475 : for (auto s : make_range(elem.n_sides()))
200 : {
201 : // Can we make a swap here? With what neighbor, with what
202 : // far node? Use a negative tolerance to avoid swapping
203 : // back and forth.
204 100090 : auto [can_swap, nn] =
205 3553195 : can_delaunay_swap(elem, s, -TOLERANCE*TOLERANCE);
206 3553195 : if (!can_swap)
207 3553053 : continue;
208 :
209 8 : Elem * neigh = elem.neighbor_ptr(s);
210 :
211 : // If we made it here it's time to diagonal swap
212 142 : const unsigned short n = (s+2)%3;
213 :
214 8 : const std::array<Node *,4> nodes {elem.node_ptr(n),
215 142 : elem.node_ptr((n+1)%3), neigh->node_ptr(nn),
216 142 : elem.node_ptr((n+2)%3)};
217 :
218 : // If we have to swap then either we or any of our neighbors
219 : // might no longer be Delaunay
220 568 : for (auto ds : make_range(3))
221 : {
222 438 : if (elem.neighbor_ptr(ds))
223 355 : check_delaunay_on.insert(elem.neighbor_ptr(ds));
224 438 : if (neigh->neighbor_ptr(ds))
225 355 : check_delaunay_on.insert(neigh->neighbor_ptr(ds));
226 : }
227 :
228 : // An interior boundary between two newly triangulated
229 : // triangles shouldn't have any bcids
230 4 : libmesh_assert(!boundary_info.n_boundary_ids(neigh, (nn+1)%3));
231 4 : libmesh_assert(!boundary_info.n_boundary_ids(&elem, (n+1)%3));
232 :
233 : // The two changing boundary sides might have bcids
234 8 : std::vector<boundary_id_type> bcids;
235 142 : boundary_info.boundary_ids(&elem, (n+2)%3, bcids);
236 142 : if (!bcids.empty())
237 : {
238 71 : boundary_info.remove_side(&elem, (n+2)%3);
239 71 : boundary_info.add_side(neigh, (nn+1)%3, bcids);
240 : }
241 :
242 142 : boundary_info.boundary_ids(neigh, (nn+2)%3, bcids);
243 142 : if (!bcids.empty())
244 : {
245 71 : boundary_info.remove_side(neigh, (nn+2)%3);
246 71 : boundary_info.add_side(&elem, (n+1)%3, bcids);
247 : }
248 :
249 142 : elem.set_node((n+2)%3, nodes[2]);
250 142 : neigh->set_node((nn+2)%3, nodes[0]);
251 :
252 : // No need for a temporary array to swap these around, if we
253 : // do it in the right order.
254 : //
255 : // Watch me neigh->neigh
256 8 : Elem * neighneigh = neigh->neighbor_ptr((nn+2)%3);
257 142 : if (neighneigh)
258 : {
259 71 : unsigned int snn = neighneigh->which_neighbor_am_i(neigh);
260 4 : neighneigh->set_neighbor(snn, &elem);
261 : }
262 :
263 8 : Elem * elemoldneigh = elem.neighbor_ptr((n+2)%3);
264 142 : if (elemoldneigh)
265 : {
266 71 : unsigned int seon = elemoldneigh->which_neighbor_am_i(&elem);
267 4 : elemoldneigh->set_neighbor(seon, neigh);
268 : }
269 :
270 12 : elem.set_neighbor((n+1)%3, neigh->neighbor_ptr((nn+2)%3));
271 142 : neigh->set_neighbor((nn+1)%3, elem.neighbor_ptr((n+2)%3));
272 4 : elem.set_neighbor((n+2)%3, neigh);
273 4 : neigh->set_neighbor((nn+2)%3, &elem);
274 :
275 : // Start over after this much change, don't just loop to the
276 : // next neighbor
277 4 : break;
278 : }
279 : }
280 83638 : }
281 :
282 :
283 16188 : unsigned int segment_intersection(const Elem & elem,
284 : Point & source,
285 : const Point & target,
286 : unsigned int source_side)
287 : {
288 456 : libmesh_assert_equal_to(elem.dim(), 2);
289 :
290 16188 : const auto ns = elem.n_sides();
291 :
292 34861 : for (auto s : make_range(ns))
293 : {
294 : // Don't go backwards just because some FP roundoff said to
295 34861 : if (s == source_side)
296 18147 : continue;
297 :
298 34435 : const Point v0 = elem.point(s);
299 34435 : const Point v1 = elem.point((s+1)%ns);
300 :
301 : // Calculate intersection parameters (fractions of the distance
302 : // along each segment)
303 34435 : const Real raydx = target(0)-source(0),
304 34435 : raydy = target(1)-source(1),
305 34435 : edgedx = v1(0)-v0(0),
306 34435 : edgedy = v1(1)-v0(1);
307 34435 : const Real denom = edgedx * raydy - edgedy * raydx;
308 :
309 : // divide-by-zero means the segments are parallel
310 34435 : if (denom == 0)
311 0 : continue;
312 :
313 34435 : const Real one_over_denom = 1 / denom;
314 :
315 34435 : const Real targetsdx = v1(0)-target(0),
316 34435 : targetsdy = v1(1)-target(1);
317 :
318 36375 : const Real t_num = targetsdx * raydy -
319 34435 : targetsdy * raydx;
320 34435 : const Real t = t_num * one_over_denom;
321 :
322 34435 : if (t < -TOLERANCE*TOLERANCE || t > 1 + TOLERANCE*TOLERANCE)
323 6624 : continue;
324 :
325 27619 : const Real u_num = targetsdx * edgedy - targetsdy * edgedx;
326 27619 : const Real u = u_num * one_over_denom;
327 :
328 27619 : if (u < -TOLERANCE*TOLERANCE || u > 1 + TOLERANCE*TOLERANCE)
329 11109 : continue;
330 :
331 : /*
332 : // Partial workaround for an old poly2tri bug (issue #39): if we
333 : // end up with boundary points that are nearly-collinear but
334 : // infinitesimally concave, p2t::CDT::Triangulate throws a "null
335 : // triangle" exception. So let's try to be infinitesimally
336 : // convex instead.
337 : const Real ray_fraction = (1-u) * (1+TOLERANCE*TOLERANCE);
338 : */
339 16188 : const Real ray_fraction = (1-u);
340 :
341 16188 : source(0) += raydx * ray_fraction;
342 16188 : source(1) += raydy * ray_fraction;
343 15732 : return s;
344 : }
345 :
346 0 : return libMesh::invalid_uint;
347 : }
348 :
349 : }
350 :
351 : namespace libMesh
352 : {
353 : //
354 : // Function definitions for the Poly2TriTriangulator class
355 : //
356 :
357 : // Constructor
358 2059 : Poly2TriTriangulator::Poly2TriTriangulator(UnstructuredMesh & mesh,
359 2059 : dof_id_type n_boundary_nodes)
360 : : TriangulatorInterface(mesh),
361 1943 : _n_boundary_nodes(n_boundary_nodes),
362 2059 : _refine_bdy_allowed(true)
363 : {
364 2059 : }
365 :
366 :
367 3886 : Poly2TriTriangulator::~Poly2TriTriangulator() = default;
368 :
369 :
370 : // Primary function responsible for performing the triangulation
371 1988 : void Poly2TriTriangulator::triangulate()
372 : {
373 112 : LOG_SCOPE("triangulate()", "Poly2TriTriangulator");
374 :
375 : // We only operate on serialized meshes. And it's not safe to
376 : // serialize earlier, because it would then be possible for the user
377 : // to re-parallelize the mesh in between there and here.
378 2092 : MeshSerializer serializer(_mesh);
379 :
380 : // We don't yet support every set of Triangulator options in the
381 : // poly2tri implementation
382 :
383 : // We don't support convex hull triangulation, only triangulation of
384 : // (implicitly defined, by node ordering) polygons (with holes if
385 : // requested)
386 1988 : if (_triangulation_type != PSLG)
387 0 : libmesh_not_implemented();
388 :
389 : // We currently don't handle region specifications
390 1988 : if (_regions)
391 0 : libmesh_not_implemented();
392 :
393 : // We won't support quads any time soon, or 1D/3D in this interface
394 : // ever.
395 2044 : if (_elem_type != TRI3 &&
396 1940 : _elem_type != TRI6 &&
397 2 : _elem_type != TRI7)
398 0 : libmesh_not_implemented();
399 :
400 : // If we have no explicit segments defined, we may get them from
401 : // mesh elements
402 1988 : this->elems_to_segments();
403 :
404 : // If we *still* have no explicit segments defined, we get them from
405 : // the order of nodes.
406 1775 : this->nodes_to_segments(_n_boundary_nodes);
407 :
408 : // Insert additional new points in between existing boundary points,
409 : // if that is requested and reasonable
410 1775 : this->insert_any_extra_boundary_points();
411 :
412 : // Triangulate the points we have, then see if we need to add more;
413 : // repeat until we don't need to add more.
414 : //
415 : // This is currently done redundantly in parallel; make sure no
416 : // processor quits before the others.
417 130 : do
418 : {
419 180 : libmesh_parallel_only(_mesh.comm());
420 6390 : this->triangulate_current_points();
421 : }
422 6390 : while (this->insert_refinement_points());
423 :
424 50 : libmesh_parallel_only(_mesh.comm());
425 :
426 : // Okay, we really do need to support boundary ids soon, but we
427 : // don't yet
428 1775 : if (_markers)
429 0 : libmesh_not_implemented();
430 :
431 1775 : _mesh.set_mesh_dimension(2);
432 :
433 : // To the naked eye, a few smoothing iterations usually looks better,
434 : // so we do this by default unless the user says not to.
435 1775 : if (this->_smooth_after_generating)
436 8 : LaplaceMeshSmoother(_mesh, 2).smooth();
437 :
438 : // The user might have requested TRI6 or higher instead of TRI3. We
439 : // can do this before prepare_for_use() because all we need for it
440 : // is find_neighbors(), which we did in insert_refinement_points()
441 1775 : this->increase_triangle_order();
442 :
443 : // Prepare the mesh for use before returning. This ensures (among
444 : // other things) that it is partitioned and therefore users can
445 : // iterate over local elements, etc.
446 1704 : _mesh.prepare_for_use();
447 1972 : }
448 :
449 :
450 710 : void Poly2TriTriangulator::set_desired_area_function
451 : (FunctionBase<Real> * desired)
452 : {
453 710 : if (desired)
454 280 : _desired_area_func = desired->clone();
455 : else
456 16 : _desired_area_func.reset();
457 710 : }
458 :
459 :
460 763960 : FunctionBase<Real> * Poly2TriTriangulator::get_desired_area_function ()
461 : {
462 763960 : return _desired_area_func.get();
463 : }
464 :
465 :
466 22365 : bool Poly2TriTriangulator::is_refine_boundary_allowed
467 : (const BoundaryInfo & boundary_info,
468 : const Elem & elem,
469 : unsigned int side)
470 : {
471 : // We should only be calling this on a boundary side
472 630 : libmesh_assert(!elem.neighbor_ptr(side));
473 :
474 1260 : std::vector<boundary_id_type> bcids;
475 22365 : boundary_info.boundary_ids(&elem, side, bcids);
476 :
477 : // We should have one bcid on every boundary side.
478 630 : libmesh_assert_equal_to(bcids.size(), 1);
479 :
480 22365 : if (bcids[0] == 0)
481 19525 : return this->refine_boundary_allowed();
482 :
483 : // If we're not on an outer boundary side we'd better be on a hole
484 : // side
485 80 : libmesh_assert(this->_holes);
486 :
487 2840 : const boundary_id_type hole_num = bcids[0]-1;
488 80 : libmesh_assert_less(hole_num, this->_holes->size());
489 2840 : const Hole * hole = (*this->_holes)[hole_num];
490 2840 : return hole->refine_boundary_allowed();
491 : }
492 :
493 :
494 6390 : void Poly2TriTriangulator::triangulate_current_points()
495 : {
496 360 : LOG_SCOPE("triangulate_current_points()", "Poly2TriTriangulator");
497 :
498 : // Will the triangulation have holes?
499 6390 : const std::size_t n_holes = _holes != nullptr ? _holes->size() : 0;
500 :
501 : // Mapping from Poly2Tri points to libMesh nodes, so we can get the
502 : // connectivity translated back later.
503 360 : std::map<const p2t::Point, Node *, P2TPointCompare> point_node_map;
504 :
505 : // Poly2Tri data structures
506 : // Poly2Tri takes vectors of pointers-to-Point for some reason, but
507 : // we'll just make those shims to vectors of Point rather than
508 : // individually/manually heap allocating everything.
509 540 : std::vector<p2t::Point> outer_boundary_points;
510 6750 : std::vector<std::vector<p2t::Point>> inner_hole_points(n_holes);
511 :
512 6390 : dof_id_type nn = _mesh.max_node_id();
513 6390 : libmesh_error_msg_if
514 : (!nn, "Poly2TriTriangulator cannot triangulate an empty mesh!");
515 :
516 : // Unless we're using an explicit segments list, we assume node ids
517 : // are contiguous here.
518 6390 : if (this->segments.empty())
519 0 : libmesh_error_msg_if
520 : (_mesh.n_nodes() != nn,
521 : "Poly2TriTriangulator needs contiguous node ids or explicit segments!");
522 :
523 : // And if we have more nodes than outer boundary points, the rest
524 : // may be interior "Steiner points". We use a set here so we can
525 : // cheaply search and erase from it later, when we're identifying
526 : // hole points.
527 360 : std::set<p2t::Point, P2TPointCompare> steiner_points;
528 :
529 : // If we were asked to use all mesh nodes as boundary nodes, now's
530 : // the time to see how many that is.
531 6390 : if (_n_boundary_nodes == DofObject::invalid_id)
532 : {
533 1775 : _n_boundary_nodes = _mesh.n_nodes();
534 50 : libmesh_assert_equal_to(std::ptrdiff_t(_n_boundary_nodes),
535 : std::distance(_mesh.nodes_begin(),
536 : _mesh.nodes_end()));
537 :
538 : }
539 : else
540 130 : libmesh_assert_less_equal(_n_boundary_nodes,
541 : _mesh.n_nodes());
542 :
543 : // Prepare poly2tri points for our nodes, sorted into outer boundary
544 : // points and interior Steiner points.
545 :
546 6390 : if (this->segments.empty())
547 : {
548 : // If we have no segments even after taking elems into account,
549 : // the nodal id ordering defines our outer polyline ordering
550 0 : for (auto & node : _mesh.node_ptr_range())
551 : {
552 0 : const p2t::Point pt = to_p2t(*node);
553 :
554 : // If we're out of boundary nodes, the rest are going to be
555 : // Steiner points or hole points
556 0 : if (node->id() < _n_boundary_nodes)
557 0 : outer_boundary_points.push_back(pt);
558 : else
559 0 : steiner_points.insert(pt);
560 :
561 : // We're not going to support overlapping nodes on the boundary
562 0 : if (point_node_map.count(pt))
563 0 : libmesh_not_implemented();
564 :
565 0 : point_node_map.emplace(pt, node);
566 0 : }
567 : }
568 : // If we have explicit segments defined, that's our outer polyline
569 : // ordering:
570 : else
571 : {
572 : // Let's make sure our segments are in order
573 180 : dof_id_type last_id = DofObject::invalid_id;
574 :
575 : // Add nodes from every segment, in order, to the outer polyline
576 142852 : for (auto [segment_start, segment_end] : this->segments)
577 : {
578 136462 : if (last_id != DofObject::invalid_id)
579 130072 : libmesh_error_msg_if(segment_start != last_id,
580 : "Disconnected triangulator segments");
581 130696 : last_id = segment_end;
582 :
583 136462 : Node * node = _mesh.query_node_ptr(segment_start);
584 :
585 136462 : libmesh_error_msg_if(!node,
586 : "Triangulator segments reference nonexistent node id " <<
587 : segment_start);
588 :
589 136462 : outer_boundary_points.emplace_back(double((*node)(0)), double((*node)(1)));
590 3844 : p2t::Point * pt = &outer_boundary_points.back();
591 :
592 : // We're not going to support overlapping nodes on the boundary
593 3844 : if (point_node_map.count(*pt))
594 0 : libmesh_not_implemented_msg
595 : ("Triangulating overlapping boundary nodes is unsupported");
596 :
597 132618 : point_node_map.emplace(*pt, node);
598 : }
599 :
600 6390 : libmesh_error_msg_if(last_id != this->segments[0].first,
601 : "Non-closed-loop triangulator segments");
602 :
603 : // If we have points that aren't in any segments, those will be
604 : // Steiner points
605 987432 : for (auto & node : _mesh.node_ptr_range())
606 : {
607 515672 : const p2t::Point pt = to_p2t(*node);
608 501544 : if (const auto it = point_node_map.find(pt);
609 14128 : it == point_node_map.end())
610 : {
611 354798 : steiner_points.insert(pt);
612 354798 : point_node_map.emplace(pt, node);
613 : }
614 : else
615 3844 : libmesh_assert_equal_to(it->second, node);
616 6030 : }
617 : }
618 :
619 : // If we have any elements from a previous triangulation, we're
620 : // going to replace them with our own. If we have any elements that
621 : // were used to create our segments, we're done creating and we no
622 : // longer need them.
623 6390 : _mesh.clear_elems();
624 :
625 : // Keep track of what boundary ids we want to assign to each new
626 : // triangle. We'll give the outer boundary BC 0, and give holes ids
627 : // starting from 1. We've already got the point_node_map to find
628 : // nodes, so we can just key on pairs of node ids to identify a side.
629 : std::unordered_map<std::pair<dof_id_type,dof_id_type>,
630 360 : boundary_id_type, libMesh::hash> side_boundary_id;
631 :
632 6390 : const boundary_id_type outer_bcid = 0;
633 360 : const std::size_t n_outer = outer_boundary_points.size();
634 :
635 142852 : for (auto i : make_range(n_outer))
636 : {
637 : const Node * node1 =
638 140306 : libmesh_map_find(point_node_map, outer_boundary_points[i]),
639 : * node2 =
640 140306 : libmesh_map_find(point_node_map, outer_boundary_points[(i+1)%n_outer]);
641 :
642 132618 : side_boundary_id.emplace(std::make_pair(node1->id(),
643 11532 : node2->id()),
644 3844 : outer_bcid);
645 : }
646 :
647 : // Create poly2tri triangulator with our mesh points
648 6570 : std::vector<p2t::Point *> outer_boundary_pointers(n_outer);
649 : std::transform(outer_boundary_points.begin(),
650 : outer_boundary_points.end(),
651 : outer_boundary_pointers.begin(),
652 4024 : [](p2t::Point & p) { return &p; });
653 :
654 :
655 : // Make sure shims for holes last as long as the CDT does; the
656 : // poly2tri headers don't make clear whether or not they're hanging
657 : // on to references to these vectors, and it would be reasonable for
658 : // them to do so.
659 6750 : std::vector<std::vector<p2t::Point *>> inner_hole_pointers(n_holes);
660 :
661 6750 : p2t::CDT cdt{outer_boundary_pointers};
662 :
663 : // Add any holes
664 12780 : for (auto h : make_range(n_holes))
665 : {
666 6570 : const Hole * initial_hole = (*_holes)[h];
667 180 : auto it = replaced_holes.find(initial_hole);
668 : const Hole & our_hole =
669 6432 : (it == replaced_holes.end()) ?
670 42 : *initial_hole : *it->second;
671 360 : auto & poly2tri_hole = inner_hole_points[h];
672 :
673 53179 : for (auto i : make_range(our_hole.n_points()))
674 : {
675 46789 : Point p = our_hole.point(i);
676 92260 : poly2tri_hole.emplace_back(to_p2t(p));
677 :
678 1318 : const auto & pt = poly2tri_hole.back();
679 :
680 : // This won't be a steiner point.
681 1318 : steiner_points.erase(pt);
682 :
683 : // If we see a hole point already in the mesh, we'll share
684 : // that node. This might be a problem if it's a boundary
685 : // node, but it might just be the same hole point already
686 : // added during a previous triangulation refinement step.
687 1318 : if (point_node_map.count(pt))
688 : {
689 1166 : libmesh_assert_equal_to
690 : (point_node_map[pt],
691 : _mesh.query_node_ptr(point_node_map[pt]->id()));
692 : }
693 : else
694 : {
695 5396 : Node * node = _mesh.add_point(p, nn++);
696 5396 : point_node_map[pt] = node;
697 : }
698 : }
699 :
700 6390 : const boundary_id_type inner_bcid = h+1;
701 360 : const std::size_t n_inner = poly2tri_hole.size();
702 :
703 53179 : for (auto i : make_range(n_inner))
704 : {
705 : const Node * node1 =
706 48107 : libmesh_map_find(point_node_map, poly2tri_hole[i]),
707 : * node2 =
708 48107 : libmesh_map_find(point_node_map, poly2tri_hole[(i+1)%n_inner]);
709 :
710 45471 : side_boundary_id.emplace(std::make_pair(node1->id(),
711 3954 : node2->id()),
712 1318 : inner_bcid);
713 : }
714 :
715 360 : auto & poly2tri_ptrs = inner_hole_pointers[h];
716 6390 : poly2tri_ptrs.resize(n_inner);
717 :
718 : std::transform(poly2tri_hole.begin(),
719 : poly2tri_hole.end(),
720 : poly2tri_ptrs.begin(),
721 1498 : [](p2t::Point & p) { return &p; });
722 :
723 6390 : cdt.AddHole(poly2tri_ptrs);
724 : }
725 :
726 : // Add any steiner points. We had them in a set, but post-C++11
727 : // that won't give us non-const element access (even if we
728 : // pinky-promise not to change the elements in any way that affects
729 : // our Comparator), and Poly2Tri wants non-const elements (to store
730 : // edge data?), so we need to move them here.
731 6750 : std::vector<p2t::Point> steiner_vector(steiner_points.begin(), steiner_points.end());
732 180 : steiner_points.clear();
733 330079 : for (auto & p : steiner_vector)
734 323689 : cdt.AddPoint(&p);
735 :
736 : // Triangulate!
737 6390 : cdt.Triangulate();
738 :
739 : // Get poly2tri triangles, turn them into libMesh triangles
740 6570 : std::vector<p2t::Triangle *> triangles = cdt.GetTriangles();
741 :
742 : // Do our own numbering, even on DistributedMesh
743 180 : dof_id_type next_id = 0;
744 :
745 6390 : BoundaryInfo & boundary_info = _mesh.get_boundary_info();
746 6390 : boundary_info.clear();
747 :
748 : // Add the triangles to our Mesh data structure.
749 837019 : for (auto ptri_ptr : triangles)
750 : {
751 23398 : p2t::Triangle & ptri = *ptri_ptr;
752 :
753 : // We always use TRI3 here, since that's what we have nodes for;
754 : // if we need a higher order we can convert at the end.
755 854027 : auto elem = Elem::build_with_id(TRI3, next_id++);
756 3322516 : for (auto v : make_range(3))
757 : {
758 70194 : const p2t::Point & vertex = *ptri.GetPoint(v);
759 :
760 2491887 : Node * node = libmesh_map_find(point_node_map, vertex);
761 70194 : libmesh_assert(node);
762 2491887 : elem->set_node(v, node);
763 : }
764 :
765 : // We expect a consistent triangle orientation
766 23398 : libmesh_assert(!elem->is_flipped());
767 :
768 877425 : Elem * added_elem = _mesh.add_elem(std::move(elem));
769 :
770 3322516 : for (auto v : make_range(3))
771 : {
772 140388 : const Node & node1 = added_elem->node_ref(v),
773 2491887 : & node2 = added_elem->node_ref((v+1)%3);
774 :
775 2491887 : auto it = side_boundary_id.find(std::make_pair(node1.id(), node2.id()));
776 2491887 : if (it == side_boundary_id.end())
777 2355425 : it = side_boundary_id.find(std::make_pair(node2.id(), node1.id()));
778 2491887 : if (it != side_boundary_id.end())
779 183251 : boundary_info.add_side(added_elem, v, it->second);
780 : }
781 783833 : }
782 18450 : }
783 :
784 :
785 :
786 6390 : bool Poly2TriTriangulator::insert_refinement_points()
787 : {
788 360 : LOG_SCOPE("insert_refinement_points()", "Poly2TriTriangulator");
789 :
790 6390 : if (this->minimum_angle() != 0)
791 0 : libmesh_not_implemented();
792 :
793 : // We need neighbor pointers for ray casting and cavity finding
794 6390 : UnstructuredMesh & mesh = dynamic_cast<UnstructuredMesh &>(this->_mesh);
795 6390 : mesh.find_neighbors();
796 :
797 1974 : if (this->desired_area() == 0 &&
798 6461 : this->get_desired_area_function() == nullptr &&
799 2 : !this->has_auto_area_function())
800 2 : return false;
801 :
802 6319 : BoundaryInfo & boundary_info = _mesh.get_boundary_info();
803 :
804 : // We won't immediately add these, lest we invalidate iterators on a
805 : // ReplicatedMesh. They'll still be in the mesh neighbor topology
806 : // for the purpose of doing Delaunay cavity stuff, so we need to
807 : // manage memory here, but there's no point in adding them to the
808 : // Mesh just to remove them again afterward when we hit up poly2tri.
809 :
810 : // We'll need to be able to remove new elems from new_elems, in
811 : // cases where a later refinement insertion has a not-yet-added
812 : // element in its cavity, so we'll use a map here to make searching
813 : // possible.
814 : //
815 : // For parallel consistency, we can't order a container we plan to
816 : // iterate through based on Elem * or a hash of it. We'll be doing
817 : // Delaunay swaps so we can't iterate based on geometry. These are
818 : // not-yet-added elements so we can't iterate based on proper
819 : // element ids ... but we can set a temporary element id to use for
820 : // the purpose.
821 : struct cmp {
822 687618 : bool operator()(Elem * a, Elem * b) const {
823 687618 : libmesh_assert(a == b || a->id() != b->id());
824 1949161 : return (a->id() < b->id());
825 : }
826 : } comp;
827 :
828 356 : std::map<Elem *, std::unique_ptr<Elem>, decltype(comp)> new_elems(comp);
829 :
830 : // We should already be Delaunay when we get here, otherwise we
831 : // won't be able to stay Delaunay later. But we're *not* always
832 : // Delaunay when we get here? What the hell, poly2tri? Fixing this
833 : // is expensive!
834 : {
835 : // restore_delaunay should get to the same Delaunay triangulation up to
836 : // isomorphism regardless of ordering ... but we actually care
837 : // about the isomorphisms! If a triangle's nodes are permuted on
838 : // one processor vs another that's an issue. So sort our input
839 : // carefully.
840 : std::set<Elem *, decltype(comp)> all_elems
841 6853 : { mesh.elements_begin(), mesh.elements_end(), comp };
842 :
843 6319 : restore_delaunay(all_elems, boundary_info);
844 :
845 178 : libmesh_assert_delaunay(mesh, new_elems);
846 : }
847 :
848 : // Map of which points follow which in the boundary polylines. If
849 : // we have to add new boundary points, we'll use this to construct
850 : // an updated this->segments to retriangulate with. If we have to
851 : // add new hole points, we'll use this to insert points into an
852 : // ArbitraryHole.
853 356 : std::unordered_map<Point, Node *> next_boundary_node;
854 :
855 : // In cases where we've been working with contiguous node id ranges;
856 : // let's keep it that way.
857 6319 : dof_id_type nn = _mesh.max_node_id();
858 6319 : dof_id_type ne = _mesh.max_elem_id();
859 :
860 : // We can't handle duplicated nodes. We shouldn't ever create one,
861 : // but let's make sure of that.
862 : #ifdef DEBUG
863 178 : std::unordered_set<Point> mesh_points;
864 14450 : for (const Node * node : mesh.node_ptr_range())
865 : {
866 14272 : libmesh_assert(!mesh_points.count(*node));
867 14272 : mesh_points.insert(*node);
868 : }
869 : #endif
870 :
871 72963 : auto add_point = [&mesh,
872 : #ifdef DEBUG
873 : &mesh_points,
874 : #endif
875 86031 : &nn](const Point & p)
876 : {
877 : #ifdef DEBUG
878 2178 : libmesh_assert(!mesh_points.count(p));
879 2178 : mesh_points.insert(p);
880 : #endif
881 76721 : return mesh.add_point(p, nn++);
882 6141 : };
883 :
884 1493752 : for (auto & elem : mesh.element_ptr_range())
885 : {
886 : // element_ptr_range skips deleted elements ... right?
887 21468 : libmesh_assert(elem);
888 21468 : libmesh_assert(elem->valid_id());
889 :
890 : // We only handle triangles in our triangulation
891 21468 : libmesh_assert_equal_to(elem->level(), 0u);
892 21468 : libmesh_assert_equal_to(elem->type(), TRI3);
893 :
894 : // If this triangle is as small as we desire, move along
895 762114 : if (!should_refine_elem(*elem))
896 684795 : continue;
897 :
898 : // Otherwise add a Steiner point. We'd like to add the
899 : // circumcenter ...
900 77319 : Point new_pt = elem->quasicircumcenter();
901 :
902 : // And to give it a node;
903 77319 : Node * new_node = nullptr;
904 :
905 : // But that might be outside our triangle, or even outside the
906 : // boundary. We'll find a triangle that should contain our new
907 : // point
908 77319 : Elem * cavity_elem = elem; // Start looking at elem anyway
909 :
910 : // We'll refine a boundary later if necessary.
911 20033 : auto boundary_refine = [this, &next_boundary_node,
912 : &cavity_elem, &new_node]
913 89700 : (unsigned int side)
914 : {
915 598 : libmesh_ignore(this); // Only used in dbg/devel
916 598 : libmesh_assert(new_node);
917 598 : libmesh_assert(new_node->valid_id());
918 :
919 21229 : Node * old_segment_start = cavity_elem->node_ptr(side),
920 21229 : * old_segment_end = cavity_elem->node_ptr((side+1)%3);
921 598 : libmesh_assert(old_segment_start);
922 598 : libmesh_assert(old_segment_start->valid_id());
923 598 : libmesh_assert(old_segment_end);
924 598 : libmesh_assert(old_segment_end->valid_id());
925 :
926 21229 : if (auto it = next_boundary_node.find(*old_segment_start);
927 598 : it != next_boundary_node.end())
928 : {
929 0 : libmesh_assert(it->second == old_segment_end);
930 0 : it->second = new_node;
931 : }
932 : else
933 : {
934 : // This would be an O(N) sanity check if we already
935 : // have a segments vector or any holes. :-P
936 598 : libmesh_assert(!this->segments.empty() ||
937 : (_holes && !_holes->empty()) ||
938 : (old_segment_end->id() ==
939 : old_segment_start->id() + 1));
940 21229 : next_boundary_node[*old_segment_start] = new_node;
941 : }
942 :
943 21229 : next_boundary_node[*new_node] = old_segment_end;
944 96370 : };
945 :
946 : // Let's find a triangle containing our new point, or at least
947 : // containing the end of a ray leading from our current triangle
948 : // to the new point.
949 77319 : Point ray_start = elem->vertex_average();
950 :
951 : // What side are we coming from, and what side are we going to?
952 2178 : unsigned int source_side = invalid_uint;
953 2178 : unsigned int side = invalid_uint;
954 :
955 85768 : while (!cavity_elem->contains_point(new_pt))
956 : {
957 16188 : side = segment_intersection(*cavity_elem, ray_start, new_pt, source_side);
958 :
959 456 : libmesh_assert_not_equal_to (side, invalid_uint);
960 :
961 16188 : Elem * neigh = cavity_elem->neighbor_ptr(side);
962 : // If we're on a boundary, stop there. Refine the boundary
963 : // if we're allowed, the boundary element otherwise.
964 16188 : if (!neigh)
965 : {
966 7739 : if (this->is_refine_boundary_allowed(boundary_info,
967 : *cavity_elem,
968 : side))
969 : {
970 7100 : new_pt = ray_start;
971 7100 : new_node = add_point(new_pt);
972 7100 : boundary_refine(side);
973 : }
974 : else
975 : {
976 : // Should we just add the vertex average of the
977 : // boundary element, to minimize the number of
978 : // slivers created?
979 : //
980 : // new_pt = cavity_elem->vertex_average();
981 : //
982 : // That works for a while, but it
983 : // seems to be able to "run away" and leave us with
984 : // crazy slivers on boundaries if we push interior
985 : // refinement too far while disabling boundary
986 : // refinement.
987 : //
988 : // Let's go back to refining our original problem
989 : // element.
990 639 : cavity_elem = elem;
991 639 : new_pt = cavity_elem->vertex_average();
992 639 : new_node = add_point(new_pt);
993 :
994 : // This was going to be a side refinement but it's
995 : // now an internal refinement
996 18 : side = invalid_uint;
997 : }
998 :
999 218 : break;
1000 : }
1001 :
1002 8449 : source_side = neigh->which_neighbor_am_i(cavity_elem);
1003 8449 : cavity_elem = neigh;
1004 238 : side = invalid_uint;
1005 : }
1006 :
1007 : // If we're ready to create a new node and we're not on a
1008 : // boundary ... should we be? We don't want to create any
1009 : // sliver elements or confuse poly2tri or anything.
1010 77319 : if (side == invalid_uint && !new_node)
1011 : {
1012 1960 : unsigned int worst_side = libMesh::invalid_uint;
1013 1960 : Real worst_cos = 1;
1014 278320 : for (auto s : make_range(3u))
1015 : {
1016 : // We never snap to a non-domain-boundary
1017 214620 : if (cavity_elem->neighbor_ptr(s))
1018 177951 : continue;
1019 :
1020 25631 : Real ax = cavity_elem->point(s)(0) - new_pt(0),
1021 25631 : ay = cavity_elem->point(s)(1) - new_pt(1),
1022 25631 : bx = cavity_elem->point((s+1)%3)(0) - new_pt(0),
1023 25631 : by = cavity_elem->point((s+1)%3)(1) - new_pt(1);
1024 27075 : const Real my_cos = (ax*bx+ay*by) /
1025 25631 : std::sqrt(ax*ax+ay*ay) /
1026 25631 : std::sqrt(bx*bx+by*by);
1027 :
1028 25631 : if (my_cos < worst_cos)
1029 : {
1030 700 : worst_side = s;
1031 700 : worst_cos = my_cos;
1032 : }
1033 : }
1034 :
1035 : // If we'd create a sliver element on the side, let's just
1036 : // refine the side instead, if we're allowed.
1037 69580 : if (worst_cos < -0.6) // -0.5 is the best we could enforce?
1038 : {
1039 412 : side = worst_side;
1040 :
1041 14626 : if (this->is_refine_boundary_allowed(boundary_info,
1042 : *cavity_elem,
1043 : side))
1044 : {
1045 : // Let's just try bisecting for now
1046 14129 : new_pt = (cavity_elem->point(side) +
1047 14129 : cavity_elem->point((side+1)%3)) / 2;
1048 14129 : new_node = add_point(new_pt);
1049 14129 : boundary_refine(side);
1050 : }
1051 : else // Do the best we can under these restrictions
1052 : {
1053 497 : new_pt = cavity_elem->vertex_average();
1054 497 : new_node = add_point(new_pt);
1055 :
1056 : // This was going to be a side refinement but it's
1057 : // now an internal refinement
1058 14 : side = invalid_uint;
1059 : }
1060 : }
1061 : else
1062 55366 : new_node = add_point(new_pt);
1063 : }
1064 : else
1065 218 : libmesh_assert(new_node);
1066 :
1067 : // Find the Delaunay cavity around the new point.
1068 4356 : std::set<Elem *, decltype(comp)> cavity(comp);
1069 :
1070 79497 : std::set<Elem *, decltype(comp)> unchecked_cavity ({cavity_elem}, comp);
1071 242607 : while (!unchecked_cavity.empty())
1072 : {
1073 9312 : std::set<Elem *, decltype(comp)> checking_cavity(comp);
1074 4656 : checking_cavity.swap(unchecked_cavity);
1075 385530 : for (Elem * checking_elem : checking_cavity)
1076 : {
1077 880968 : for (auto s : make_range(3u))
1078 : {
1079 660726 : Elem * neigh = checking_elem->neighbor_ptr(s);
1080 660726 : if (!neigh || checking_cavity.count(neigh) || cavity.count(neigh))
1081 204693 : continue;
1082 :
1083 456033 : if (in_circumcircle(*neigh, new_pt, TOLERANCE*TOLERANCE))
1084 138897 : unchecked_cavity.insert(neigh);
1085 : }
1086 : }
1087 :
1088 4656 : libmesh_merge_move(cavity, checking_cavity);
1089 : }
1090 :
1091 : // Retriangulate the Delaunay cavity.
1092 : // Each of our cavity triangle edges that are exterior to the
1093 : // cavity will be a source of one new triangle.
1094 :
1095 : // Set of elements that might need Delaunay swaps
1096 4356 : std::set<Elem *, decltype(comp)> check_delaunay_on(comp);
1097 :
1098 : // Keep maps for doing neighbor pointer assignment. Not going
1099 : // to iterate through these so hashing pointers is fine.
1100 : std::unordered_map<Node *, std::pair<Elem *, boundary_id_type>>
1101 4356 : neighbors_CCW, neighbors_CW;
1102 :
1103 297561 : for (Elem * old_elem : cavity)
1104 : {
1105 880968 : for (auto s : make_range(3u))
1106 : {
1107 660726 : Elem * neigh = old_elem->neighbor_ptr(s);
1108 660726 : if (!neigh || !cavity.count(neigh))
1109 : {
1110 374880 : Node * node_CW = old_elem->node_ptr(s),
1111 374880 : * node_CCW = old_elem->node_ptr((s+1)%3);
1112 :
1113 : auto set_neighbors =
1114 353760 : [&neighbors_CW, &neighbors_CCW, &node_CW,
1115 : &node_CCW, &boundary_info]
1116 1178636 : (Elem * new_neigh, boundary_id_type bcid)
1117 : {
1118 : // Set clockwise neighbor and vice-versa if possible
1119 374880 : if (const auto CW_it = neighbors_CW.find(node_CW);
1120 10560 : CW_it == neighbors_CW.end())
1121 : {
1122 5414 : libmesh_assert(!neighbors_CCW.count(node_CW));
1123 15974 : neighbors_CCW[node_CW] = std::make_pair(new_neigh, bcid);
1124 : }
1125 : else
1126 : {
1127 182683 : Elem * neigh_CW = CW_it->second.first;
1128 182683 : if (new_neigh)
1129 : {
1130 9688 : new_neigh->set_neighbor(0, neigh_CW);
1131 171962 : boundary_id_type bcid_CW = CW_it->second.second;
1132 171962 : if (bcid_CW != BoundaryInfo::invalid_id)
1133 15123 : boundary_info.add_side(new_neigh, 0, bcid_CW);
1134 :
1135 : }
1136 182683 : if (neigh_CW)
1137 : {
1138 9440 : neigh_CW->set_neighbor(2, new_neigh);
1139 167560 : if (bcid != BoundaryInfo::invalid_id)
1140 10721 : boundary_info.add_side(neigh_CW, 2, bcid);
1141 : }
1142 5146 : neighbors_CW.erase(CW_it);
1143 : }
1144 :
1145 : // Set counter-CW neighbor and vice-versa if possible
1146 374880 : if (const auto CCW_it = neighbors_CCW.find(node_CCW);
1147 10560 : CCW_it == neighbors_CCW.end())
1148 : {
1149 5146 : libmesh_assert(!neighbors_CW.count(node_CCW));
1150 5146 : neighbors_CW[node_CCW] = std::make_pair(new_neigh, bcid);
1151 : }
1152 : else
1153 : {
1154 192197 : Elem * neigh_CCW = CCW_it->second.first;
1155 192197 : if (new_neigh)
1156 : {
1157 186091 : boundary_id_type bcid_CCW = CCW_it->second.second;
1158 10484 : new_neigh->set_neighbor(2, neigh_CCW);
1159 186091 : if (bcid_CCW != BoundaryInfo::invalid_id)
1160 10508 : boundary_info.add_side(new_neigh, 2, bcid_CCW);
1161 : }
1162 192197 : if (neigh_CCW)
1163 : {
1164 10236 : neigh_CCW->set_neighbor(0, new_neigh);
1165 181689 : if (bcid != BoundaryInfo::invalid_id)
1166 6106 : boundary_info.add_side(neigh_CCW, 0, bcid);
1167 : }
1168 5414 : neighbors_CCW.erase(CCW_it);
1169 : }
1170 739200 : };
1171 :
1172 : // We aren't going to try to add a sliver element if we
1173 : // have a new boundary node here. We do need to
1174 : // keep track of other elements' neighbors, though.
1175 374880 : if (old_elem == cavity_elem &&
1176 3040 : s == side)
1177 : {
1178 1196 : std::vector<boundary_id_type> bcids;
1179 21229 : boundary_info.boundary_ids(old_elem, s, bcids);
1180 598 : libmesh_assert_equal_to(bcids.size(), 1);
1181 21229 : set_neighbors(nullptr, bcids[0]);
1182 598 : continue;
1183 : }
1184 :
1185 363613 : auto new_elem = Elem::build_with_id(TRI3, ne++);
1186 353651 : new_elem->set_node(0, new_node);
1187 353651 : new_elem->set_node(1, node_CW);
1188 353651 : new_elem->set_node(2, node_CCW);
1189 9962 : libmesh_assert(!new_elem->is_flipped());
1190 :
1191 : // Set in-and-out-of-cavity neighbor pointers
1192 19924 : new_elem->set_neighbor(1, neigh);
1193 353651 : if (neigh)
1194 : {
1195 : const unsigned int neigh_s =
1196 313110 : neigh->which_neighbor_am_i(old_elem);
1197 26460 : neigh->set_neighbor(neigh_s, new_elem.get());
1198 : }
1199 : else
1200 : {
1201 2284 : std::vector<boundary_id_type> bcids;
1202 40541 : boundary_info.boundary_ids(old_elem, s, bcids);
1203 40541 : boundary_info.add_side(new_elem.get(), 1, bcids);
1204 : }
1205 :
1206 : // Set in-cavity neighbors' neighbor pointers
1207 363613 : set_neighbors(new_elem.get(), BoundaryInfo::invalid_id);
1208 :
1209 : // C++ allows function argument evaluation in any
1210 : // order, but we need get() to precede move
1211 353651 : Elem * new_elem_ptr = new_elem.get();
1212 343689 : new_elems.emplace(new_elem_ptr, std::move(new_elem));
1213 :
1214 343689 : check_delaunay_on.insert(new_elem_ptr);
1215 333727 : }
1216 : }
1217 :
1218 220242 : boundary_info.remove(old_elem);
1219 : }
1220 :
1221 : // Now that we're done using our cavity elems (including with a
1222 : // cavity.find() that used a comparator that dereferences the
1223 : // elements!) it's safe to delete them.
1224 297561 : for (Elem * old_elem : cavity)
1225 : {
1226 220242 : if (const auto it = new_elems.find(old_elem);
1227 6204 : it == new_elems.end())
1228 154709 : mesh.delete_elem(old_elem);
1229 : else
1230 63687 : new_elems.erase(it);
1231 : }
1232 :
1233 : // Everybody found their match?
1234 2178 : libmesh_assert(neighbors_CW.empty());
1235 2178 : libmesh_assert(neighbors_CCW.empty());
1236 :
1237 : // Because we're preserving boundaries here, our naive cavity
1238 : // triangulation might not be a Delaunay triangulation. Let's
1239 : // check and if necessary fix that; we depend on it when doing
1240 : // future point insertions.
1241 77319 : restore_delaunay(check_delaunay_on, boundary_info);
1242 :
1243 : // This is too expensive to do on every cavity in devel mode
1244 : #ifdef DEBUG
1245 2178 : libmesh_assert_delaunay(mesh, new_elems);
1246 : #endif
1247 5963 : }
1248 :
1249 : // If we added any new boundary nodes, we're going to need to keep
1250 : // track of the changes they made to the outer polyline and/or to
1251 : // any holes.
1252 6319 : if (!next_boundary_node.empty())
1253 : {
1254 : auto checked_emplace = [this](dof_id_type new_first,
1255 11398 : dof_id_type new_second)
1256 : {
1257 : #ifdef DEBUG
1258 70128 : for (auto [first, second] : this->segments)
1259 : {
1260 67250 : libmesh_assert_not_equal_to(first, new_first);
1261 67250 : libmesh_assert_not_equal_to(second, new_second);
1262 : }
1263 2878 : if (!this->segments.empty())
1264 2764 : libmesh_assert_equal_to(this->segments.back().second, new_first);
1265 : #endif
1266 2878 : libmesh_assert_not_equal_to(new_first, new_second);
1267 :
1268 99291 : this->segments.emplace_back (new_first, new_second);
1269 99291 : };
1270 :
1271 : // Keep track of the outer polyline
1272 4047 : if (this->segments.empty())
1273 : {
1274 0 : dof_id_type last_id = DofObject::invalid_id;
1275 :
1276 : // Custom loop because we increment node_it 1+ times inside
1277 0 : for (auto node_it = _mesh.nodes_begin(),
1278 0 : node_end = _mesh.nodes_end();
1279 0 : node_it != node_end;)
1280 : {
1281 0 : Node & node = **node_it;
1282 0 : ++node_it;
1283 :
1284 0 : const dof_id_type node_id = node.id();
1285 :
1286 : // Don't add Steiner points
1287 0 : if (node_id >= _n_boundary_nodes)
1288 0 : break;
1289 :
1290 : // Connect up the previous node, if we didn't already
1291 : // connect it after some newly inserted nodes
1292 0 : if (!this->segments.empty())
1293 0 : last_id = this->segments.back().second;
1294 :
1295 0 : if (last_id != DofObject::invalid_id &&
1296 0 : last_id != node_id)
1297 0 : checked_emplace(last_id, node_id);
1298 :
1299 0 : last_id = node_id;
1300 :
1301 : // Connect to any newly inserted nodes
1302 0 : Node * this_node = &node;
1303 0 : auto it = next_boundary_node.find(*this_node);
1304 0 : while (it != next_boundary_node.end())
1305 : {
1306 0 : libmesh_assert(this_node->valid_id());
1307 0 : Node * next_node = it->second;
1308 0 : libmesh_assert(next_node->valid_id());
1309 :
1310 0 : if (node_it != node_end &&
1311 0 : next_node == *node_it)
1312 0 : ++node_it;
1313 :
1314 0 : checked_emplace(this_node->id(), next_node->id());
1315 :
1316 0 : this_node = next_node;
1317 0 : if (this_node->id() == this->segments.front().first)
1318 0 : break;
1319 :
1320 0 : it = next_boundary_node.find(*this_node);
1321 : }
1322 : }
1323 :
1324 : // We expect a closed loop here
1325 0 : if (this->segments.back().second != this->segments.front().first)
1326 0 : checked_emplace(this->segments.back().second,
1327 0 : this->segments.front().first);
1328 : }
1329 : else
1330 : {
1331 228 : std::vector<std::pair<unsigned int, unsigned int>> old_segments;
1332 4047 : old_segments.swap(this->segments);
1333 :
1334 114 : auto old_it = old_segments.begin();
1335 :
1336 4047 : const Node * node = _mesh.node_ptr(old_it->first);
1337 114 : const Node * const first_node = node;
1338 :
1339 2764 : do
1340 : {
1341 5756 : const dof_id_type node_id = node->id();
1342 102169 : if (const auto it = next_boundary_node.find(*node);
1343 2878 : it == next_boundary_node.end())
1344 : {
1345 139897 : while (node_id != old_it->first)
1346 : {
1347 2124 : ++old_it;
1348 2124 : libmesh_assert(old_it != old_segments.end());
1349 : }
1350 64539 : node = mesh.node_ptr(old_it->second);
1351 : }
1352 : else
1353 : {
1354 37630 : node = it->second;
1355 : }
1356 :
1357 102169 : checked_emplace(node_id, node->id());
1358 : }
1359 102169 : while (node != first_node);
1360 : }
1361 :
1362 : // Keep track of any holes
1363 4047 : if (this->_holes)
1364 : {
1365 : // Do we have any holes that need to be newly replaced?
1366 4686 : for (const Hole * hole : *this->_holes)
1367 : {
1368 1071 : if (this->replaced_holes.count(hole))
1369 1065 : continue;
1370 :
1371 36 : bool hole_point_insertion = false;
1372 5822 : for (auto p : make_range(hole->n_points()))
1373 9240 : if (next_boundary_node.count(hole->point(p)))
1374 : {
1375 4 : hole_point_insertion = true;
1376 4 : break;
1377 : }
1378 1278 : if (hole_point_insertion)
1379 : this->replaced_holes.emplace
1380 146 : (hole, std::make_unique<ArbitraryHole>(*hole));
1381 : }
1382 :
1383 : // If we have any holes that are being replaced, make sure
1384 : // their replacements are up to date.
1385 4686 : for (const Hole * hole : *this->_holes)
1386 : {
1387 66 : auto hole_it = replaced_holes.find(hole);
1388 2343 : if (hole_it == replaced_holes.end())
1389 1704 : continue;
1390 :
1391 34 : ArbitraryHole & arb = *hole_it->second;
1392 :
1393 : // We only need to update a replacement that's just had
1394 : // new points inserted
1395 34 : bool point_inserted = false;
1396 12496 : for (const Point & point : arb.get_points())
1397 672 : if (next_boundary_node.count(point))
1398 : {
1399 18 : point_inserted = true;
1400 18 : break;
1401 : }
1402 :
1403 1207 : if (!point_inserted)
1404 552 : continue;
1405 :
1406 : // Find all points in the replacement hole
1407 18 : std::vector<Point> new_points;
1408 :
1409 : // Our outer polyline is expected to have points in
1410 : // counter-clockwise order, so it proceeds "to the left"
1411 : // from the point of view of rays inside the domain
1412 : // pointing outward, and our next_boundary_node ordering
1413 : // was filled accordingly.
1414 : //
1415 : // Our inner holes are expected to have points in
1416 : // counter-clockwise order, but for holes "to the left
1417 : // as viewed from the hole interior is the *opposite* of
1418 : // "to the left as viewed from the domain interior". We
1419 : // need to build the updated hole ordering "backwards".
1420 :
1421 : // We should never see duplicate points when we add one
1422 : // to a hole; if we do then we did something wrong.
1423 816 : auto push_back_new_point = [&new_points](const Point & p) {
1424 : // O(1) assert in devel
1425 278 : libmesh_assert(new_points.empty() ||
1426 : new_points.back() != p);
1427 : #ifdef DEBUG
1428 : // O(N) asserts in dbg
1429 2620 : for (auto old_p : new_points)
1430 2342 : libmesh_assert_not_equal_to(old_p, p);
1431 : #endif
1432 9869 : new_points.push_back(p);
1433 9591 : };
1434 :
1435 326 : for (auto point_it = arb.get_points().rbegin(),
1436 18 : point_end = arb.get_points().rend();
1437 6106 : point_it != point_end;)
1438 : {
1439 5467 : Point point = *point_it;
1440 154 : ++point_it;
1441 :
1442 5603 : if (new_points.empty() ||
1443 136 : (point != new_points.back() &&
1444 136 : point != new_points.front()))
1445 154 : push_back_new_point(point);
1446 :
1447 154 : auto it = next_boundary_node.find(point);
1448 9869 : while (it != next_boundary_node.end())
1449 : {
1450 4828 : point = *it->second;
1451 136 : if (point == new_points.front())
1452 12 : break;
1453 4514 : if (point_it != point_end &&
1454 112 : point == *point_it)
1455 56 : ++point_it;
1456 124 : push_back_new_point(point);
1457 124 : it = next_boundary_node.find(point);
1458 : }
1459 : }
1460 :
1461 621 : std::reverse(new_points.begin(), new_points.end());
1462 :
1463 1242 : arb.set_points(std::move(new_points));
1464 : }
1465 : }
1466 : }
1467 :
1468 : // Okay, *now* we can add the new elements.
1469 294437 : for (auto & [raw_elem, unique_elem] : new_elems)
1470 : {
1471 8116 : libmesh_assert_equal_to(raw_elem, unique_elem.get());
1472 8116 : libmesh_assert(!raw_elem->is_flipped());
1473 8116 : libmesh_ignore(raw_elem); // Old gcc warns "unused variable"
1474 304350 : mesh.add_elem(std::move(unique_elem));
1475 : }
1476 :
1477 : // Did we add anything?
1478 6319 : return !new_elems.empty();
1479 : }
1480 :
1481 :
1482 762114 : bool Poly2TriTriangulator::should_refine_elem(Elem & elem)
1483 : {
1484 762114 : const Real min_area_target = this->desired_area();
1485 762114 : FunctionBase<Real> *area_func = this->has_auto_area_function() ? this->get_auto_area_function() : this->get_desired_area_function();
1486 :
1487 : // If this isn't a question, why are we here?
1488 21468 : libmesh_assert(min_area_target > 0 ||
1489 : area_func != nullptr ||
1490 : this->has_auto_area_function());
1491 :
1492 762114 : const Real area = elem.volume();
1493 :
1494 : // If we don't have position-dependent area targets we can make a
1495 : // decision quickly
1496 762114 : if (!area_func && !this->has_auto_area_function())
1497 181192 : return (area > min_area_target);
1498 16364 : else if(area_func && this->has_auto_area_function())
1499 : libmesh_warning("WARNING: both desired are function and automatic area function are set. Using automatic area function.");
1500 :
1501 : // If we do?
1502 : //
1503 : // See if we're meeting the local area target at all the elem
1504 : // vertices first
1505 2180552 : for (auto v : make_range(elem.n_vertices()))
1506 : {
1507 : // If we have an auto area function, we'll use it and override other area options
1508 1696520 : const Real local_area_target = (*area_func)(elem.point(v));
1509 1650040 : libmesh_error_msg_if
1510 : (local_area_target <= 0,
1511 : "Non-positive desired element areas are unachievable");
1512 1650040 : if (area > local_area_target)
1513 1420 : return true;
1514 : }
1515 :
1516 : // If our vertices are happy, it's still possible that our interior
1517 : // isn't. Are we allowed not to bother checking it?
1518 530512 : if (!min_area_target)
1519 14944 : return false;
1520 :
1521 0 : libmesh_not_implemented_msg
1522 : ("Combining a minimum desired_area with an area function isn't yet supported.");
1523 : }
1524 :
1525 :
1526 : } // namespace libMesh
1527 :
1528 :
1529 : // Unnecessary at end of file, *except* maybe we'll do a unity build
1530 : // someday, and in the meantime test coverage is nice.
1531 : #include "libmesh/restore_ieee754.h"
1532 :
1533 :
1534 : #endif // LIBMESH_HAVE_POLY2TRI
|
