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