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