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