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 : // Local includes
22 : #include "libmesh/mesh_triangle_holes.h"
23 :
24 : #include "libmesh/boundary_info.h"
25 : #include "libmesh/elem.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/mesh_base.h"
28 : #include "libmesh/mesh_serializer.h"
29 : #include "libmesh/node.h"
30 : #include "libmesh/parallel_algebra.h" // Packing<Point>
31 : #include "libmesh/simple_range.h"
32 :
33 : // TIMPI includes
34 : #include "timpi/parallel_implementation.h" // broadcast
35 :
36 : // C++ includes
37 : #include <algorithm>
38 :
39 : namespace
40 : {
41 : using namespace libMesh;
42 :
43 24804 : int signof(Real val) {
44 840811 : return (0 < val) - (val < 0);
45 : }
46 :
47 : // Return 1 iff counter-clockwise turn
48 : // Return -1 iff clockwise turn
49 : // Return 0 iff collinear
50 617344 : int orientation(const Point & p0,
51 : const Point & p1,
52 : const Point & p2)
53 : {
54 636110 : const Real detleft = (p0(0)-p2(0))*(p1(1)-p2(1));
55 636110 : const Real detright = (p0(1)-p2(1))*(p1(0)-p2(0));
56 :
57 636110 : return signof(detleft - detright);
58 : }
59 :
60 : // Same, but for the ray target as it goes to infinity
61 198663 : int ray_orientation(const Point & p0,
62 : const Point & p1,
63 : const Point & source,
64 : const Point & ray_target)
65 : {
66 6038 : const Point rayvec = ray_target - source;
67 6038 : const Point edgevec = p1 - p0;
68 204701 : const Real det = edgevec(0)*rayvec(1)-edgevec(1)*rayvec(0);
69 :
70 204701 : return signof(det);
71 : }
72 :
73 : // Return true iff the ray from source toward ray_target crosses the
74 : // edge from edge_pt0 (non-inclusive) to edge_pt1 (inclusive) in the
75 : // forward ray direction. edge_pt2 is the vertex following edge_pt1
76 : // along the polygon; it is used to classify a crossing that lands
77 : // exactly on the shared vertex edge_pt1.
78 : //
79 : // The decision is made purely from orientation predicates on the
80 : // vertices, never from a parametric edge coordinate compared against
81 : // a tolerance. This is what guarantees consistency: a shared vertex
82 : // yields the exact same orientation when viewed from either of its
83 : // two edges, so a ray passing through a vertex is counted exactly
84 : // once for a "through" crossing and not at all for a "tangent"
85 : // (glancing) touch, regardless of floating-point round-off.
86 204701 : bool is_intersection(const Point & source,
87 : const Point & ray_target,
88 : const Point & edge_pt0,
89 : const Point & edge_pt1,
90 : const Point & edge_pt2)
91 : {
92 198663 : const int orient_st0 = orientation(source, ray_target, edge_pt0);
93 198663 : const int orient_st1 = orientation(source, ray_target, edge_pt1);
94 :
95 : // Which side of the edge line the source lies on, and which way the
96 : // ray turns off the edge, together tell us whether the crossing is
97 : // ahead of the source rather than behind it.
98 198663 : const int orient_edge_s = orientation(edge_pt0, edge_pt1, source);
99 198663 : const int orient_edge_t = ray_orientation(edge_pt0, edge_pt1, source, ray_target);
100 204701 : const bool forward = (orient_edge_s != orient_edge_t);
101 :
102 : // The ray line passes exactly through edge_pt1, the inclusive
103 : // endpoint of this edge. Count it only as a genuine "through"
104 : // crossing: the previous vertex (edge_pt0) and the next vertex
105 : // (edge_pt2) must lie on opposite sides of the ray line. A tangent
106 : // touch (both neighbors on the same side) or a collinear edge
107 : // (edge_pt0 also on the line) does not count.
108 204701 : if (orient_st1 == 0)
109 : {
110 21355 : const int orient_st2 = orientation(source, ray_target, edge_pt2);
111 39975 : return forward && (orient_st0 != 0) && (orient_st0 == -orient_st2);
112 : }
113 :
114 : // The ray line passes through edge_pt0, the non-inclusive endpoint;
115 : // that vertex is counted (if at all) by the previous edge.
116 : //
117 : // A run of two or more vertices all lying exactly on the ray line
118 : // is therefore never counted, even when the boundary arrives at the
119 : // run from one side and leaves from the other, which ought to count
120 : // once. Deciding that needs a look further along the polygon than
121 : // edge_pt2. This is long-standing behavior, not something the
122 : // orientation predicates changed, and it takes exact collinearity
123 : // to reach.
124 182694 : if (orient_st0 == 0)
125 524 : return false;
126 :
127 : // Ordinary case: count it iff the edge endpoints straddle the ray
128 : // line.
129 164963 : return forward && (orient_st0 == -orient_st1);
130 : }
131 :
132 : // Returns a non-negative distance iff the ray from source in the
133 : // direction of ray_target intersects the edge from pt0
134 : // (non-inclusive) to pt1 (inclusive), -1 otherwise.
135 : //
136 : // If the intersection is a "glancing" one at a corner, return -1.
137 : //
138 : // edge_pt2 is the vertex following edge_pt1, needed to tell a
139 : // "glancing" corner from one the ray passes through.
140 204701 : Real find_intersection(const Point & source,
141 : const Point & ray_target,
142 : const Point & edge_pt0,
143 : const Point & edge_pt1,
144 : const Point & edge_pt2)
145 : {
146 : // Decide whether the ray crosses this edge using only orientation
147 : // predicates on the edge vertices, so that shared vertices are
148 : // classified consistently between neighboring edges.
149 204701 : if (!is_intersection(source, ray_target, edge_pt0, edge_pt1, edge_pt2))
150 5352 : return -1;
151 :
152 : // Now that the crossing decision is settled, use parametric math
153 : // only to recover the distance to the intersection.
154 22544 : const Real raydx = ray_target(0)-source(0),
155 22544 : raydy = ray_target(1)-source(1),
156 22544 : edgedx = edge_pt1(0)-edge_pt0(0),
157 22544 : edgedy = edge_pt1(1)-edge_pt0(1);
158 22544 : const Real denom = edgedx * raydy - edgedy * raydx;
159 :
160 : // divide-by-zero means the segments are parallel; is_intersection
161 : // rejects that case, but guard against it anyway.
162 22544 : if (denom == 0)
163 0 : return -1;
164 :
165 22544 : const Real one_over_denom = 1 / denom;
166 :
167 22544 : const Real targetsdx = edge_pt1(0)-ray_target(0),
168 22544 : targetsdy = edge_pt1(1)-ray_target(1);
169 :
170 22544 : const Real u_num = targetsdx * edgedy - targetsdy * edgedx;
171 22544 : const Real u = u_num * one_over_denom;
172 22544 : const Real ray_fraction = (1-u);
173 :
174 : // is_intersection guarantees a forward crossing, so clamp away any
175 : // floating-point noise that would otherwise make a vertex hit read
176 : // as a tiny negative distance.
177 : const Real distance =
178 22544 : ray_fraction * std::sqrt(raydx*raydx + raydy*raydy);
179 44055 : return std::max(Real(0), distance);
180 : }
181 : }
182 :
183 :
184 : namespace libMesh
185 : {
186 :
187 : //
188 : // Hole member functions
189 : //
190 1144 : Real TriangulatorInterface::Hole::area() const
191 : {
192 1144 : return this->areavec().norm() / 2;
193 : }
194 :
195 :
196 1144 : RealGradient TriangulatorInterface::Hole::areavec() const
197 : {
198 1144 : const unsigned int np = this->n_points();
199 :
200 1144 : if (np < 3)
201 2 : return 0;
202 :
203 1073 : const Point p0 = this->point(0);
204 :
205 : // Every segment (p_{i-1},p_i) from i=2 on defines a triangle w.r.t.
206 : // p_0. Add up the cross products of those triangles. We'll save
207 : // the division by 2 and the norm for the end.
208 : //
209 : // Your hole points had best be coplanar, but this should work
210 : // regardless of which plane they're in. If you're in the XY plane,
211 : // then the standard counter-clockwise hole point ordering gives you
212 : // a positive areavec(2);
213 :
214 34 : RealGradient areavec = 0;
215 :
216 3661 : for (unsigned int i=2; i != np; ++i)
217 : {
218 2588 : const Point e_0im = this->point(i-1) - p0,
219 2588 : e_0i = this->point(i) - p0;
220 :
221 2500 : areavec += e_0i.cross(e_0im);
222 : }
223 :
224 1073 : return areavec;
225 : }
226 :
227 :
228 :
229 : std::vector<Real>
230 36630 : TriangulatorInterface::Hole::find_ray_intersections(Point ray_start,
231 : Point ray_target) const
232 : {
233 36630 : const auto np = this->n_points();
234 :
235 1096 : std::vector<Real> intersection_distances;
236 :
237 241331 : for (auto i : make_range(np))
238 : {
239 204701 : const Point & p0 = this->point(i),
240 204701 : & p1 = this->point((i+1)%np),
241 204701 : & p2 = this->point((i+2)%np);
242 : const Real intersection_distance =
243 204701 : find_intersection(ray_start, ray_target, p0, p1, p2);
244 204701 : if (intersection_distance >= 0)
245 : intersection_distances.push_back
246 22544 : (intersection_distance);
247 : }
248 :
249 36630 : return intersection_distances;
250 : }
251 :
252 :
253 :
254 1582 : Point TriangulatorInterface::Hole::calculate_inside_point() const
255 : {
256 : // Start with the vertex average
257 :
258 : // Turns out "I'm a fully compliant C++17 compiler!" doesn't
259 : // mean "I have a full C++17 standard library!"
260 : // inside = std::reduce(points.begin(), points.end());
261 54 : Point inside = 0;
262 16675 : for (auto i : make_range(this->n_points()))
263 15093 : inside += this->point(i);
264 :
265 1582 : inside /= this->n_points();
266 :
267 : // Count the number of intersections with a ray to the right,
268 : // keep track of how far they are
269 54 : Point ray_target = inside + Point(1);
270 : std::vector<Real> intersection_distances =
271 1636 : this->find_ray_intersections(inside, ray_target);
272 :
273 : // The vertex average isn't on the interior, and we found no
274 : // intersections to the right? Try looking to the left.
275 1582 : if (!intersection_distances.size())
276 : {
277 0 : ray_target = inside - Point(1);
278 : intersection_distances =
279 0 : this->find_ray_intersections(inside, ray_target);
280 : }
281 :
282 : // I'd make this an assert, but I'm not 100% confident we can't
283 : // get here via some kind of FP error on a weird hole shape.
284 1636 : libmesh_error_msg_if
285 : (!intersection_distances.size(),
286 : "Can't find a center for a MeshedHole!");
287 :
288 1582 : if (intersection_distances.size() % 2)
289 54 : return inside;
290 :
291 : // The vertex average is outside. So go from the vertex average to
292 : // the closest edge intersection, then halfway to the next-closest.
293 :
294 : // Find the nearest first.
295 0 : Real min_distance = std::numeric_limits<Real>::max(),
296 0 : second_distance = std::numeric_limits<Real>::max();
297 0 : for (Real d : intersection_distances)
298 0 : if (d < min_distance)
299 : {
300 0 : second_distance = min_distance;
301 0 : min_distance = d;
302 : }
303 :
304 0 : const Point ray = ray_target - inside;
305 0 : inside += ray * (min_distance + second_distance)/2;
306 :
307 0 : return inside;
308 : }
309 :
310 :
311 35048 : bool TriangulatorInterface::Hole::contains(Point p) const
312 : {
313 : // Count the number of intersections with a ray to the right,
314 : // keep track of how far they are
315 1042 : Point ray_target = p + Point(1);
316 : std::vector<Real> intersection_distances =
317 35048 : this->find_ray_intersections(p, ray_target);
318 :
319 : // Odd number of intersections == we're inside
320 : // Even number == we're outside
321 37132 : return intersection_distances.size() % 2;
322 : }
323 :
324 :
325 :
326 : //
327 : // PolygonHole member functions
328 : //
329 1365 : TriangulatorInterface::PolygonHole::PolygonHole(const Point & center,
330 : Real radius,
331 1365 : unsigned int n_points_in) :
332 1273 : _center(center),
333 1273 : _radius(radius),
334 1365 : _n_points(n_points_in)
335 1365 : {}
336 :
337 :
338 37148 : unsigned int TriangulatorInterface::PolygonHole::n_points() const
339 : {
340 37148 : return _n_points;
341 : }
342 :
343 :
344 334887 : Point TriangulatorInterface::PolygonHole::point(const unsigned int n) const
345 : {
346 : // The nth point lies at the angle theta = 2 * pi * n / _n_points
347 334887 : const Real theta = static_cast<Real>(n) * 2.0 * libMesh::pi / static_cast<Real>(_n_points);
348 :
349 344445 : return Point(_center(0) + _radius*std::cos(theta), // x=r*cos(theta)
350 334887 : _center(1) + _radius*std::sin(theta), // y=r*sin(theta)
351 354003 : 0.);
352 : }
353 :
354 :
355 154 : Point TriangulatorInterface::PolygonHole::inside() const
356 : {
357 : // The center of the hole is definitely inside.
358 154 : return _center;
359 : }
360 :
361 :
362 : //
363 : // AffineHole member functions
364 : //
365 284355 : Point TriangulatorInterface::AffineHole::point(const unsigned int n) const
366 : {
367 284355 : return this->transform(_underlying.point(n));
368 : }
369 :
370 :
371 0 : Point TriangulatorInterface::AffineHole::inside() const
372 : {
373 0 : return this->transform(_underlying.inside());
374 : }
375 :
376 :
377 284355 : Point TriangulatorInterface::AffineHole::transform(const Point & p) const
378 : {
379 284355 : const Real cos_a = std::cos(_angle);
380 284355 : const Real sin_a = std::sin(_angle);
381 284355 : return Point(p(0)*cos_a-p(1)*sin_a + _shift(0),
382 308385 : p(1)*cos_a+p(1)*sin_a + _shift(1));
383 : }
384 :
385 :
386 : //
387 : // ArbitraryHole member functions
388 : //
389 288 : TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Point & center,
390 288 : std::vector<Point> points)
391 268 : : _center(center),
392 298 : _points(std::move(points))
393 : {
394 288 : _segment_indices.push_back(0);
395 298 : _segment_indices.push_back(_points.size());
396 288 : }
397 :
398 :
399 0 : TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Point & center,
400 : std::vector<Point> points,
401 0 : std::vector<unsigned int> segment_indices)
402 0 : : _center(center),
403 0 : _points(std::move(points)),
404 0 : _segment_indices(std::move(segment_indices))
405 0 : {}
406 :
407 :
408 793 : TriangulatorInterface::ArbitraryHole::ArbitraryHole(std::vector<Point> points)
409 849 : : _points(std::move(points))
410 : {
411 793 : _segment_indices.push_back(0);
412 821 : _segment_indices.push_back(_points.size());
413 793 : _center = this->calculate_inside_point();
414 793 : }
415 :
416 :
417 142 : TriangulatorInterface::ArbitraryHole::ArbitraryHole(const Hole & orig)
418 142 : : _center(orig.inside())
419 : {
420 142 : const unsigned int np = orig.n_points();
421 142 : _points.reserve(np);
422 710 : for (auto i : make_range(np))
423 1120 : _points.push_back(orig.point(i));
424 142 : }
425 :
426 :
427 16101 : unsigned int TriangulatorInterface::ArbitraryHole::n_points() const
428 : {
429 16613 : return _points.size();
430 : }
431 :
432 :
433 350162 : Point TriangulatorInterface::ArbitraryHole::point(const unsigned int n) const
434 : {
435 10596 : libmesh_assert_less (n, _points.size());
436 360758 : return _points[n];
437 : }
438 :
439 :
440 4 : Point TriangulatorInterface::ArbitraryHole::inside() const
441 : {
442 4 : return _center;
443 : }
444 :
445 :
446 16 : std::vector<unsigned int> TriangulatorInterface::ArbitraryHole::segment_indices() const
447 : {
448 16 : return _segment_indices;
449 : }
450 :
451 :
452 : //
453 : // MeshedHole member functions
454 : //
455 :
456 :
457 943 : TriangulatorInterface::MeshedHole::MeshedHole(const MeshBase & mesh,
458 943 : std::set<std::size_t> ids)
459 1033 : : _center(std::numeric_limits<Real>::max())
460 : {
461 : // We'll want to do this on one processor and broadcast to the rest;
462 : // otherwise we can get out of sync by doing things like using
463 : // pointers as keys.
464 36 : libmesh_parallel_only(mesh.comm());
465 :
466 : MeshSerializer serial(const_cast<MeshBase &>(mesh),
467 970 : /* serial */ true, /* only proc 0 */ true);
468 :
469 : // Try to keep in sync even if we throw an error on proc 0, so we
470 : // can examine errors in our unit tests in parallel too.
471 42 : std::string error_reported;
472 :
473 78 : auto report_error = [&mesh, &error_reported](std::string er) {
474 36 : error_reported = std::move(er);
475 36 : mesh.comm().broadcast(error_reported);
476 102 : libmesh_error_msg(error_reported);
477 943 : };
478 :
479 979 : if (mesh.processor_id() != 0)
480 : {
481 : // Make sure proc 0 didn't just fail
482 777 : mesh.comm().broadcast(error_reported);
483 1125 : libmesh_error_msg_if(!error_reported.empty(), error_reported);
484 :
485 : // Receive the points proc 0 will send later
486 600 : mesh.comm().broadcast(_points);
487 600 : mesh.comm().broadcast(_midpoints);
488 15 : return;
489 : }
490 :
491 : // We'll find all the line segments first, then stitch them together
492 : // afterward. If the line segments come from 2D element sides then
493 : // we'll label their edge_type as "1" for clockwise orientation
494 : // around the element or "2" for CCW, to make it easier to detect
495 : // and scream about cases where we have a disconnected outer
496 : // boundary.
497 : std::multimap<const Node *,
498 36 : std::pair<const Node *, int>> hole_edge_map;
499 :
500 : // If we're looking at higher-order elements, we have mid-edge edge
501 : // nodes to worry about. hole_midpoint_map[{m,n}][i] should give us
502 : // the ith mid-edge node traveling from vertex m to vertex n
503 : std::map<std::pair<const Node *, const Node *>,
504 36 : std::vector<const Node *>> hole_midpoint_map;
505 :
506 36 : std::vector<boundary_id_type> bcids;
507 :
508 18 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
509 :
510 2068 : for (const auto & elem : mesh.active_element_ptr_range())
511 : {
512 994 : if (elem->dim() == 1)
513 : {
514 634 : if (ids.empty() || ids.count(elem->subdomain_id()))
515 : {
516 504 : hole_edge_map.emplace(elem->node_ptr(0),
517 452 : std::make_pair(elem->node_ptr(1),
518 156 : /*edge*/ 0));
519 504 : hole_edge_map.emplace(elem->node_ptr(1),
520 452 : std::make_pair(elem->node_ptr(0),
521 156 : /*edge*/ 0));
522 504 : if (elem->type() == EDGE3)
523 : {
524 213 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(0),
525 259 : elem->node_ptr(1)),
526 282 : std::vector<const Node *>{elem->node_ptr(2)});
527 213 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(1),
528 259 : elem->node_ptr(0)),
529 495 : std::vector<const Node *>{elem->node_ptr(2)});
530 : }
531 268 : else if (elem->type() == EDGE4)
532 : {
533 0 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(0),
534 0 : elem->node_ptr(1)),
535 0 : std::vector<const Node *>{elem->node_ptr(2),
536 0 : elem->node_ptr(3)});
537 0 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(1),
538 0 : elem->node_ptr(0)),
539 0 : std::vector<const Node *>{elem->node_ptr(3),
540 0 : elem->node_ptr(2)});
541 : }
542 : else
543 29 : libmesh_assert_equal_to(elem->default_side_order(), 1);
544 : }
545 634 : continue;
546 : }
547 :
548 360 : if (elem->dim() == 2)
549 : {
550 360 : const auto ns = elem->n_sides();
551 1776 : for (auto s : make_range(ns))
552 : {
553 1416 : boundary_info.boundary_ids(elem, s, bcids);
554 :
555 198 : bool add_edge = false;
556 1614 : if (!elem->neighbor_ptr(s) && ids.empty())
557 38 : add_edge = true;
558 :
559 198 : if (!add_edge)
560 1120 : for (auto b : bcids)
561 0 : if (ids.count(b))
562 0 : add_edge = true;
563 :
564 1196 : if (add_edge)
565 : {
566 296 : hole_edge_map.emplace(elem->node_ptr(s),
567 258 : std::make_pair(elem->node_ptr((s+1)%ns),
568 114 : /*counter-CW*/ 2));
569 : // Do we really need to support flipped 2D elements?
570 296 : hole_edge_map.emplace(elem->node_ptr((s+1)%ns),
571 258 : std::make_pair(elem->node_ptr(s),
572 114 : /*clockwise*/ 1));
573 :
574 296 : if (elem->default_side_order() == 2)
575 : {
576 96 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr(s),
577 128 : elem->node_ptr((s+1)%ns)),
578 144 : std::vector<const Node *>{elem->node_ptr(s+ns)});
579 96 : hole_midpoint_map.emplace(std::make_pair(elem->node_ptr((s+1)%ns),
580 128 : elem->node_ptr(s)),
581 240 : std::vector<const Node *>{elem->node_ptr(s+ns)});
582 : }
583 : else
584 22 : libmesh_assert_equal_to(elem->default_side_order(), 1);
585 :
586 296 : continue;
587 : }
588 : }
589 : }
590 130 : }
591 :
592 166 : if (hole_edge_map.empty())
593 30 : report_error("No valid hole edges found in mesh!");
594 :
595 : // Function to pull a vector of points out of the map; a loop of
596 : // edges connecting these points defines a hole boundary. If the
597 : // mesh has multiple boundaries (e.g. because it had holes itself),
598 : // then a random vector will be extracted; this function will be
599 : // called multiple times so that the various options can be
600 : // compared. We choose the largest option.
601 : auto extract_edge_vector =
602 1848 : [&report_error, &hole_edge_map, &hole_midpoint_map]() {
603 : std::tuple<std::vector<const Node *>, std::vector<const Node *>, int>
604 : hole_points_and_edge_type
605 230 : {{hole_edge_map.begin()->first, hole_edge_map.begin()->second.first},
606 400 : {}, hole_edge_map.begin()->second.second};
607 :
608 20 : auto & hole_points = std::get<0>(hole_points_and_edge_type);
609 20 : auto & midpoint_points = std::get<1>(hole_points_and_edge_type);
610 20 : int & edge_type = std::get<2>(hole_points_and_edge_type);
611 :
612 : // We won't be needing to search for this edge
613 20 : hole_edge_map.erase(hole_points.front());
614 :
615 : // Sort the remaining edges into a connected order
616 190 : for (const Node * last = hole_points.front(),
617 190 : * n = hole_points.back();
618 790 : n != hole_points.front();
619 138 : last = n,
620 600 : n = hole_points.back())
621 : {
622 69 : auto [next_it_begin, next_it_end] = hole_edge_map.equal_range(n);
623 :
624 612 : if (std::distance(next_it_begin, next_it_end) != 2)
625 15 : report_error("Bad edge topology found by MeshedHole");
626 :
627 600 : const Node * next = nullptr;
628 1800 : for (const auto & [key, val] : as_range(next_it_begin, next_it_end))
629 : {
630 136 : libmesh_assert_equal_to(key, n);
631 136 : libmesh_ignore(key);
632 136 : libmesh_assert_not_equal_to(val.first, n);
633 :
634 : // Don't go backwards on the edge we just traversed
635 1200 : if (val.first == last)
636 530 : continue;
637 :
638 : // We can support mixes of Edge and Tri-side edges, but we
639 : // can't do proper error detection on flipped triangles.
640 600 : if (val.second != edge_type &&
641 0 : val.second != 0)
642 : {
643 0 : if (!edge_type)
644 0 : edge_type = val.second;
645 : else
646 0 : report_error("MeshedHole sees inconsistent triangle orientations on boundary");
647 : }
648 600 : next = val.first;
649 : }
650 :
651 : // We should never hit the same n twice!
652 68 : hole_edge_map.erase(next_it_begin, next_it_end);
653 :
654 600 : hole_points.push_back(next);
655 : }
656 :
657 961 : for (auto i : make_range(hole_points.size()-1))
658 : {
659 938 : const auto & midpoints = hole_midpoint_map[{hole_points[i],hole_points[i+1]}];
660 687 : midpoint_points.insert(midpoint_points.end(),
661 349 : midpoints.begin(), midpoints.end());
662 : }
663 :
664 19 : hole_points.pop_back();
665 :
666 178 : return hole_points_and_edge_type;
667 166 : };
668 :
669 : /*
670 : * If it's not obvious which loop we find is really the loop we
671 : * want, then we should die with a nice error message.
672 : */
673 18 : int n_negative_areas = 0,
674 18 : n_positive_areas = 0,
675 18 : n_edgeelem_loops = 0;
676 :
677 39 : std::vector<const Node *> outer_hole_points, outer_mid_points;
678 18 : int outer_edge_type = -1;
679 18 : Real twice_outer_area = 0,
680 18 : abs_twice_outer_area = 0;
681 :
682 : #ifdef DEBUG
683 : // Area and edge type, for error reporting
684 36 : std::vector<std::pair<Real, int>> areas;
685 : #endif
686 :
687 332 : while (!hole_edge_map.empty()) {
688 209 : auto [hole_points, mid_points, edge_type] = extract_edge_vector();
689 :
690 178 : if (edge_type == 0)
691 : {
692 126 : ++n_edgeelem_loops;
693 126 : if (n_edgeelem_loops > 1)
694 15 : report_error("MeshedHole is confused by multiple loops of Edge elements");
695 114 : if (n_positive_areas || n_negative_areas)
696 0 : report_error("MeshedHole is confused by meshes with both Edge and 2D-side boundaries");
697 : }
698 :
699 36 : const std::size_t n_hole_points = hole_points.size();
700 166 : if (n_hole_points < 3)
701 11 : report_error("Loop with only " + std::to_string(n_hole_points) +
702 : " hole edges found in mesh!");
703 :
704 18 : Real twice_this_area = 0;
705 166 : const Point p0 = *hole_points[0];
706 562 : for (unsigned int i=2; i != n_hole_points; ++i)
707 : {
708 444 : const Point e_0im = *hole_points[i-1] - p0,
709 396 : e_0i = *hole_points[i] - p0;
710 :
711 396 : twice_this_area += e_0i.cross(e_0im)(2);
712 : }
713 :
714 18 : auto abs_twice_this_area = std::abs(twice_this_area);
715 :
716 166 : if (((twice_this_area > 0) && edge_type == 2) ||
717 77 : ((twice_this_area < 0) && edge_type == 1))
718 0 : ++n_positive_areas;
719 166 : else if (edge_type != 0)
720 52 : ++n_negative_areas;
721 :
722 : #ifdef DEBUG
723 18 : areas.push_back({twice_this_area/2,edge_type});
724 : #endif
725 :
726 166 : if (abs_twice_this_area > abs_twice_outer_area)
727 : {
728 17 : twice_outer_area = twice_this_area;
729 17 : abs_twice_outer_area = abs_twice_this_area;
730 17 : outer_hole_points = std::move(hole_points);
731 17 : outer_mid_points = std::move(mid_points);
732 154 : outer_edge_type = edge_type;
733 : }
734 : }
735 :
736 158 : _points.resize(outer_hole_points.size());
737 : std::transform(outer_hole_points.begin(),
738 : outer_hole_points.end(),
739 : _points.begin(),
740 672 : [](const Node * n){ return Point(*n); });
741 158 : _midpoints.resize(outer_mid_points.size());
742 : std::transform(outer_mid_points.begin(),
743 : outer_mid_points.end(),
744 : _midpoints.begin(),
745 364 : [](const Node * n){ return Point(*n); });
746 :
747 142 : if (!twice_outer_area)
748 0 : report_error("Zero-area MeshedHoles are not currently supported");
749 :
750 : // We ordered ourselves counter-clockwise? But a hole is expected
751 : // to be clockwise, so use the reverse order.
752 142 : if (twice_outer_area > 0)
753 : {
754 65 : std::reverse(_points.begin(), _points.end());
755 :
756 : // Our midpoints are numbered e.g.
757 : // (01a)(01b)(12a)(12b)(23a)(23b)(30a)(30b) for points 0123, but
758 : // if we reverse to get 3210 then we want our midpoints to be
759 : // (23b)(23a)(12b)(12a)(01b)(01a)(30b)(30a)
760 91 : const unsigned int n_midpoints = _midpoints.size() / _points.size();
761 6 : auto split_it = _midpoints.end() - n_midpoints;
762 65 : std::reverse(_midpoints.begin(), split_it);
763 65 : std::reverse(split_it, _midpoints.end());
764 : }
765 :
766 : #ifdef DEBUG
767 1 : auto print_areas = [areas](){
768 1 : libMesh::out << "Found boundary areas:\n";
769 4 : static const std::vector<std::string> edgenames {"E","CW","CCW"};
770 3 : for (auto area : areas)
771 2 : libMesh::out << '(' << edgenames[area.second] << ' ' <<
772 2 : area.first << ')';
773 1 : libMesh::out << std::endl;
774 33 : };
775 : #else
776 : auto print_areas = [](){};
777 : #endif
778 :
779 142 : if (((twice_outer_area > 0) && outer_edge_type == 2) ||
780 64 : ((twice_outer_area < 0) && outer_edge_type == 1))
781 : {
782 0 : if (n_positive_areas > 1)
783 : {
784 0 : print_areas();
785 0 : report_error("MeshedHole found " +
786 0 : std::to_string(n_positive_areas) +
787 : " counter-clockwise boundaries and cannot choose one!");
788 : }
789 :
790 : }
791 142 : else if (outer_edge_type != 0)
792 : {
793 40 : if (n_negative_areas > 1)
794 : {
795 1 : print_areas();
796 27 : report_error("MeshedHole found " +
797 51 : std::to_string(n_negative_areas) +
798 : " clockwise boundaries and cannot choose one!");
799 : }
800 :
801 : }
802 :
803 : // Hey, no errors! Broadcast that empty string.
804 130 : mesh.comm().broadcast(error_reported);
805 130 : mesh.comm().broadcast(_points);
806 130 : mesh.comm().broadcast(_midpoints);
807 871 : }
808 :
809 :
810 3627 : unsigned int TriangulatorInterface::MeshedHole::n_points() const
811 : {
812 3793 : return _points.size();
813 : }
814 :
815 :
816 4364 : unsigned int TriangulatorInterface::MeshedHole::n_midpoints() const
817 : {
818 172 : libmesh_assert (!(_midpoints.size() % _points.size()));
819 4708 : return _midpoints.size() / _points.size();
820 : }
821 :
822 :
823 36709 : Point TriangulatorInterface::MeshedHole::point(const unsigned int n) const
824 : {
825 1504 : libmesh_assert_less (n, _points.size());
826 38213 : return _points[n];
827 : }
828 :
829 :
830 1965 : Point TriangulatorInterface::MeshedHole::midpoint(const unsigned int m,
831 : const unsigned int n) const
832 : {
833 1965 : const unsigned int n_mid = this->n_midpoints();
834 78 : libmesh_assert_less (m, n_mid);
835 78 : libmesh_assert_less (n, _points.size());
836 2043 : return _midpoints[n*n_mid+m];
837 : }
838 :
839 :
840 158 : Point TriangulatorInterface::MeshedHole::inside() const
841 : {
842 : // This is expensive to compute, so only do it when we first need it
843 158 : if (_center(0) == std::numeric_limits<Real>::max())
844 150 : _center = this->calculate_inside_point();
845 :
846 158 : return _center;
847 : }
848 :
849 :
850 :
851 : } // namespace libMesh
|
