libMesh
Loading...
Searching...
No Matches
mesh_triangulation.C
Go to the documentation of this file.
1#include <libmesh/boundary_info.h>
2#include <libmesh/elem.h>
3#include <libmesh/mesh.h>
4#include <libmesh/mesh_generation.h>
5#include <libmesh/mesh_modification.h>
6#include <libmesh/mesh_triangle_holes.h>
7#include <libmesh/mesh_triangle_interface.h>
8#include <libmesh/parallel_implementation.h> // max()
9#include <libmesh/parsed_function.h>
10#include <libmesh/point.h>
11#include <libmesh/poly2tri_triangulator.h>
12
13#include "test_comm.h"
14#include "libmesh_cppunit.h"
15
16#include <algorithm>
17#include <cmath>
18#include <regex>
19
20
21using namespace libMesh;
22
23class MeshTriangulationTest : public CppUnit::TestCase
24{
29public:
30 LIBMESH_CPPUNIT_TEST_SUITE( MeshTriangulationTest );
31
32 CPPUNIT_TEST( testTriangleHoleArea );
33 CPPUNIT_TEST( testTriangleHoleContains );
34 CPPUNIT_TEST( testHoleContainsRayDegeneracies );
35
36#ifdef LIBMESH_HAVE_POLY2TRI
37 CPPUNIT_TEST( testPoly2Tri );
38 CPPUNIT_TEST( testPoly2TriHalfDomain );
39 CPPUNIT_TEST( testPoly2TriHalfDomainEdge3 );
40 CPPUNIT_TEST( testPoly2TriInterp );
41 CPPUNIT_TEST( testPoly2TriInterp2 );
42 CPPUNIT_TEST( testPoly2TriHoles );
43 CPPUNIT_TEST( testPoly2TriMeshedHoles );
44 CPPUNIT_TEST( testPoly2TriEdge3ToTri6 );
45 CPPUNIT_TEST( testPoly2TriEdge3ToTri7CenterFixup );
46 CPPUNIT_TEST( testPoly2TriEdge3ToTri6FlipThrows );
47# ifdef LIBMESH_ENABLE_AMR
48 CPPUNIT_TEST( testPoly2TriRoundHole );
49# endif
50 CPPUNIT_TEST( testPoly2TriEdges );
51 CPPUNIT_TEST( testPoly2TriEdge3s );
52 CPPUNIT_TEST( testPoly2TriBadEdges );
53 CPPUNIT_TEST( testPoly2TriBad1DMultiBoundary );
54 CPPUNIT_TEST( testPoly2TriBad2DMultiBoundary );
55 CPPUNIT_TEST( testPoly2TriEdgesRefined );
56 CPPUNIT_TEST( testPoly2TriSegments );
57 CPPUNIT_TEST( testPoly2TriRefined );
58 CPPUNIT_TEST( testPoly2TriNonRefined );
59 CPPUNIT_TEST( testPoly2TriExtraRefined );
60 CPPUNIT_TEST( testPoly2TriHolesRefined );
61 CPPUNIT_TEST( testPoly2TriHolesInterpRefined );
62 CPPUNIT_TEST( testPoly2TriHolesInteriorRefined );
63 CPPUNIT_TEST( testPoly2TriHolesInteriorExtraRefined );
64 // This covers an old poly2tri collinearity-tolerance bug
65 CPPUNIT_TEST( testPoly2TriHolesExtraRefined );
66
67 // These cover more recent tolerance issues when verifying holes
68 CPPUNIT_TEST( testPoly2TriHolePerturbed );
69 CPPUNIT_TEST( testPoly2TriHoleTangentPerturbed );
70
71 CPPUNIT_TEST( testPoly2TriNonUniformRefined );
72 CPPUNIT_TEST( testPoly2TriHolesNonUniformRefined );
73#endif
74
75#ifdef LIBMESH_HAVE_TRIANGLE
76 CPPUNIT_TEST( testTriangle );
77 CPPUNIT_TEST( testTriangleHalfDomain );
78 CPPUNIT_TEST( testTriangleInterp );
79 CPPUNIT_TEST( testTriangleInterp2 );
80 CPPUNIT_TEST( testTriangleHoles );
81 CPPUNIT_TEST( testTriangleMeshedHoles );
82# ifdef LIBMESH_ENABLE_AMR
83 CPPUNIT_TEST( testTriangleRoundHole );
84# endif
85 CPPUNIT_TEST( testTriangleEdges );
86 CPPUNIT_TEST( testTriangleSegments );
87 CPPUNIT_TEST( testTriangleEdge3ToTri6 );
88#endif
89
90 CPPUNIT_TEST_SUITE_END();
91
92public:
93 void setUp() {}
94
95 void tearDown() {}
96
97 void commonSettings (TriangulatorInterface & triangulator)
98 {
99 // Use the point order to define the boundary, because our
100 // Poly2Tri implementation doesn't do convex hulls yet, even when
101 // that would give the same answer.
102 triangulator.triangulation_type() = TriangulatorInterface::PSLG;
103
104 // Don't try to insert points unless we're requested to later
105 triangulator.desired_area() = 1000;
106 triangulator.minimum_angle() = 0;
107 triangulator.smooth_after_generating() = false;
108 triangulator.set_verify_hole_boundaries(true);
109 }
110
111 void testExceptionBase(MeshBase & mesh,
112 TriangulatorInterface & triangulator,
113 const char * re)
114 {
115#ifdef LIBMESH_ENABLE_EXCEPTIONS
116 // We can't just CPPUNIT_ASSERT_THROW, because we want to make
117 // sure we were thrown from the right place with the right error
118 // message!
119 bool threw_desired_exception = false;
120 try {
121 this->testTriangulatorBase(mesh, triangulator);
122 }
123 catch (libMesh::LogicError & e) {
124 std::regex msg_regex(re);
125 CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
126 threw_desired_exception = true;
127 }
128 catch (CppUnit::Exception & e) {
129 throw e;
130 }
131 catch (...) {
132 CPPUNIT_ASSERT_MESSAGE("Unexpected exception type thrown", false);
133 }
134 CPPUNIT_ASSERT(threw_desired_exception);
135#endif
136 }
137
138
139 void testTriangleHoleArea()
140 {
141 LOG_UNIT_TEST;
142
143 // Using center=(1,0), radius=2 for the heck of it
144 Point center{1};
145 Real radius = 2;
146 std::vector<TriangulatorInterface::PolygonHole> polyholes;
147
148 // Line
149 polyholes.emplace_back(center, radius, 2);
150 // Triangle
151 polyholes.emplace_back(center, radius, 3);
152 // Square
153 polyholes.emplace_back(center, radius, 4);
154 // Pentagon
155 polyholes.emplace_back(center, radius, 5);
156 // Hexagon
157 polyholes.emplace_back(center, radius, 6);
158
159 for (int i=0; i != 5; ++i)
160 {
161 const int n_sides = i+2;
162 const TriangulatorInterface::Hole & hole = polyholes[i];
163
164 const Real computed_area = hole.area();
165 const Real theta = pi/n_sides;
166 const Real half_side_length = radius*std::cos(theta);
167 const Real apothem = radius*std::sin(theta);
168 const Real area = n_sides * apothem * half_side_length;
169
170 LIBMESH_ASSERT_FP_EQUAL(computed_area, area, TOLERANCE*TOLERANCE);
171 }
172
173 TriangulatorInterface::ArbitraryHole arbhole {center, {{0,-1},{2,-1},{2,1},{0,2}}};
174 LIBMESH_ASSERT_FP_EQUAL(arbhole.area(), Real(5), TOLERANCE*TOLERANCE);
175
176#ifdef LIBMESH_HAVE_TRIANGLE
177 // Make sure we're compatible with the old naming structure too
178 TriangleInterface::PolygonHole square(center, radius, 4);
179 LIBMESH_ASSERT_FP_EQUAL(square.area(), 2*radius*radius, TOLERANCE*TOLERANCE);
180#endif
181 }
182
183
184 void testTriangleHoleContains()
185 {
186 LOG_UNIT_TEST;
187
188 // Using center=(1,2), radius=2 for the heck of it
189 Point center{1,2};
190 Real radius = 2;
191 std::vector<TriangulatorInterface::PolygonHole> polyholes;
192
193 auto check_corners = [center, radius]
194 (const TriangulatorInterface::Hole & hole,
195 unsigned int np)
196 {
197 CPPUNIT_ASSERT_EQUAL(np, hole.n_points());
198 for (auto i : make_range(np))
199 {
200 const Real theta = i * 2 * libMesh::pi / np;
201 const Real xin = center(0) + radius * .99 * std::cos(theta);
202 const Real xout = center(0) + radius * 1.01 * std::cos(theta);
203 const Real yin = center(1) + radius * .99 * std::sin(theta);
204 const Real yout = center(1) + radius * 1.01 * std::sin(theta);
205
206 CPPUNIT_ASSERT(hole.contains(Point(xin, yin)));
207 CPPUNIT_ASSERT(!hole.contains(Point(xout, yout)));
208 }
209 };
210
211 const TriangulatorInterface::PolygonHole triangle(center, radius, 3);
212 check_corners(triangle, 3);
213 const TriangulatorInterface::PolygonHole diamond(center, radius, 4);
214 check_corners(diamond, 4);
215 const TriangulatorInterface::PolygonHole hexagon(center, radius, 6);
216 check_corners(hexagon, 6);
217
218 TriangulatorInterface::ArbitraryHole jaggy
219 {{1,0}, {{0,-1},{2,-1},{2,1},{1.75,-.5},{1.5,1},{1.25,-.5},
220 {1,1},{.75,-.5},{.5,1},{.25,-.5},{0,1}}};
221
222 CPPUNIT_ASSERT(jaggy.contains({.1,-.3}));
223 CPPUNIT_ASSERT(jaggy.contains({.5,.9}));
224 CPPUNIT_ASSERT(jaggy.contains({.9,-.3}));
225 CPPUNIT_ASSERT(jaggy.contains({1,0}));
226 CPPUNIT_ASSERT(jaggy.contains({1.1,-.4}));
227 CPPUNIT_ASSERT(jaggy.contains({1.5,.9}));
228 CPPUNIT_ASSERT(jaggy.contains({1.9,-.3}));
229
230 CPPUNIT_ASSERT(jaggy.contains({1.9,-.5}));
231 CPPUNIT_ASSERT(jaggy.contains({1.6,-.5}));
232 CPPUNIT_ASSERT(jaggy.contains({1.1,-.5}));
233 CPPUNIT_ASSERT(jaggy.contains({.5,-.5}));
234 CPPUNIT_ASSERT(jaggy.contains({.2,-.5}));
235
236 CPPUNIT_ASSERT(!jaggy.contains({.1,.7}));
237 CPPUNIT_ASSERT(!jaggy.contains({.5,1.1}));
238 CPPUNIT_ASSERT(!jaggy.contains({.9,.7}));
239 CPPUNIT_ASSERT(!jaggy.contains({1,-1.1}));
240 CPPUNIT_ASSERT(!jaggy.contains({1.1,.8}));
241 CPPUNIT_ASSERT(!jaggy.contains({1.5,1.1}));
242 CPPUNIT_ASSERT(!jaggy.contains({1.9,.9}));
243
244 CPPUNIT_ASSERT(!jaggy.contains({1.9,1}));
245 CPPUNIT_ASSERT(!jaggy.contains({1.4,1}));
246 CPPUNIT_ASSERT(!jaggy.contains({.9,1}));
247 CPPUNIT_ASSERT(!jaggy.contains({.4,1}));
248 CPPUNIT_ASSERT(!jaggy.contains({-.2,1}));
249 CPPUNIT_ASSERT(!jaggy.contains({-.2,0}));
250 CPPUNIT_ASSERT(!jaggy.contains({1.2,0}));
251
252 TriangulatorInterface::ArbitraryHole square_jaggy
253 {{1,0}, {{-.25,-1},{2,-1},{2,1},{1.75,1},{1.75,-.5},{1.5,-.5},
254 {1.5,1},{1.25,1},{1.25,-.5},{1,-.5},{1,1},{.75,1},
255 {.75,-.5},{.5,-.5},{.5,1},{.25,1},{.25,-.5},{0,-.5},
256 {0,1},{-.25,1}}};
257
258 CPPUNIT_ASSERT(square_jaggy.contains({-.1,-.3}));
259 CPPUNIT_ASSERT(square_jaggy.contains({.4,.9}));
260 CPPUNIT_ASSERT(square_jaggy.contains({.9,-.3}));
261 CPPUNIT_ASSERT(square_jaggy.contains({.9,0}));
262 CPPUNIT_ASSERT(square_jaggy.contains({1.1,-.6}));
263 CPPUNIT_ASSERT(square_jaggy.contains({1.4,.9}));
264 CPPUNIT_ASSERT(square_jaggy.contains({1.9,-.3}));
265
266 CPPUNIT_ASSERT(!square_jaggy.contains({.1,-.3}));
267 CPPUNIT_ASSERT(!square_jaggy.contains({.6,.9}));
268 CPPUNIT_ASSERT(!square_jaggy.contains({1.1,-.3}));
269 CPPUNIT_ASSERT(!square_jaggy.contains({1.1,0}));
270 CPPUNIT_ASSERT(!square_jaggy.contains({1.1,-1.6}));
271 CPPUNIT_ASSERT(!square_jaggy.contains({1.6,.9}));
272 CPPUNIT_ASSERT(!square_jaggy.contains({2.1,-.3}));
273
274 CPPUNIT_ASSERT(square_jaggy.contains({-.1,-.5}));
275 CPPUNIT_ASSERT(square_jaggy.contains({.3,-.5}));
276 CPPUNIT_ASSERT(square_jaggy.contains({.9,-.5}));
277 CPPUNIT_ASSERT(square_jaggy.contains({1.3,-.5}));
278 CPPUNIT_ASSERT(square_jaggy.contains({1.9,-.5}));
279
280 CPPUNIT_ASSERT(!square_jaggy.contains({-.3,1}));
281 CPPUNIT_ASSERT(!square_jaggy.contains({.2,1}));
282 CPPUNIT_ASSERT(!square_jaggy.contains({.6,1}));
283 CPPUNIT_ASSERT(!square_jaggy.contains({1.1,1}));
284 CPPUNIT_ASSERT(!square_jaggy.contains({1.6,1}));
285 CPPUNIT_ASSERT(!square_jaggy.contains({2.1,1}));
286 }
287
288
289 // The case from libMesh issue #3496, with no triangulator in the way.
290 // contains() casts its ray in the +x direction, so a hole vertex that
291 // very nearly shares the query point's y is the degenerate case.
292 //
293 // Rays that hit a vertex exactly - passing through one, grazing one,
294 // or running along a whole edge - are already covered by the holes in
295 // testTriangleHoleContains. What is special here is the scale: the
296 // ray misses the boundary's vertices by a few ulps, which is where
297 // comparing a parametric edge coordinate against a tolerance used to
298 // classify the same vertex differently from each of its two edges.
299 void testHoleContainsRayDegeneracies()
300 {
301 LOG_UNIT_TEST;
302
303 // Every vertex of a hole has to test as inside a boundary that
304 // encloses it, however nearly collinear the ray from it is with the
305 // boundary's own vertices.
306 const TriangulatorInterface::PolygonHole
307 perturbed(Point(0,4.e-16), std::sqrt(2_R)/2, 4);
308 const TriangulatorInterface::ArbitraryHole outer_bdy
309 {{0,0}, {{100,0},{100,100},{0,100},{-100,100},
310 {-100,0},{-100,-100},{0,-100},{100,-100}}};
311
312 for (auto i : make_range(perturbed.n_points()))
313 CPPUNIT_ASSERT(outer_bdy.contains(perturbed.point(i)));
314 }
315
316
317 void testTriangulatorBase(MeshBase & mesh,
318 TriangulatorInterface & triangulator)
319 {
320 commonSettings(triangulator);
321
322 triangulator.triangulate();
323
324 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(2));
325 for (const auto & elem : mesh.element_ptr_range())
326 {
327 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
328
329 // Make sure we're not getting any inverted elements
330 auto cross_prod =
331 (elem->point(1) - elem->point(0)).cross
332 (elem->point(2) - elem->point(0));
333
334 CPPUNIT_ASSERT_GREATER(Real(0), cross_prod(2));
335
336 bool found_triangle = false;
337 for (const auto & node : elem->node_ref_range())
338 {
339 const Point & point = node;
340 if (point == Point(0,0))
341 {
342 found_triangle = true;
343 CPPUNIT_ASSERT((elem->point(0) == Point(0,0) &&
344 elem->point(1) == Point(1,0) &&
345 elem->point(2) == Point(0,1)) ||
346 (elem->point(1) == Point(0,0) &&
347 elem->point(2) == Point(1,0) &&
348 elem->point(0) == Point(0,1)) ||
349 (elem->point(2) == Point(0,0) &&
350 elem->point(0) == Point(1,0) &&
351 elem->point(1) == Point(0,1)));
352 }
353 if (point == Point(1,2))
354 {
355 found_triangle = true;
356 CPPUNIT_ASSERT((elem->point(0) == Point(0,1) &&
357 elem->point(1) == Point(1,0) &&
358 elem->point(2) == Point(1,2)) ||
359 (elem->point(1) == Point(0,1) &&
360 elem->point(2) == Point(1,0) &&
361 elem->point(0) == Point(1,2)) ||
362 (elem->point(2) == Point(0,1) &&
363 elem->point(0) == Point(1,0) &&
364 elem->point(1) == Point(1,2)));
365 }
366 }
367 CPPUNIT_ASSERT(found_triangle);
368 }
369 }
370
371 void testEdge3ToTri6Base(MeshBase & mesh,
372 TriangulatorInterface & triangulator)
373 {
374 commonSettings(triangulator);
375
376 triangulator.elem_type() = TRI6;
377 triangulator.triangulate();
378 const Real hsq2 = std::sqrt(2) / 2.0;
379
380 // Check element number
381 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), (dof_id_type)2);
382 for (const auto & elem : mesh.element_ptr_range())
383 {
384 // Check element type
385 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI6);
386
387 for (const auto & i_side : make_range(elem->n_sides()))
388 {
389 // We only want to check the sides on the external boundary
390 if (elem->neighbor_ptr(i_side) == nullptr)
391 {
392 const Point & side_pt_0 = *(elem->node_ptr(i_side));
393 const Point & side_pt_1 = *(elem->node_ptr((i_side + 1) % 3));
394 const Point & side_pt_2 = *(elem->node_ptr(i_side + 3));
395 const bool x_sign = (side_pt_0(0) == 1 || side_pt_1(0) == 1);
396 const bool y_sign = (side_pt_0(1) == 1 || side_pt_1(1) == 1);
397 const Point ref_side_pt_2 = Point (x_sign ? hsq2 : -hsq2, y_sign ? hsq2 : -hsq2);
398 CPPUNIT_ASSERT_EQUAL(ref_side_pt_2, side_pt_2);
399 }
400 }
401 }
402 }
403
404 void testTriangulator(MeshBase & mesh,
405 TriangulatorInterface & triangulator)
406 {
407 // A non-square quad, so we don't have ambiguity about which
408 // diagonal a Delaunay algorithm will pick.
409 // Manually-numbered points, so we can use the point numbering as
410 // a segment ordering even on DistributedMesh.
411 mesh.add_point(Point(0,0), 0);
412 mesh.add_point(Point(1,0), 1);
413 mesh.add_point(Point(1,2), 2);
414 mesh.add_point(Point(0,1), 3);
415
416 this->testTriangulatorBase(mesh, triangulator);
417 }
418
419
420 void testTriangulatorTrapMesh (UnstructuredMesh & mesh)
421 {
422 // A non-square quad, so we don't have ambiguity about which
423 // diagonal a Delaunay algorithm will pick.
424 // Manually-numbered points, so we can use the point numbering as
425 // a segment ordering even on DistributedMesh.
426 mesh.add_point(Point(0,0), 0);
427 mesh.add_point(Point(1,0), 1);
428 mesh.add_point(Point(1,2), 2);
429 mesh.add_point(Point(0,1), 3);
430 }
431
432
433 void testTriangulatorInterp(UnstructuredMesh & mesh,
434 TriangulatorInterface & triangulator,
435 int interpolate_boundary_points,
436 dof_id_type n_expected_elem,
437 Real expected_total_area = 1.5,
438 Real desired_area = 1000)
439 {
440 commonSettings(triangulator);
441
442 if (!mesh.n_nodes())
443 testTriangulatorTrapMesh(mesh);
444
445 // Interpolate points!
446 triangulator.set_interpolate_boundary_points(interpolate_boundary_points);
447
448 // Try to insert points?
449 triangulator.desired_area() = desired_area;
450
451 triangulator.triangulate();
452
453 if (n_expected_elem)
454 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_expected_elem);
455
456 Real area = 0;
457 for (const auto & elem : mesh.active_local_element_ptr_range())
458 {
459 CPPUNIT_ASSERT_EQUAL(elem->level(), 0u);
460 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
461
462 area += elem->volume();
463 }
464
465 mesh.comm().sum(area);
466
467 LIBMESH_ASSERT_FP_EQUAL(area, expected_total_area, TOLERANCE*TOLERANCE);
468 }
469
470
471 void testFoundCenters(const MeshBase & mesh,
472 const std::vector<Point> & expected_centers)
473 {
474 std::vector<bool> found_centers(expected_centers.size(), false);
475
476 for (const auto & elem : mesh.element_ptr_range())
477 {
478 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
479
480 // Make sure we're not getting any inverted elements
481 auto cross_prod =
482 (elem->point(1) - elem->point(0)).cross
483 (elem->point(2) - elem->point(0));
484
485 CPPUNIT_ASSERT_GREATER(Real(0), cross_prod(2));
486
487 // Make sure we're finding all the elements we expect
488 Point center = elem->vertex_average();
489
490 bool found_mine = false;
491 for (auto i : index_range(expected_centers))
492 {
493 Point possible = expected_centers[i];
494
495 if (possible.absolute_fuzzy_equals(center, TOLERANCE*TOLERANCE))
496 {
497 found_mine = true;
498 found_centers[i] = true;
499 }
500 }
501 CPPUNIT_ASSERT(found_mine);
502 }
503
504 mesh.comm().max(found_centers);
505
506 for (auto found_it : found_centers)
507 CPPUNIT_ASSERT(found_it);
508 }
509
510
511 void testTriangulatorHoles(MeshBase & mesh,
512 TriangulatorInterface & triangulator)
513 {
514 // A square quad; we'll put a diamond hole in the middle to make
515 // the Delaunay selection unambiguous.
516 mesh.add_point(Point(-1,-1), 0);
517 mesh.add_point(Point(1,-1), 1);
518 mesh.add_point(Point(1,1), 2);
519 mesh.add_point(Point(-1,1), 3);
520
521 commonSettings(triangulator);
522
523 // Add a diamond hole in the center
524 TriangulatorInterface::PolygonHole diamond(Point(0), std::sqrt(2_R)/2, 4);
525 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
526 triangulator.attach_hole_list(&holes);
527
528 triangulator.triangulate();
529
530 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(8));
531
532 // Center coordinates for all the elements we expect
533 Real r2p2o6 = (std::sqrt(2_R)+2)/6;
534 Real r2p4o6 = (std::sqrt(2_R)+4)/6;
535
536 std::vector <Point> expected_centers
537 { {r2p2o6,r2p2o6}, {r2p2o6,-r2p2o6},
538 {-r2p2o6,r2p2o6}, {-r2p2o6,-r2p2o6},
539 {0,r2p4o6}, {r2p4o6, 0},
540 {0,-r2p4o6}, {-r2p4o6, 0}
541 };
542
543 testFoundCenters(mesh, expected_centers);
544 }
545
546
547 void testTriangulatorHolePerturbed(MeshBase & mesh,
548 TriangulatorInterface & triangulator)
549 {
550 // Points based on a simplification of a hole verification failure
551 // case
552 mesh.add_point(Point(100,0), 0);
553 mesh.add_point(Point(100,100), 1);
554 mesh.add_point(Point(0,100), 2);
555 mesh.add_point(Point(-100,100), 3);
556 mesh.add_point(Point(-100,0), 4);
557 mesh.add_point(Point(-100,-100), 5);
558 mesh.add_point(Point(0,-100), 6);
559 mesh.add_point(Point(100,-100), 7);
560
561 commonSettings(triangulator);
562
563 // Don't insert points here; we want the triangulation determined
564 // purely by the boundary and hole points
565 triangulator.desired_area() = 1e16;
566
567 // Add a diamond hole in *almost* the center; the tiny perturbation
568 // puts a hole vertex ray-collinear with a boundary vertex, which
569 // used to trip up the hole verification
570 TriangulatorInterface::PolygonHole
571 diamond(Point(0,4.e-16),
572 std::sqrt(2_R)/2, 4);
573 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
574 triangulator.attach_hole_list(&holes);
575
576 triangulator.triangulate();
577
578 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(12));
579
580 // Center coordinates for all the elements we expect
581 const Real r2p200o6 = (std::sqrt(2_R)+200)/6,
582 r2p400o6 = (std::sqrt(2_R)+400)/6;
583
584 std::vector <Point> expected_centers
585 { {r2p400o6,100./3}, {-r2p400o6,100./3},
586 {r2p400o6,-100./3}, {-r2p400o6,-100./3},
587 {100./3,r2p400o6}, {-100./3,r2p400o6},
588 {100./3,-r2p400o6}, {-100./3,-r2p400o6},
589 {r2p200o6,r2p200o6}, {-r2p200o6,r2p200o6},
590 {r2p200o6,-r2p200o6}, {-r2p200o6,-r2p200o6},
591 };
592
593 testFoundCenters(mesh, expected_centers);
594 }
595
596
597 void testTriangulatorHoleTangentPerturbed(MeshBase & mesh,
598 TriangulatorInterface & triangulator)
599 {
600 // Points based on a simplification of a hole verification failure
601 // case
602 mesh.add_point(Point(200,0), 0);
603 mesh.add_point(Point(200,100), 1);
604 mesh.add_point(Point(100,100), 2);
605 mesh.add_point(Point(0,100), 3);
606 mesh.add_point(Point(-100,100), 4);
607 mesh.add_point(Point(-200,100), 5);
608 mesh.add_point(Point(-200,0), 6);
609 mesh.add_point(Point(-200,-100), 7);
610 mesh.add_point(Point(-100,-100), 8);
611 mesh.add_point(Point(0,-100), 9);
612 mesh.add_point(Point(100,-100), 10);
613 mesh.add_point(Point(200,-100), 11);
614
615 commonSettings(triangulator);
616
617 // Don't insert points here; we want the triangulation determined
618 // purely by the boundary and hole points
619 triangulator.desired_area() = 1e16;
620
621 // Two diamond holes, in *almost* the center of each half of that
622 // rectangle. The tiny perturbations make a ray from one hole
623 // vertex "tangent" to a vertex of the other hole, which used to
624 // trip up the hole verification double-intersection check.
625 TriangulatorInterface::PolygonHole
626 left_diamond(Point(-100,4.e-16),
627 std::sqrt(2_R)/2, 4),
628 right_diamond(Point(100,-4.e-16),
629 std::sqrt(2_R)/2, 4);
630 const std::vector<TriangulatorInterface::Hole*> holes
631 { &left_diamond, &right_diamond };
632 triangulator.attach_hole_list(&holes);
633
634 triangulator.triangulate();
635
636 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(22));
637
638 // Center coordinates for all the elements we expect
639 const Real r2p200o6 = (std::sqrt(2_R)+200)/6,
640 r2p400o6 = (std::sqrt(2_R)+400)/6;
641
642 std::vector <Point> expected_centers
643 { // Right half, around the diamond hole near x = 100
644 {100+r2p400o6,100./3}, {100+r2p400o6,-100./3},
645 {100+100./3,r2p400o6}, {100-100./3,r2p400o6},
646 {100+100./3,-r2p400o6}, {100-100./3,-r2p400o6},
647 {100+r2p200o6,r2p200o6}, {100-r2p200o6,r2p200o6},
648 {100+r2p200o6,-r2p200o6}, {100-r2p200o6,-r2p200o6},
649 // Left half, mirror of the right around the hole near x = -100
650 {-100-r2p400o6,100./3}, {-100-r2p400o6,-100./3},
651 {-100-100./3,r2p400o6}, {-100+100./3,r2p400o6},
652 {-100-100./3,-r2p400o6}, {-100+100./3,-r2p400o6},
653 {-100-r2p200o6,r2p200o6}, {-100+r2p200o6,r2p200o6},
654 {-100-r2p200o6,-r2p200o6}, {-100+r2p200o6,-r2p200o6},
655 // The two elements meeting the shared wall at x = 0
656 {0,100./3}, {0,-100./3}
657 };
658
659 testFoundCenters(mesh, expected_centers);
660 }
661
662
663 void testTriangulatorMeshedHoles(MeshBase & mesh,
664 TriangulatorInterface & triangulator)
665 {
666 // A square quad; we'll put a square hole in the middle. Offset
667 // this to catch a potential bug I missed on the first try.
668 mesh.add_point(Point(19,19), 0);
669 mesh.add_point(Point(21,19), 1);
670 mesh.add_point(Point(21,21), 2);
671 mesh.add_point(Point(19,21), 3);
672
673 commonSettings(triangulator);
674
675 // Add a square meshed hole in the center
676 Mesh centermesh { mesh.comm() };
677 MeshTools::Generation::build_square (centermesh, 2, 2, 19.5, 20.5, 19.5, 20.5, QUAD4);
678
679 TriangulatorInterface::MeshedHole centerhole { centermesh };
680
681 CPPUNIT_ASSERT_EQUAL(centerhole.n_points(), 8u);
682 CPPUNIT_ASSERT_EQUAL(centerhole.area(), Real(1));
683 Point inside = centerhole.inside();
684 CPPUNIT_ASSERT_EQUAL(inside(0), Real(20));
685 CPPUNIT_ASSERT_EQUAL(inside(1), Real(20));
686
687 const std::vector<TriangulatorInterface::Hole*> holes { &centerhole };
688 triangulator.attach_hole_list(&holes);
689
690 triangulator.triangulate();
691
692 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(12));
693
694 // Center coordinates for all the elements we expect
695 std::vector <Point> expected_centers
696 { {-0.5, Real(2)/3}, {0, Real(5)/6},
697 {0.5, Real(2)/3},
698 {Real(2)/3, -0.5}, {Real(5)/6, 0},
699 {Real(2)/3, 0.5},
700 {-0.5, -Real(2)/3}, {0, -Real(5)/6},
701 {0.5, -Real(2)/3},
702 {-Real(2)/3, -0.5}, {-Real(5)/6, 0},
703 {-Real(2)/3, 0.5} };
704
705 // With the same offset
706 for (auto & p : expected_centers)
707 p += Point(20,20);
708
709 testFoundCenters(mesh, expected_centers);
710 }
711
712
713#ifdef LIBMESH_ENABLE_AMR
714 void testTriangulatorRoundHole(MeshBase & mesh,
715 TriangulatorInterface & triangulator)
716 {
717 // A square quad; we'll put a round hole in the middle.
718 mesh.add_point(Point(19,19), 0);
719 mesh.add_point(Point(21,19), 1);
720 mesh.add_point(Point(21,21), 2);
721 mesh.add_point(Point(19,21), 3);
722
723 commonSettings(triangulator);
724
725 // Add a square meshed hole in the center
726 const Real radius = 0.5;
727 const Point center{20,20};
728
729 Mesh centermesh { mesh.comm() };
730 MeshTools::Generation::build_sphere(centermesh, radius, 1, QUAD9, 0);
731 MeshTools::Modification::translate(centermesh, center(0), center(1));
732
733 TriangulatorInterface::MeshedHole centerhole { centermesh };
734
735 CPPUNIT_ASSERT_EQUAL(centerhole.n_points(), 8u);
736 Point inside = centerhole.inside();
737 CPPUNIT_ASSERT_EQUAL(inside(0), Real(20));
738 CPPUNIT_ASSERT_EQUAL(inside(1), Real(20));
739
740 const std::vector<TriangulatorInterface::Hole*> holes { &centerhole };
741 triangulator.attach_hole_list(&holes);
742 triangulator.elem_type() = TRI6;
743
744 triangulator.triangulate();
745
746 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(12));
747 CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), static_cast<dof_id_type>(36));
748
749 // Make sure we didn't screw up the outer sides. We should
750 // have exact values for the outer vertices, so we can use
751 // those for a map.
752 std::map<std::pair<Real, Real>, Point> outer_midpoints
753 {{{19, 19}, {20, 19}},
754 {{21, 19}, {21, 20}},
755 {{21, 21}, {20, 21}},
756 {{19, 21}, {19, 20}}
757 };
758
759 for (const auto & elem : mesh.active_local_element_ptr_range())
760 for (const auto n : make_range(elem->n_sides()))
761 {
762 if (elem->neighbor_ptr(n))
763 continue;
764
765 auto it =
766 outer_midpoints.find(std::make_pair(elem->point(n)(0),
767 elem->point(n)(1)));
768 if (it != outer_midpoints.end())
769 {
770 const Point error = it->second - elem->point(n+3);
771 CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
772 error.norm_sq());
773 }
774 else
775 {
776 const Point radius1 = elem->point(n) - center;
777 CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
778 std::abs(radius1.norm()-radius));
779
780 const Point radius2 = elem->point((n+1)%3) - center;
781 CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
782 std::abs(radius2.norm()-radius));
783
784 const Point radius3 = elem->point(n+3) - center;
785 CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
786 std::abs(radius3.norm()-radius));
787 }
788 }
789 }
790#endif // LIBMESH_ENABLE_AMR
791
792
793 void testTriangulatorSegments(MeshBase & mesh,
794 TriangulatorInterface & triangulator)
795 {
796 // The same quad as testTriangulator, but out of order
797 mesh.add_point(Point(0,0), 0);
798 mesh.add_point(Point(1,2), 1);
799 mesh.add_point(Point(1,0), 2);
800 mesh.add_point(Point(0,1), 3);
801
802 // Segments to put them in order
803 triangulator.segments = {{0,2},{2,1},{1,3},{3,0}};
804
805 this->testTriangulatorBase(mesh, triangulator);
806 }
807
808
809 void testEdgesMesh(MeshBase & mesh, ElemType elem_type)
810 {
811 // The same quad as testTriangulator, but out of order
812 auto node0 = mesh.add_point(Point(0,0), 0);
813 auto node1 = mesh.add_point(Point(1,2), 1);
814 auto node2 = mesh.add_point(Point(1,0), 2);
815 auto node3 = mesh.add_point(Point(0,1), 3);
816
817 // Edges, also out of order, but enough to put them in order
818 auto edge13 = mesh.add_elem(Elem::build(elem_type));
819 edge13->set_node(0, node1);
820 edge13->set_node(1, node3);
821 auto edge02 = mesh.add_elem(Elem::build(elem_type));
822 edge02->set_node(0, node0);
823 edge02->set_node(1, node2);
824 auto edge30 = mesh.add_elem(Elem::build(elem_type));
825 edge30->set_node(0, node3);
826 edge30->set_node(1, node0);
827 auto edge21 = mesh.add_elem(Elem::build(elem_type));
828 edge21->set_node(0, node2);
829 edge21->set_node(1, node1);
830
831 // Add mid-edge nodes if asked to
832 if (elem_type == EDGE3)
833 {
834 auto node4 = mesh.add_point(Point(.5,1.5), 4);
835 edge13->set_node(2, node4);
836 auto node5 = mesh.add_point(Point(.5,0), 5);
837 edge02->set_node(2, node5);
838 auto node6 = mesh.add_point(Point(0,.5), 6);
839 edge30->set_node(2, node6);
840 auto node7 = mesh.add_point(Point(1,1), 7);
841 edge21->set_node(2, node7);
842 }
843 else
844 libmesh_assert(elem_type == EDGE2);
845
846 mesh.prepare_for_use();
847 }
848
849 void testEdge3Mesh(MeshBase & mesh)
850 {
851 const Real hsq2 = std::sqrt(2) / 2.0;
852 auto node0 = mesh.add_point(Point(1,0), 0);
853 auto node1 = mesh.add_point(Point(hsq2,hsq2), 1);
854 auto node2 = mesh.add_point(Point(0,1), 2);
855 auto node3 = mesh.add_point(Point(-hsq2,hsq2), 3);
856 auto node4 = mesh.add_point(Point(-1,0), 4);
857 auto node5 = mesh.add_point(Point(-hsq2,-hsq2), 5);
858 auto node6 = mesh.add_point(Point(0,-1), 6);
859 auto node7 = mesh.add_point(Point(hsq2,-hsq2), 7);
860
861 auto edge021 = mesh.add_elem(Elem::build(EDGE3));
862 edge021->set_node(0, node0);
863 edge021->set_node(1, node2);
864 edge021->set_node(2, node1);
865 auto edge243 = mesh.add_elem(Elem::build(EDGE3));
866 edge243->set_node(0, node2);
867 edge243->set_node(1, node4);
868 edge243->set_node(2, node3);
869 auto edge465 = mesh.add_elem(Elem::build(EDGE3));
870 edge465->set_node(0, node4);
871 edge465->set_node(1, node6);
872 edge465->set_node(2, node5);
873 auto edge607 = mesh.add_elem(Elem::build(EDGE3));
874 edge607->set_node(0, node6);
875 edge607->set_node(1, node0);
876 edge607->set_node(2, node7);
877
878 mesh.prepare_for_use();
879 }
880
881
882 void testHalfDomain(MeshBase & mesh,
883 TriangulatorInterface & triangulator,
884 ElemType elem_type)
885 {
886 // We might have 4 or 8 nodes on the outer boundary, depending on
887 // whether it has mid-edge nodes
888 const dof_id_type off = (elem_type == EDGE3)*4;
889
890 // A pentagon we'll avoid via subdomain ids
891 auto node4 = mesh.add_point(Point(2,0), 4+off);
892 auto node5 = mesh.add_point(Point(3,0), 5+off);
893 auto node6 = mesh.add_point(Point(3,2), 6+off);
894 auto node7 = mesh.add_point(Point(2,2), 7+off);
895 auto node8 = mesh.add_point(Point(2,1), 8+off);
896
897 auto edge45 = mesh.add_elem(Elem::build(elem_type));
898 edge45->set_node(0, node4);
899 edge45->set_node(1, node5);
900 edge45->subdomain_id() = 1;
901 auto edge56 = mesh.add_elem(Elem::build(elem_type));
902 edge56->set_node(0, node5);
903 edge56->set_node(1, node6);
904 edge56->subdomain_id() = 1;
905 auto edge67 = mesh.add_elem(Elem::build(elem_type));
906 edge67->set_node(0, node6);
907 edge67->set_node(1, node7);
908 edge67->subdomain_id() = 1;
909 auto edge78 = mesh.add_elem(Elem::build(elem_type));
910 edge78->set_node(0, node7);
911 edge78->set_node(1, node8);
912 edge78->subdomain_id() = 1;
913 auto edge84 = mesh.add_elem(Elem::build(elem_type));
914 edge84->set_node(0, node8);
915 edge84->set_node(1, node4);
916 edge84->subdomain_id() = 1;
917
918 if (elem_type == EDGE3)
919 {
920 auto node9 = mesh.add_point(Point(2.5,0), 9+off);
921 edge45->set_node(2, node9);
922 auto node10 = mesh.add_point(Point(3,1), 10+off);
923 edge56->set_node(2, node10);
924 auto node11 = mesh.add_point(Point(2.5,2), 11+off);
925 edge67->set_node(2, node11);
926 auto node12 = mesh.add_point(Point(2,1.5), 12+off);
927 edge78->set_node(2, node12);
928 auto node13 = mesh.add_point(Point(2,.5), 13+off);
929 edge84->set_node(2, node13);
930 }
931 else
932 libmesh_assert(elem_type == EDGE2);
933
934 testEdgesMesh(mesh, elem_type);
935
936 std::set<std::size_t> bdy_ids {0};
937 triangulator.set_outer_boundary_ids(bdy_ids);
938
939 this->testTriangulatorBase(mesh, triangulator);
940 }
941
942
943#ifdef LIBMESH_HAVE_TRIANGLE
944 void testTriangle()
945 {
946 LOG_UNIT_TEST;
947
948 Mesh mesh(*TestCommWorld);
949 TriangleInterface triangle(mesh);
950 testTriangulator(mesh, triangle);
951 }
952
953
954 void testTriangleHalfDomain()
955 {
956 LOG_UNIT_TEST;
957
958 Mesh mesh(*TestCommWorld);
959 TriangleInterface triangle(mesh);
960 testHalfDomain(mesh, triangle, EDGE2);
961 }
962
963
964 void testTriangleInterp()
965 {
966 LOG_UNIT_TEST;
967
968 Mesh mesh(*TestCommWorld);
969 TriangleInterface triangle(mesh);
970 testTriangulatorInterp(mesh, triangle, 1, 6);
971 }
972
973
974 void testTriangleInterp2()
975 {
976 LOG_UNIT_TEST;
977
978 Mesh mesh(*TestCommWorld);
979 TriangleInterface triangle(mesh);
980 testTriangulatorInterp(mesh, triangle, 2, 10);
981 }
982
983
984 void testTriangleHoles()
985 {
986 LOG_UNIT_TEST;
987
988 Mesh mesh(*TestCommWorld);
989 TriangleInterface triangle(mesh);
990 testTriangulatorHoles(mesh, triangle);
991 }
992
993
994 void testTriangleMeshedHoles()
995 {
996 LOG_UNIT_TEST;
997
998 Mesh mesh(*TestCommWorld);
999 TriangleInterface triangle(mesh);
1000 testTriangulatorMeshedHoles(mesh, triangle);
1001 }
1002
1003
1004#ifdef LIBMESH_ENABLE_AMR
1005 void testTriangleRoundHole()
1006 {
1007 LOG_UNIT_TEST;
1008
1009 Mesh mesh(*TestCommWorld);
1010 TriangleInterface triangle(mesh);
1011 testTriangulatorRoundHole(mesh, triangle);
1012 }
1013#endif // LIBMESH_ENABLE_AMR
1014
1015
1016 void testTriangleEdges()
1017 {
1018 LOG_UNIT_TEST;
1019
1020 Mesh mesh(*TestCommWorld);
1021 TriangleInterface triangulator(mesh);
1022 testEdgesMesh(mesh, EDGE2);
1023
1024 this->testTriangulatorBase(mesh, triangulator);
1025 }
1026
1027 void testTriangleSegments()
1028 {
1029 LOG_UNIT_TEST;
1030
1031 Mesh mesh(*TestCommWorld);
1032 TriangleInterface triangle(mesh);
1033 testTriangulatorSegments(mesh, triangle);
1034 }
1035
1036 void testTriangleEdge3ToTri6()
1037 {
1038 LOG_UNIT_TEST;
1039
1040 Mesh mesh(*TestCommWorld);
1041 TriangleInterface triangulator(mesh);
1042
1043 testEdge3Mesh(mesh);
1044
1045 this->testEdge3ToTri6Base(mesh, triangulator);
1046 }
1047
1048#endif // LIBMESH_HAVE_TRIANGLE
1049
1050
1051#ifdef LIBMESH_HAVE_POLY2TRI
1052 void testPoly2Tri()
1053 {
1054 LOG_UNIT_TEST;
1055
1056 Mesh mesh(*TestCommWorld);
1057 Poly2TriTriangulator p2t_tri(mesh);
1058 testTriangulator(mesh, p2t_tri);
1059 }
1060
1061
1062 void testPoly2TriHalfDomain()
1063 {
1064 LOG_UNIT_TEST;
1065
1066 Mesh mesh(*TestCommWorld);
1067 Poly2TriTriangulator p2t_tri(mesh);
1068 testHalfDomain(mesh, p2t_tri, EDGE2);
1069 }
1070
1071
1072 void testPoly2TriHalfDomainEdge3()
1073 {
1074 LOG_UNIT_TEST;
1075
1076 Mesh mesh(*TestCommWorld);
1077 Poly2TriTriangulator p2t_tri(mesh);
1078 testHalfDomain(mesh, p2t_tri, EDGE3);
1079 }
1080
1081
1082 void testPoly2TriInterp()
1083 {
1084 LOG_UNIT_TEST;
1085
1086 Mesh mesh(*TestCommWorld);
1087 Poly2TriTriangulator p2t_tri(mesh);
1088 testTriangulatorInterp(mesh, p2t_tri, 1, 6);
1089 }
1090
1091
1092 void testPoly2TriInterp2()
1093 {
1094 LOG_UNIT_TEST;
1095
1096 Mesh mesh(*TestCommWorld);
1097 Poly2TriTriangulator p2t_tri(mesh);
1098 testTriangulatorInterp(mesh, p2t_tri, 2, 10);
1099 }
1100
1101
1102 void testPoly2TriHoles()
1103 {
1104 LOG_UNIT_TEST;
1105
1106 Mesh mesh(*TestCommWorld);
1107 Poly2TriTriangulator p2t_tri(mesh);
1108 testTriangulatorHoles(mesh, p2t_tri);
1109 }
1110
1111
1112 void testPoly2TriHolePerturbed()
1113 {
1114 LOG_UNIT_TEST;
1115
1116 Mesh mesh(*TestCommWorld);
1117 Poly2TriTriangulator p2t_tri(mesh);
1118 testTriangulatorHolePerturbed(mesh, p2t_tri);
1119 }
1120
1121
1122 void testPoly2TriHoleTangentPerturbed()
1123 {
1124 LOG_UNIT_TEST;
1125
1126 Mesh mesh(*TestCommWorld);
1127 Poly2TriTriangulator p2t_tri(mesh);
1128 testTriangulatorHoleTangentPerturbed(mesh, p2t_tri);
1129 }
1130
1131
1132 void testPoly2TriMeshedHoles()
1133 {
1134 LOG_UNIT_TEST;
1135
1136 Mesh mesh(*TestCommWorld);
1137 Poly2TriTriangulator p2t_tri(mesh);
1138 testTriangulatorMeshedHoles(mesh, p2t_tri);
1139 }
1140
1141 void testPoly2TriEdge3ToTri6()
1142 {
1143 LOG_UNIT_TEST;
1144
1145 Mesh mesh(*TestCommWorld);
1146 Poly2TriTriangulator triangulator(mesh);
1147
1148 testEdge3Mesh(mesh);
1149
1150 this->testEdge3ToTri6Base(mesh, triangulator);
1151 }
1152
1153
1154 void testPoly2TriEdge3ToTri7CenterFixup()
1155 {
1156 LOG_UNIT_TEST;
1157
1158 Mesh mesh(*TestCommWorld);
1159 Poly2TriTriangulator triangulator(mesh);
1160
1161 // Re-use the curved (unit-circle) Edge3 boundary from testEdge3Mesh:
1162 // four arcs whose midpoints lie on the unit circle. After
1163 // triangulation with TRI7 + boundary-midpoint snap, every TRI7
1164 // element's interior node should land at the curved Tri6 image of
1165 // the reference centroid (1/3, 1/3) -- the formula applied by
1166 // TriangulatorInterface::fixup_tri7_center_nodes() -- *not* at the
1167 // straight-edge vertex centroid (P0+P1+P2)/3 that
1168 // MeshBase::all_complete_order() initially placed there.
1169 testEdge3Mesh(mesh);
1170
1171 commonSettings(triangulator);
1172 triangulator.elem_type() = TRI7;
1173 // The interior-node relocation is opt-in (off by default), so
1174 // enable it here -- this is the behavior under test.
1175 triangulator.set_fixup_tri7_center_nodes(true);
1176 triangulator.triangulate();
1177
1178 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), (dof_id_type)2);
1179
1180 bool saw_curved_element = false;
1181 for (const auto & elem : mesh.element_ptr_range())
1182 {
1183 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI7);
1184
1185 const Point & x0 = elem->point(0);
1186 const Point & x1 = elem->point(1);
1187 const Point & x2 = elem->point(2);
1188 const Point & x3 = elem->point(3);
1189 const Point & x4 = elem->point(4);
1190 const Point & x5 = elem->point(5);
1191 const Point & x6 = elem->point(6);
1192
1193 const Point expected =
1194 (Real(-1)/9) * (x0 + x1 + x2) +
1195 (Real( 4)/9) * (x3 + x4 + x5);
1196 const Point straight = (x0 + x1 + x2) / Real(3);
1197
1198 // Centroid must match the curved-mapping formula exactly.
1199 CPPUNIT_ASSERT(x6.relative_fuzzy_equals(expected, TOLERANCE));
1200
1201 // On the inscribed-square triangulation of the unit circle
1202 // each triangle has at least one curved external side, so the
1203 // curved-mapping centroid is strictly different from the
1204 // straight-edge vertex average -- otherwise the fixup would
1205 // be a no-op and we wouldn't actually be testing it.
1206 if (!x6.relative_fuzzy_equals(straight, TOLERANCE))
1207 saw_curved_element = true;
1208 }
1209
1210 mesh.comm().max(saw_curved_element);
1211 CPPUNIT_ASSERT(saw_curved_element);
1212 }
1213
1214
1215 void testPoly2TriEdge3ToTri6FlipThrows()
1216 {
1217#ifdef LIBMESH_ENABLE_EXCEPTIONS
1218 LOG_UNIT_TEST;
1219
1220 Mesh mesh(*TestCommWorld);
1221
1222 // Build a triangular domain (0,0)-(4,0)-(2,2) using three Edge3
1223 // segments. The bottom segment's mid-edge node is placed at
1224 // (2, 1.95) -- almost on top of the apex (2, 2) -- so that
1225 // snapping the resulting TRI6 mid-edge node to that location
1226 // tangles the element (the curved bottom edge folds across the
1227 // apex). TriangulatorInterface::verify_quadratic_elements()
1228 // should detect the non-positive Jacobian and throw.
1229 auto vA = mesh.add_point(Point(0, 0), 0);
1230 auto vB = mesh.add_point(Point(4, 0), 1);
1231 auto vC = mesh.add_point(Point(2, 2), 2);
1232 auto mAB = mesh.add_point(Point(2, 1.95), 3); // pathological
1233 auto mBC = mesh.add_point(Point(3, 1), 4);
1234 auto mCA = mesh.add_point(Point(1, 1), 5);
1235
1236 auto edgeAB = mesh.add_elem(Elem::build(EDGE3));
1237 edgeAB->set_node(0, vA);
1238 edgeAB->set_node(1, vB);
1239 edgeAB->set_node(2, mAB);
1240
1241 auto edgeBC = mesh.add_elem(Elem::build(EDGE3));
1242 edgeBC->set_node(0, vB);
1243 edgeBC->set_node(1, vC);
1244 edgeBC->set_node(2, mBC);
1245
1246 auto edgeCA = mesh.add_elem(Elem::build(EDGE3));
1247 edgeCA->set_node(0, vC);
1248 edgeCA->set_node(1, vA);
1249 edgeCA->set_node(2, mCA);
1250
1251 mesh.prepare_for_use();
1252
1253 Poly2TriTriangulator triangulator(mesh);
1254 commonSettings(triangulator);
1255 triangulator.elem_type() = TRI6;
1256
1257 bool threw_desired_exception = false;
1258 try {
1259 triangulator.triangulate();
1260 }
1261 catch (libMesh::LogicError & e) {
1262 const std::regex msg_regex("tangled quadratic triangle");
1263 CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
1264 threw_desired_exception = true;
1265 }
1266 catch (CppUnit::Exception & e) {
1267 throw e;
1268 }
1269 catch (...) {
1270 CPPUNIT_ASSERT_MESSAGE("Unexpected exception type thrown", false);
1271 }
1272 CPPUNIT_ASSERT(threw_desired_exception);
1273#endif
1274 }
1275
1276
1277#ifdef LIBMESH_ENABLE_AMR
1278 void testPoly2TriRoundHole()
1279 {
1280 LOG_UNIT_TEST;
1281
1282 Mesh mesh(*TestCommWorld);
1283 Poly2TriTriangulator p2t_tri(mesh);
1284 testTriangulatorRoundHole(mesh, p2t_tri);
1285 }
1286#endif // LIBMESH_ENABLE_AMR
1287
1288
1289 void testPoly2TriEdges()
1290 {
1291 LOG_UNIT_TEST;
1292
1293 Mesh mesh(*TestCommWorld);
1294 Poly2TriTriangulator triangulator(mesh);
1295 testEdgesMesh(mesh, EDGE2);
1296
1297 this->testTriangulatorBase(mesh, triangulator);
1298 }
1299
1300
1301 void testPoly2TriEdge3s()
1302 {
1303 LOG_UNIT_TEST;
1304
1305 Mesh mesh(*TestCommWorld);
1306 Poly2TriTriangulator triangulator(mesh);
1307 testEdgesMesh(mesh, EDGE3);
1308
1309 this->testTriangulatorBase(mesh, triangulator);
1310 }
1311
1312
1313 void testPoly2TriBadEdges()
1314 {
1315 LOG_UNIT_TEST;
1316
1317 Mesh mesh(*TestCommWorld);
1318 Poly2TriTriangulator triangulator(mesh);
1319
1320 // The same quad as testTriangulator, but out of order
1321 auto node0 = mesh.add_point(Point(0,0), 0);
1322 auto node1 = mesh.add_point(Point(1,2), 1);
1323 auto node2 = mesh.add_point(Point(1,0), 2);
1324 auto node3 = mesh.add_point(Point(0,1), 3);
1325
1326 // Edges, but not enough to complete the quad
1327 auto edge13 = mesh.add_elem(Elem::build(EDGE2));
1328 edge13->set_node(0, node1);
1329 edge13->set_node(1, node3);
1330 auto edge02 = mesh.add_elem(Elem::build(EDGE2));
1331 edge02->set_node(0, node0);
1332 edge02->set_node(1, node2);
1333 auto edge30 = mesh.add_elem(Elem::build(EDGE2));
1334 edge30->set_node(0, node3);
1335 edge30->set_node(1, node0);
1336
1337 mesh.prepare_for_use();
1338
1339 testExceptionBase(mesh, triangulator, "Bad edge topology");
1340 }
1341
1342
1343 void testPoly2TriBad1DMultiBoundary()
1344 {
1345 LOG_UNIT_TEST;
1346
1347 Mesh mesh(*TestCommWorld);
1348 Poly2TriTriangulator triangulator(mesh);
1349
1350 // Two separate triangles
1351 auto node0 = mesh.add_point(Point(0,0), 0);
1352 auto node1 = mesh.add_point(Point(0,1), 1);
1353 auto node2 = mesh.add_point(Point(1,0), 2);
1354
1355 auto node3 = mesh.add_point(Point(2,0), 3);
1356 auto node4 = mesh.add_point(Point(2,1), 4);
1357 auto node5 = mesh.add_point(Point(3,0), 5);
1358
1359 auto edge01 = mesh.add_elem(Elem::build(EDGE2));
1360 edge01->set_node(0, node0);
1361 edge01->set_node(1, node1);
1362 auto edge12 = mesh.add_elem(Elem::build(EDGE2));
1363 edge12->set_node(0, node1);
1364 edge12->set_node(1, node2);
1365 auto edge20 = mesh.add_elem(Elem::build(EDGE2));
1366 edge20->set_node(0, node2);
1367 edge20->set_node(1, node0);
1368
1369 auto edge34 = mesh.add_elem(Elem::build(EDGE2));
1370 edge34->set_node(0, node3);
1371 edge34->set_node(1, node4);
1372 auto edge45 = mesh.add_elem(Elem::build(EDGE2));
1373 edge45->set_node(0, node4);
1374 edge45->set_node(1, node5);
1375 auto edge53 = mesh.add_elem(Elem::build(EDGE2));
1376 edge53->set_node(0, node5);
1377 edge53->set_node(1, node3);
1378
1379 mesh.prepare_for_use();
1380
1381 testExceptionBase(mesh, triangulator, "multiple loops of Edge");
1382 }
1383
1384 void testPoly2TriBad2DMultiBoundary()
1385 {
1386 LOG_UNIT_TEST;
1387
1388 Mesh mesh(*TestCommWorld);
1389 Poly2TriTriangulator triangulator(mesh);
1390
1391 // Two separate triangles
1392 auto node0 = mesh.add_point(Point(0,0), 0);
1393 auto node1 = mesh.add_point(Point(1,0), 1);
1394 auto node2 = mesh.add_point(Point(0,1), 2);
1395
1396 auto node3 = mesh.add_point(Point(2,0), 3);
1397 auto node4 = mesh.add_point(Point(3,0), 4);
1398 auto node5 = mesh.add_point(Point(2,1), 5);
1399
1400 auto tri012 = mesh.add_elem(Elem::build(TRI3));
1401 tri012->set_node(0, node0);
1402 tri012->set_node(1, node1);
1403 tri012->set_node(2, node2);
1404 auto tri345 = mesh.add_elem(Elem::build(TRI3));
1405 tri345->set_node(0, node3);
1406 tri345->set_node(1, node4);
1407 tri345->set_node(2, node5);
1408
1409 mesh.prepare_for_use();
1410
1411 testExceptionBase(mesh, triangulator, "cannot choose one");
1412 }
1413
1414
1415 void testPoly2TriEdgesRefined()
1416 {
1417 LOG_UNIT_TEST;
1418
1419 Mesh mesh(*TestCommWorld);
1420 Poly2TriTriangulator triangulator(mesh);
1421 testEdgesMesh(mesh, EDGE2);
1422 testPoly2TriRefinementBase(mesh, nullptr, 1.5, 14);
1423 }
1424
1425
1426 void testPoly2TriSegments()
1427 {
1428 LOG_UNIT_TEST;
1429
1430 Mesh mesh(*TestCommWorld);
1431 Poly2TriTriangulator p2t_tri(mesh);
1432 testTriangulatorSegments(mesh, p2t_tri);
1433 }
1434
1435 void testPoly2TriRefinementBase
1436 (UnstructuredMesh & mesh,
1437 const std::vector<TriangulatorInterface::Hole*> * holes,
1438 Real expected_total_area,
1439 dof_id_type n_original_elem,
1440 Real desired_area = 0.1,
1441 FunctionBase<Real> * area_func = nullptr)
1442 {
1443 Poly2TriTriangulator triangulator(mesh);
1444
1445 commonSettings(triangulator);
1446
1447 if (holes)
1448 triangulator.attach_hole_list(holes);
1449
1450 // Try to insert points!
1451 triangulator.desired_area() = desired_area;
1452 triangulator.set_desired_area_function(area_func);
1453
1454 triangulator.triangulate();
1455
1456 // If refinement should have increased our element count, check it
1457 if (desired_area || area_func)
1458 CPPUNIT_ASSERT_GREATER(n_original_elem, mesh.n_elem()); // n_elem+++
1459 else
1460 CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_original_elem);
1461
1462 Real area = 0;
1463 for (const auto & elem : mesh.active_local_element_ptr_range())
1464 {
1465 CPPUNIT_ASSERT_EQUAL(elem->level(), 0u);
1466 CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
1467
1468 const Real my_area = elem->volume();
1469
1470 // my_area <= desired_area, wow this macro ordering hurts
1471 if (desired_area != 0)
1472 CPPUNIT_ASSERT_LESSEQUAL(desired_area, my_area);
1473
1474 if (area_func != nullptr)
1475 for (auto v : make_range(elem->n_vertices()))
1476 {
1477 const Real local_desired_area =
1478 (*area_func)(elem->point(v));
1479 CPPUNIT_ASSERT_LESSEQUAL(local_desired_area, my_area);
1480 }
1481
1482 area += my_area;
1483 }
1484
1485 mesh.comm().sum(area);
1486
1487 LIBMESH_ASSERT_FP_EQUAL(area, expected_total_area, TOLERANCE*TOLERANCE);
1488 }
1489
1490 void testPoly2TriRefined()
1491 {
1492 LOG_UNIT_TEST;
1493
1494 Mesh mesh(*TestCommWorld);
1495 testTriangulatorTrapMesh(mesh);
1496 testPoly2TriRefinementBase(mesh, nullptr, 1.5, 15);
1497 }
1498
1499 void testPoly2TriNonRefined()
1500 {
1501 LOG_UNIT_TEST;
1502
1503 Mesh mesh(*TestCommWorld);
1504 testTriangulatorTrapMesh(mesh);
1505 // Make sure we see 0 as "don't refine", not "infinitely refine"
1506 testPoly2TriRefinementBase(mesh, nullptr, 1.5, 2, 0);
1507 }
1508
1509
1510 void testPoly2TriExtraRefined()
1511 {
1512 LOG_UNIT_TEST;
1513
1514 Mesh mesh(*TestCommWorld);
1515 testTriangulatorTrapMesh(mesh);
1516 testPoly2TriRefinementBase(mesh, nullptr, 1.5, 150, 0.01);
1517 }
1518
1519 void testPoly2TriNonUniformRefined()
1520 {
1521#ifdef LIBMESH_HAVE_FPARSER
1522 ParsedFunction<Real> var_area {"0.002*(1+2*x)*(1+2*y)"};
1523 Mesh mesh(*TestCommWorld);
1524 testTriangulatorTrapMesh(mesh);
1525 testPoly2TriRefinementBase(mesh, nullptr, 1.5, 150, 0, &var_area);
1526#endif // LIBMESH_HAVE_FPARSER
1527 }
1528
1529 void testPoly2TriHolesRefined()
1530 {
1531 LOG_UNIT_TEST;
1532
1533 // Add a diamond hole
1534 TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1535 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1536
1537 Mesh mesh(*TestCommWorld);
1538 testTriangulatorTrapMesh(mesh);
1539 testPoly2TriRefinementBase(mesh, &holes, 1.25, 13);
1540 }
1541
1542 void testPoly2TriHolesInterpRefined()
1543 {
1544 LOG_UNIT_TEST;
1545
1546 Mesh mesh(*TestCommWorld);
1547 testTriangulatorTrapMesh(mesh);
1548 Poly2TriTriangulator p2t_tri(mesh);
1549
1550 Real total_area = 1.5;
1551
1552 // Try a narrower point
1553 mesh.node_ref(2)(1) = 3;
1554 total_area += 0.5;
1555
1556 // Add a bunch of tiny diamond holes
1557 const int N=3, M=3;
1558 TriangulatorInterface::PolygonHole diamond(Point(), std::sqrt(2_R)/20, 4);
1559
1560 std::vector<TriangulatorInterface::AffineHole> hole_data;
1561 // Reserve so we don't invalidate pointers
1562 hole_data.reserve(M*N);
1563 std::vector<TriangulatorInterface::Hole*> holes;
1564 for (int i : make_range(M))
1565 for (int j : make_range(N))
1566 {
1567 Point shift(Real(i+1)/(M+1),Real(j+1)/(N+1));
1568 hole_data.emplace_back(diamond, 0, shift);
1569 holes.push_back(&hole_data.back());
1570 total_area -= hole_data.back().area();
1571 }
1572
1573 p2t_tri.attach_hole_list(&holes);
1574
1575 p2t_tri.set_refine_boundary_allowed(false);
1576
1577 testTriangulatorInterp(mesh, p2t_tri, 4, 0, total_area, 0.03_R);
1578 }
1579
1580 void testPoly2TriHolesInteriorRefinedBase
1581 (dof_id_type n_original_elem,
1582 Real desired_area)
1583 {
1584 // Add a diamond hole, disallowing refinement of it
1585 TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1586
1587 CPPUNIT_ASSERT_EQUAL(diamond.refine_boundary_allowed(), true);
1588
1589 diamond.set_refine_boundary_allowed(false);
1590
1591 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1592
1593 // Doing extra refinement here to ensure that we had the
1594 // *opportunity* to refine the hole boundaries.
1595 Mesh mesh(*TestCommWorld);
1596 testTriangulatorTrapMesh(mesh);
1597 testPoly2TriRefinementBase(mesh, &holes, 1.25, n_original_elem, desired_area);
1598
1599 // Checking that we have more outer boundary sides than we started
1600 // with, and exactly the 4 hole boundary sides we started with.
1601 auto side_bcs = mesh.get_boundary_info().build_side_list();
1602
1603 int n_outer_sides = std::count_if(side_bcs.begin(), side_bcs.end(),
1604 [](auto t){return std::get<2>(t) == 0;});
1605 CPPUNIT_ASSERT_GREATER(4, n_outer_sides); // n_outer_sides > 4
1606 int n_hole_sides = std::count_if(side_bcs.begin(), side_bcs.end(),
1607 [](auto t){return std::get<2>(t) == 1;});
1608 CPPUNIT_ASSERT_EQUAL(n_hole_sides, 4);
1609 }
1610
1611
1612 void testPoly2TriHolesInteriorRefined()
1613 {
1614 LOG_UNIT_TEST;
1615 testPoly2TriHolesInteriorRefinedBase(25, 0.05);
1616 }
1617
1618
1619 void testPoly2TriHolesInteriorExtraRefined()
1620 {
1621 LOG_UNIT_TEST;
1622 // 0.01 creates slivers triggering a poly2tri exception for me - RHS
1623 testPoly2TriHolesInteriorRefinedBase(60, 0.02);
1624 }
1625
1626
1627 void testPoly2TriHolesExtraRefined()
1628 {
1629 // Add a diamond hole
1630 TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1631 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1632
1633 Mesh mesh(*TestCommWorld);
1634 testTriangulatorTrapMesh(mesh);
1635 testPoly2TriRefinementBase(mesh, &holes, 1.25, 125, 0.01);
1636 }
1637
1638 void testPoly2TriHolesNonUniformRefined()
1639 {
1640#ifdef LIBMESH_HAVE_FPARSER
1641 // Add a diamond hole
1642 TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1643 const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1644
1645 ParsedFunction<Real> var_area {"0.002*(0.25+2*x)*(0.25+2*y)"};
1646 Mesh mesh(*TestCommWorld);
1647 testTriangulatorTrapMesh(mesh);
1648 testPoly2TriRefinementBase(mesh, &holes, 1.25, 150, 0, &var_area);
1649#endif // LIBMESH_HAVE_FPARSER
1650 }
1651
1652
1653#endif // LIBMESH_HAVE_POLY2TRI
1654
1655};
1656
1657
1658CPPUNIT_TEST_SUITE_REGISTRATION( MeshTriangulationTest );
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriNonUniformRefined)
MeshTriangulationTest::LIBMESH_CPPUNIT_TEST_SUITE
LIBMESH_CPPUNIT_TEST_SUITE(MeshTriangulationTest)
The goal of this test is to verify proper operation of the interfaces to triangulation libraries.
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleHalfDomain)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdge3s)
MeshTriangulationTest::testEdge3Mesh
void testEdge3Mesh(MeshBase &mesh)
Definition mesh_triangulation.C:849
MeshTriangulationTest::testTriangulatorInterp
void testTriangulatorInterp(UnstructuredMesh &mesh, TriangulatorInterface &triangulator, int interpolate_boundary_points, dof_id_type n_expected_elem, Real expected_total_area=1.5, Real desired_area=1000)
Definition mesh_triangulation.C:433
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleEdges)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolePerturbed)
MeshTriangulationTest::testEdge3ToTri6Base
void testEdge3ToTri6Base(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:371
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdge3ToTri6)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriInterp)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesInteriorRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriNonRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleInterp)
MeshTriangulationTest::testEdgesMesh
void testEdgesMesh(MeshBase &mesh, ElemType elem_type)
Definition mesh_triangulation.C:809
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdge3ToTri6FlipThrows)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesInteriorExtraRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriBadEdges)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriBad1DMultiBoundary)
MeshTriangulationTest::testTriangulatorHoles
void testTriangulatorHoles(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:511
MeshTriangulationTest::testTriangulatorBase
void testTriangulatorBase(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:317
MeshTriangulationTest::testFoundCenters
void testFoundCenters(const MeshBase &mesh, const std::vector< Point > &expected_centers)
Definition mesh_triangulation.C:471
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriBad2DMultiBoundary)
MeshTriangulationTest::testTriangulatorSegments
void testTriangulatorSegments(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:793
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleRoundHole)
MeshTriangulationTest::testTriangulatorRoundHole
void testTriangulatorRoundHole(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:714
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHoles)
MeshTriangulationTest::testExceptionBase
void testExceptionBase(MeshBase &mesh, TriangulatorInterface &triangulator, const char *re)
Definition mesh_triangulation.C:111
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleSegments)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesInterpRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleEdge3ToTri6)
MeshTriangulationTest::testPoly2TriRefinementBase
void testPoly2TriRefinementBase(UnstructuredMesh &mesh, const std::vector< TriangulatorInterface::Hole * > *holes, Real expected_total_area, dof_id_type n_original_elem, Real desired_area=0.1, FunctionBase< Real > *area_func=nullptr)
Definition mesh_triangulation.C:1436
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleMeshedHoles)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesExtraRefined)
MeshTriangulationTest::testHalfDomain
void testHalfDomain(MeshBase &mesh, TriangulatorInterface &triangulator, ElemType elem_type)
Definition mesh_triangulation.C:882
MeshTriangulationTest::testTriangulator
void testTriangulator(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:404
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleInterp2)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriRefined)
MeshTriangulationTest::testTriangulatorTrapMesh
void testTriangulatorTrapMesh(UnstructuredMesh &mesh)
Definition mesh_triangulation.C:420
MeshTriangulationTest::testTriangulatorHolePerturbed
void testTriangulatorHolePerturbed(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:547
MeshTriangulationTest::testPoly2TriHolesInteriorRefinedBase
void testPoly2TriHolesInteriorRefinedBase(dof_id_type n_original_elem, Real desired_area)
Definition mesh_triangulation.C:1581
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriSegments)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2Tri)
MeshTriangulationTest::testTriangulatorHoleTangentPerturbed
void testTriangulatorHoleTangentPerturbed(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:597
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleHoleContains)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangle)
MeshTriangulationTest::commonSettings
void commonSettings(TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:97
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdgesRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleHoleArea)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriMeshedHoles)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testHoleContainsRayDegeneracies)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHolesNonUniformRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHoleTangentPerturbed)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdge3ToTri7CenterFixup)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testTriangleHoles)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriEdges)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriExtraRefined)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHalfDomainEdge3)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriRoundHole)
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriInterp2)
MeshTriangulationTest::testTriangulatorMeshedHoles
void testTriangulatorMeshedHoles(MeshBase &mesh, TriangulatorInterface &triangulator)
Definition mesh_triangulation.C:663
MeshTriangulationTest::CPPUNIT_TEST
CPPUNIT_TEST(testPoly2TriHalfDomain)
libMesh::Parallel::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::BoundaryInfo::build_side_list
std::vector< BCTuple > build_side_list(BCTupleSortBy sort_by=BCTupleSortBy::ELEM_ID) const
Definition boundary_info.C:3048
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition elem.h:2567
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition elem.C:442
libMesh::FunctionBase
Base class for functors that can be evaluated at a point and (optionally) time.
Definition function_base.h:64
libMesh::LogicError
A class to represent the internal "this should never happen" errors, to be thrown by "libmesh_error()...
Definition libmesh_exceptions.h:66
libMesh::MeshBase
This is the MeshBase class.
Definition mesh_base.h:81
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition mesh_base.h:745
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition mesh_base.h:170
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition mesh_base.C:824
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::MeshBase::add_point
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
Add a new Node at Point p to the end of the vertex array, with processor_id procid.
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::Mesh
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition mesh.h:51
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition parallel_object.h:97
libMesh::ParsedFunction
A Function generated (via FParser) by parsing a mathematical expression.
Definition parsed_function.h:61
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::Poly2TriTriangulator
A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertio...
Definition poly2tri_triangulator.h:47
libMesh::Poly2TriTriangulator::set_refine_boundary_allowed
virtual void set_refine_boundary_allowed(bool refine_bdy_allowed) override
Set whether or not the triangulation is allowed to refine the mesh boundary when refining the interio...
Definition poly2tri_triangulator.h:97
libMesh::Poly2TriTriangulator::set_desired_area_function
virtual void set_desired_area_function(FunctionBase< Real > *desired) override
Set a function giving desired triangle area as a function of position.
Definition poly2tri_triangulator.C:451
libMesh::Poly2TriTriangulator::triangulate
virtual void triangulate() override
Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points ...
Definition poly2tri_triangulator.C:371
libMesh::TriangleInterface
A C++ interface between LibMesh and the Triangle library written by J.R.
Definition mesh_triangle_interface.h:42
libMesh::TriangulatorInterface::ArbitraryHole
Another concrete instantiation of the hole, this one should be sufficiently general for most non-poly...
Definition mesh_triangle_holes.h:265
libMesh::TriangulatorInterface::Hole
An abstract class for defining a 2-dimensional hole.
Definition mesh_triangle_holes.h:49
libMesh::TriangulatorInterface::Hole::refine_boundary_allowed
virtual bool refine_boundary_allowed() const
Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interio...
Definition mesh_triangle_holes.h:140
libMesh::TriangulatorInterface::Hole::area
Real area() const
Return the area of the hole.
Definition mesh_triangle_holes.C:190
libMesh::TriangulatorInterface::Hole::set_refine_boundary_allowed
virtual void set_refine_boundary_allowed(bool refine_bdy_allowed)
Set whether or not a triangulator is allowed to refine the hole boundary when refining the mesh inter...
Definition mesh_triangle_holes.h:133
libMesh::TriangulatorInterface::MeshedHole
Another concrete instantiation of the hole, as general as ArbitraryHole, but based on an existing 1D ...
Definition mesh_triangle_holes.h:344
libMesh::TriangulatorInterface::PolygonHole
A concrete instantiation of the Hole class that describes polygonal (triangular, square,...
Definition mesh_triangle_holes.h:172
libMesh::TriangulatorInterface::PolygonHole::point
virtual Point point(const unsigned int n) const override
Return the nth point defining the hole.
Definition mesh_triangle_holes.C:344
libMesh::TriangulatorInterface::PolygonHole::n_points
virtual unsigned int n_points() const override
The number of geometric points which define the hole.
Definition mesh_triangle_holes.C:338
libMesh::TriangulatorInterface::set_outer_boundary_ids
void set_outer_boundary_ids(std::set< std::size_t > bdy_ids)
A set of ids to allow on the outer boundary loop: interpreted as boundary ids of 2D elements and/or s...
Definition triangulator_interface.h:370
libMesh::TriangulatorInterface::minimum_angle
Real & minimum_angle()
Sets and/or gets the minimum desired angle.
Definition triangulator_interface.h:200
libMesh::TriangulatorInterface::triangulation_type
TriangulationType & triangulation_type()
Sets and/or gets the desired triangulation type.
Definition triangulator_interface.h:205
libMesh::TriangulatorInterface::triangulate
virtual void triangulate()=0
This is the main public interface for this function.
libMesh::TriangulatorInterface::PSLG
@ PSLG
Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of...
Definition triangulator_interface.h:130
libMesh::TriangulatorInterface::smooth_after_generating
bool & smooth_after_generating()
Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the gri...
Definition triangulator_interface.h:244
libMesh::TriangulatorInterface::set_verify_hole_boundaries
void set_verify_hole_boundaries(bool v)
Verifying that hole boundaries don't cross the outer boundary or each other is something like O(N_bdy...
Definition triangulator_interface.h:262
libMesh::TriangulatorInterface::segments
std::vector< std::pair< unsigned int, unsigned int > > segments
When constructing a PSLG, if the node numbers do not define the desired boundary segments implicitly ...
Definition triangulator_interface.h:290
libMesh::TriangulatorInterface::attach_hole_list
void attach_hole_list(const std::vector< Hole * > *holes)
Attaches a vector of Hole* pointers which will be meshed around.
Definition triangulator_interface.h:255
libMesh::TriangulatorInterface::desired_area
Real & desired_area()
Sets and/or gets the desired triangle area.
Definition triangulator_interface.h:177
libMesh::TriangulatorInterface::elem_type
ElemType & elem_type()
Sets and/or gets the desired element type.
Definition triangulator_interface.h:164
libMesh::TriangulatorInterface::set_fixup_tri7_center_nodes
void set_fixup_tri7_center_nodes(bool v)
For TRI7 elements, after boundary mid-edge nodes have been snapped to a curved boundary in increase_t...
Definition triangulator_interface.h:273
libMesh::TriangulatorInterface::set_interpolate_boundary_points
void set_interpolate_boundary_points(int n_points)
Complicated setter, for compatibility with insert_extra_points()
Definition triangulator_interface.C:129
libMesh::TypeVector::absolute_fuzzy_equals
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:975
libMesh::TypeVector::norm_sq
auto norm_sq() const
Definition type_vector.h:928
libMesh::TypeVector::relative_fuzzy_equals
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition type_vector.h:984
libMesh::TypeVector::norm
auto norm() const
Definition type_vector.h:908
libMesh::UnstructuredMesh
The UnstructuredMesh class is derived from the MeshBase class.
Definition unstructured_mesh.h:49
TestCommWorld
Communicator * TestCommWorld
libmesh_cppunit.h
mesh
MeshBase & mesh
Definition mesh_communication.C:1460
CPPUNIT_TEST_SUITE_REGISTRATION
CPPUNIT_TEST_SUITE_REGISTRATION(MeshTriangulationTest)
libMesh::MeshTools::Generation::build_square
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
Definition mesh_generation.C:1981
libMesh::MeshTools::Generation::build_sphere
void build_sphere(UnstructuredMesh &mesh, const Real radius=1, const unsigned int n_refinements=2, const ElemType type=INVALID_ELEM, const unsigned int n_smooth=2, const bool flat=true)
Fills mesh with a mesh discretizing a ball (||x||<= radius) or sphere (||x|| = radius) domain.
Definition mesh_generation.C:2012
libMesh::MeshTools::Modification::translate
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
Translates the mesh.
Definition mesh_modification.C:328
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::pi
const Real pi
.
Definition libmesh.h:292
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151
libMesh::dof_id_type
uint8_t dof_id_type
Definition id_types.h:67
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
radius
const Real radius
Definition subdomains_ex3.C:49
test_comm.h
Generated on Fri Aug 21 2026 13:32:31 for libMesh by doxygen 1.9.8