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