libMesh
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 
21 using namespace libMesh;
22 
23 class MeshTriangulationTest : public CppUnit::TestCase
24 {
29 public:
30  LIBMESH_CPPUNIT_TEST_SUITE( MeshTriangulationTest );
31 
32  CPPUNIT_TEST( testTriangleHoleArea );
33  CPPUNIT_TEST( testTriangleHoleContains );
34 
35 #ifdef LIBMESH_HAVE_POLY2TRI
36  CPPUNIT_TEST( testPoly2Tri );
37  CPPUNIT_TEST( testPoly2TriHalfDomain );
38  CPPUNIT_TEST( testPoly2TriHalfDomainEdge3 );
39  CPPUNIT_TEST( testPoly2TriInterp );
40  CPPUNIT_TEST( testPoly2TriInterp2 );
41  CPPUNIT_TEST( testPoly2TriHoles );
42  CPPUNIT_TEST( testPoly2TriMeshedHoles );
43  CPPUNIT_TEST( testPoly2TriEdge3ToTri6 );
44 # ifdef LIBMESH_ENABLE_AMR
45  CPPUNIT_TEST( testPoly2TriRoundHole );
46 # endif
47  CPPUNIT_TEST( testPoly2TriEdges );
48  CPPUNIT_TEST( testPoly2TriEdge3s );
49  CPPUNIT_TEST( testPoly2TriBadEdges );
50  CPPUNIT_TEST( testPoly2TriBad1DMultiBoundary );
51  CPPUNIT_TEST( testPoly2TriBad2DMultiBoundary );
52  CPPUNIT_TEST( testPoly2TriEdgesRefined );
53  CPPUNIT_TEST( testPoly2TriSegments );
54  CPPUNIT_TEST( testPoly2TriRefined );
55  CPPUNIT_TEST( testPoly2TriNonRefined );
56  CPPUNIT_TEST( testPoly2TriExtraRefined );
57  CPPUNIT_TEST( testPoly2TriHolesRefined );
58  CPPUNIT_TEST( testPoly2TriHolesInterpRefined );
59  CPPUNIT_TEST( testPoly2TriHolesInteriorRefined );
60  CPPUNIT_TEST( testPoly2TriHolesInteriorExtraRefined );
61  // This covers an old poly2tri collinearity-tolerance bug
62  CPPUNIT_TEST( testPoly2TriHolesExtraRefined );
63 
64  CPPUNIT_TEST( testPoly2TriNonUniformRefined );
65  CPPUNIT_TEST( testPoly2TriHolesNonUniformRefined );
66 #endif
67 
68 #ifdef LIBMESH_HAVE_TRIANGLE
69  CPPUNIT_TEST( testTriangle );
70  CPPUNIT_TEST( testTriangleHalfDomain );
71  CPPUNIT_TEST( testTriangleInterp );
72  CPPUNIT_TEST( testTriangleInterp2 );
73  CPPUNIT_TEST( testTriangleHoles );
74  CPPUNIT_TEST( testTriangleMeshedHoles );
75 # ifdef LIBMESH_ENABLE_AMR
76  CPPUNIT_TEST( testTriangleRoundHole );
77 # endif
78  CPPUNIT_TEST( testTriangleEdges );
79  CPPUNIT_TEST( testTriangleSegments );
80  CPPUNIT_TEST( testTriangleEdge3ToTri6 );
81 #endif
82 
83  CPPUNIT_TEST_SUITE_END();
84 
85 public:
86  void setUp() {}
87 
88  void tearDown() {}
89 
90  void commonSettings (TriangulatorInterface & triangulator)
91  {
92  // Use the point order to define the boundary, because our
93  // Poly2Tri implementation doesn't do convex hulls yet, even when
94  // that would give the same answer.
96 
97  // Don't try to insert points unless we're requested to later
98  triangulator.desired_area() = 1000;
99  triangulator.minimum_angle() = 0;
100  triangulator.smooth_after_generating() = false;
101  triangulator.set_verify_hole_boundaries(true);
102  }
103 
105  TriangulatorInterface & triangulator,
106  const char * re)
107  {
108 #ifdef LIBMESH_ENABLE_EXCEPTIONS
109  // We can't just CPPUNIT_ASSERT_THROW, because we want to make
110  // sure we were thrown from the right place with the right error
111  // message!
112  bool threw_desired_exception = false;
113  try {
114  this->testTriangulatorBase(mesh, triangulator);
115  }
116  catch (libMesh::LogicError & e) {
117  std::regex msg_regex(re);
118  CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
119  threw_desired_exception = true;
120  }
121  catch (CppUnit::Exception & e) {
122  throw e;
123  }
124  catch (...) {
125  CPPUNIT_ASSERT_MESSAGE("Unexpected exception type thrown", false);
126  }
127  CPPUNIT_ASSERT(threw_desired_exception);
128 #endif
129  }
130 
131 
133  {
134  LOG_UNIT_TEST;
135 
136  // Using center=(1,0), radius=2 for the heck of it
137  Point center{1};
138  Real radius = 2;
139  std::vector<TriangulatorInterface::PolygonHole> polyholes;
140 
141  // Line
142  polyholes.emplace_back(center, radius, 2);
143  // Triangle
144  polyholes.emplace_back(center, radius, 3);
145  // Square
146  polyholes.emplace_back(center, radius, 4);
147  // Pentagon
148  polyholes.emplace_back(center, radius, 5);
149  // Hexagon
150  polyholes.emplace_back(center, radius, 6);
151 
152  for (int i=0; i != 5; ++i)
153  {
154  const int n_sides = i+2;
155  const TriangulatorInterface::Hole & hole = polyholes[i];
156 
157  const Real computed_area = hole.area();
158  const Real theta = pi/n_sides;
159  const Real half_side_length = radius*std::cos(theta);
160  const Real apothem = radius*std::sin(theta);
161  const Real area = n_sides * apothem * half_side_length;
162 
163  LIBMESH_ASSERT_FP_EQUAL(computed_area, area, TOLERANCE*TOLERANCE);
164  }
165 
166  TriangulatorInterface::ArbitraryHole arbhole {center, {{0,-1},{2,-1},{2,1},{0,2}}};
167  LIBMESH_ASSERT_FP_EQUAL(arbhole.area(), Real(5), TOLERANCE*TOLERANCE);
168 
169 #ifdef LIBMESH_HAVE_TRIANGLE
170  // Make sure we're compatible with the old naming structure too
171  TriangleInterface::PolygonHole square(center, radius, 4);
172  LIBMESH_ASSERT_FP_EQUAL(square.area(), 2*radius*radius, TOLERANCE*TOLERANCE);
173 #endif
174  }
175 
176 
178  {
179  LOG_UNIT_TEST;
180 
181  // Using center=(1,2), radius=2 for the heck of it
182  Point center{1,2};
183  Real radius = 2;
184  std::vector<TriangulatorInterface::PolygonHole> polyholes;
185 
186  auto check_corners = [center, radius]
187  (const TriangulatorInterface::Hole & hole,
188  unsigned int np)
189  {
190  CPPUNIT_ASSERT_EQUAL(np, hole.n_points());
191  for (auto i : make_range(np))
192  {
193  const Real theta = i * 2 * libMesh::pi / np;
194  const Real xin = center(0) + radius * .99 * std::cos(theta);
195  const Real xout = center(0) + radius * 1.01 * std::cos(theta);
196  const Real yin = center(1) + radius * .99 * std::sin(theta);
197  const Real yout = center(1) + radius * 1.01 * std::sin(theta);
198 
199  CPPUNIT_ASSERT(hole.contains(Point(xin, yin)));
200  CPPUNIT_ASSERT(!hole.contains(Point(xout, yout)));
201  }
202  };
203 
204  const TriangulatorInterface::PolygonHole triangle(center, radius, 3);
205  check_corners(triangle, 3);
206  const TriangulatorInterface::PolygonHole diamond(center, radius, 4);
207  check_corners(diamond, 4);
208  const TriangulatorInterface::PolygonHole hexagon(center, radius, 6);
209  check_corners(hexagon, 6);
210 
212  {{1,0}, {{0,-1},{2,-1},{2,1},{1.75,-.5},{1.5,1},{1.25,-.5},
213  {1,1},{.75,-.5},{.5,1},{.25,-.5},{0,1}}};
214 
215  CPPUNIT_ASSERT(jaggy.contains({.1,-.3}));
216  CPPUNIT_ASSERT(jaggy.contains({.5,.9}));
217  CPPUNIT_ASSERT(jaggy.contains({.9,-.3}));
218  CPPUNIT_ASSERT(jaggy.contains({1,0}));
219  CPPUNIT_ASSERT(jaggy.contains({1.1,-.4}));
220  CPPUNIT_ASSERT(jaggy.contains({1.5,.9}));
221  CPPUNIT_ASSERT(jaggy.contains({1.9,-.3}));
222 
223  CPPUNIT_ASSERT(jaggy.contains({1.9,-.5}));
224  CPPUNIT_ASSERT(jaggy.contains({1.6,-.5}));
225  CPPUNIT_ASSERT(jaggy.contains({1.1,-.5}));
226  CPPUNIT_ASSERT(jaggy.contains({.5,-.5}));
227  CPPUNIT_ASSERT(jaggy.contains({.2,-.5}));
228 
229  CPPUNIT_ASSERT(!jaggy.contains({.1,.7}));
230  CPPUNIT_ASSERT(!jaggy.contains({.5,1.1}));
231  CPPUNIT_ASSERT(!jaggy.contains({.9,.7}));
232  CPPUNIT_ASSERT(!jaggy.contains({1,-1.1}));
233  CPPUNIT_ASSERT(!jaggy.contains({1.1,.8}));
234  CPPUNIT_ASSERT(!jaggy.contains({1.5,1.1}));
235  CPPUNIT_ASSERT(!jaggy.contains({1.9,.9}));
236 
237  CPPUNIT_ASSERT(!jaggy.contains({1.9,1}));
238  CPPUNIT_ASSERT(!jaggy.contains({1.4,1}));
239  CPPUNIT_ASSERT(!jaggy.contains({.9,1}));
240  CPPUNIT_ASSERT(!jaggy.contains({.4,1}));
241  CPPUNIT_ASSERT(!jaggy.contains({-.2,1}));
242  CPPUNIT_ASSERT(!jaggy.contains({-.2,0}));
243  CPPUNIT_ASSERT(!jaggy.contains({1.2,0}));
244 
246  {{1,0}, {{-.25,-1},{2,-1},{2,1},{1.75,1},{1.75,-.5},{1.5,-.5},
247  {1.5,1},{1.25,1},{1.25,-.5},{1,-.5},{1,1},{.75,1},
248  {.75,-.5},{.5,-.5},{.5,1},{.25,1},{.25,-.5},{0,-.5},
249  {0,1},{-.25,1}}};
250 
251  CPPUNIT_ASSERT(square_jaggy.contains({-.1,-.3}));
252  CPPUNIT_ASSERT(square_jaggy.contains({.4,.9}));
253  CPPUNIT_ASSERT(square_jaggy.contains({.9,-.3}));
254  CPPUNIT_ASSERT(square_jaggy.contains({.9,0}));
255  CPPUNIT_ASSERT(square_jaggy.contains({1.1,-.6}));
256  CPPUNIT_ASSERT(square_jaggy.contains({1.4,.9}));
257  CPPUNIT_ASSERT(square_jaggy.contains({1.9,-.3}));
258 
259  CPPUNIT_ASSERT(!square_jaggy.contains({.1,-.3}));
260  CPPUNIT_ASSERT(!square_jaggy.contains({.6,.9}));
261  CPPUNIT_ASSERT(!square_jaggy.contains({1.1,-.3}));
262  CPPUNIT_ASSERT(!square_jaggy.contains({1.1,0}));
263  CPPUNIT_ASSERT(!square_jaggy.contains({1.1,-1.6}));
264  CPPUNIT_ASSERT(!square_jaggy.contains({1.6,.9}));
265  CPPUNIT_ASSERT(!square_jaggy.contains({2.1,-.3}));
266 
267  CPPUNIT_ASSERT(square_jaggy.contains({-.1,-.5}));
268  CPPUNIT_ASSERT(square_jaggy.contains({.3,-.5}));
269  CPPUNIT_ASSERT(square_jaggy.contains({.9,-.5}));
270  CPPUNIT_ASSERT(square_jaggy.contains({1.3,-.5}));
271  CPPUNIT_ASSERT(square_jaggy.contains({1.9,-.5}));
272 
273  CPPUNIT_ASSERT(!square_jaggy.contains({-.3,1}));
274  CPPUNIT_ASSERT(!square_jaggy.contains({.2,1}));
275  CPPUNIT_ASSERT(!square_jaggy.contains({.6,1}));
276  CPPUNIT_ASSERT(!square_jaggy.contains({1.1,1}));
277  CPPUNIT_ASSERT(!square_jaggy.contains({1.6,1}));
278  CPPUNIT_ASSERT(!square_jaggy.contains({2.1,1}));
279  }
280 
281 
283  TriangulatorInterface & triangulator)
284  {
285  commonSettings(triangulator);
286 
287  triangulator.triangulate();
288 
289  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(2));
290  for (const auto & elem : mesh.element_ptr_range())
291  {
292  CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
293 
294  // Make sure we're not getting any inverted elements
295  auto cross_prod =
296  (elem->point(1) - elem->point(0)).cross
297  (elem->point(2) - elem->point(0));
298 
299  CPPUNIT_ASSERT_GREATER(Real(0), cross_prod(2));
300 
301  bool found_triangle = false;
302  for (const auto & node : elem->node_ref_range())
303  {
304  const Point & point = node;
305  if (point == Point(0,0))
306  {
307  found_triangle = true;
308  CPPUNIT_ASSERT((elem->point(0) == Point(0,0) &&
309  elem->point(1) == Point(1,0) &&
310  elem->point(2) == Point(0,1)) ||
311  (elem->point(1) == Point(0,0) &&
312  elem->point(2) == Point(1,0) &&
313  elem->point(0) == Point(0,1)) ||
314  (elem->point(2) == Point(0,0) &&
315  elem->point(0) == Point(1,0) &&
316  elem->point(1) == Point(0,1)));
317  }
318  if (point == Point(1,2))
319  {
320  found_triangle = true;
321  CPPUNIT_ASSERT((elem->point(0) == Point(0,1) &&
322  elem->point(1) == Point(1,0) &&
323  elem->point(2) == Point(1,2)) ||
324  (elem->point(1) == Point(0,1) &&
325  elem->point(2) == Point(1,0) &&
326  elem->point(0) == Point(1,2)) ||
327  (elem->point(2) == Point(0,1) &&
328  elem->point(0) == Point(1,0) &&
329  elem->point(1) == Point(1,2)));
330  }
331  }
332  CPPUNIT_ASSERT(found_triangle);
333  }
334  }
335 
337  TriangulatorInterface & triangulator)
338  {
339  commonSettings(triangulator);
340 
341  triangulator.elem_type() = TRI6;
342  triangulator.triangulate();
343  const Real hsq2 = std::sqrt(2) / 2.0;
344 
345  // Check element number
346  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), (dof_id_type)2);
347  for (const auto & elem : mesh.element_ptr_range())
348  {
349  // Check element type
350  CPPUNIT_ASSERT_EQUAL(elem->type(), TRI6);
351 
352  for (const auto & i_side : make_range(elem->n_sides()))
353  {
354  // We only want to check the sides on the external boundary
355  if (elem->neighbor_ptr(i_side) == nullptr)
356  {
357  const Point & side_pt_0 = *(elem->node_ptr(i_side));
358  const Point & side_pt_1 = *(elem->node_ptr((i_side + 1) % 3));
359  const Point & side_pt_2 = *(elem->node_ptr(i_side + 3));
360  const bool x_sign = (side_pt_0(0) == 1 || side_pt_1(0) == 1);
361  const bool y_sign = (side_pt_0(1) == 1 || side_pt_1(1) == 1);
362  const Point ref_side_pt_2 = Point (x_sign ? hsq2 : -hsq2, y_sign ? hsq2 : -hsq2);
363  CPPUNIT_ASSERT_EQUAL(ref_side_pt_2, side_pt_2);
364  }
365  }
366  }
367  }
368 
370  TriangulatorInterface & triangulator)
371  {
372  // A non-square quad, so we don't have ambiguity about which
373  // diagonal a Delaunay algorithm will pick.
374  // Manually-numbered points, so we can use the point numbering as
375  // a segment ordering even on DistributedMesh.
376  mesh.add_point(Point(0,0), 0);
377  mesh.add_point(Point(1,0), 1);
378  mesh.add_point(Point(1,2), 2);
379  mesh.add_point(Point(0,1), 3);
380 
381  this->testTriangulatorBase(mesh, triangulator);
382  }
383 
384 
386  {
387  // A non-square quad, so we don't have ambiguity about which
388  // diagonal a Delaunay algorithm will pick.
389  // Manually-numbered points, so we can use the point numbering as
390  // a segment ordering even on DistributedMesh.
391  mesh.add_point(Point(0,0), 0);
392  mesh.add_point(Point(1,0), 1);
393  mesh.add_point(Point(1,2), 2);
394  mesh.add_point(Point(0,1), 3);
395  }
396 
397 
399  TriangulatorInterface & triangulator,
400  int interpolate_boundary_points,
401  dof_id_type n_expected_elem,
402  Real expected_total_area = 1.5,
403  Real desired_area = 1000)
404  {
405  commonSettings(triangulator);
406 
407  if (!mesh.n_nodes())
408  testTriangulatorTrapMesh(mesh);
409 
410  // Interpolate points!
411  triangulator.set_interpolate_boundary_points(interpolate_boundary_points);
412 
413  // Try to insert points?
414  triangulator.desired_area() = desired_area;
415 
416  triangulator.triangulate();
417 
418  if (n_expected_elem)
419  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_expected_elem);
420 
421  Real area = 0;
422  for (const auto & elem : mesh.active_local_element_ptr_range())
423  {
424  CPPUNIT_ASSERT_EQUAL(elem->level(), 0u);
425  CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
426 
427  area += elem->volume();
428  }
429 
430  mesh.comm().sum(area);
431 
432  LIBMESH_ASSERT_FP_EQUAL(area, expected_total_area, TOLERANCE*TOLERANCE);
433  }
434 
435 
437  const std::vector<Point> & expected_centers)
438  {
439  std::vector<bool> found_centers(expected_centers.size(), false);
440 
441  for (const auto & elem : mesh.element_ptr_range())
442  {
443  CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
444 
445  // Make sure we're not getting any inverted elements
446  auto cross_prod =
447  (elem->point(1) - elem->point(0)).cross
448  (elem->point(2) - elem->point(0));
449 
450  CPPUNIT_ASSERT_GREATER(Real(0), cross_prod(2));
451 
452  // Make sure we're finding all the elements we expect
453  Point center = elem->vertex_average();
454 
455  bool found_mine = false;
456  for (auto i : index_range(expected_centers))
457  {
458  Point possible = expected_centers[i];
459 
460  if (possible.absolute_fuzzy_equals(center, TOLERANCE*TOLERANCE))
461  {
462  found_mine = true;
463  found_centers[i] = true;
464  }
465  }
466  CPPUNIT_ASSERT(found_mine);
467  }
468 
469  mesh.comm().max(found_centers);
470 
471  for (auto found_it : found_centers)
472  CPPUNIT_ASSERT(found_it);
473  }
474 
475 
477  TriangulatorInterface & triangulator)
478  {
479  // A square quad; we'll put a diamond hole in the middle to make
480  // the Delaunay selection unambiguous.
481  mesh.add_point(Point(-1,-1), 0);
482  mesh.add_point(Point(1,-1), 1);
483  mesh.add_point(Point(1,1), 2);
484  mesh.add_point(Point(-1,1), 3);
485 
486  commonSettings(triangulator);
487 
488  // Add a diamond hole in the center
489  TriangulatorInterface::PolygonHole diamond(Point(0), std::sqrt(2_R)/2, 4);
490  const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
491  triangulator.attach_hole_list(&holes);
492 
493  triangulator.triangulate();
494 
495  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(8));
496 
497  // Center coordinates for all the elements we expect
498  Real r2p2o6 = (std::sqrt(2_R)+2)/6;
499  Real r2p4o6 = (std::sqrt(2_R)+4)/6;
500 
501  std::vector <Point> expected_centers
502  { {r2p2o6,r2p2o6}, {r2p2o6,-r2p2o6},
503  {-r2p2o6,r2p2o6}, {-r2p2o6,-r2p2o6},
504  {0,r2p4o6}, {r2p4o6, 0},
505  {0,-r2p4o6}, {-r2p4o6, 0}
506  };
507 
508  testFoundCenters(mesh, expected_centers);
509  }
510 
511 
513  TriangulatorInterface & triangulator)
514  {
515  // A square quad; we'll put a square hole in the middle. Offset
516  // this to catch a potential bug I missed on the first try.
517  mesh.add_point(Point(19,19), 0);
518  mesh.add_point(Point(21,19), 1);
519  mesh.add_point(Point(21,21), 2);
520  mesh.add_point(Point(19,21), 3);
521 
522  commonSettings(triangulator);
523 
524  // Add a square meshed hole in the center
525  Mesh centermesh { mesh.comm() };
526  MeshTools::Generation::build_square (centermesh, 2, 2, 19.5, 20.5, 19.5, 20.5, QUAD4);
527 
528  TriangulatorInterface::MeshedHole centerhole { centermesh };
529 
530  CPPUNIT_ASSERT_EQUAL(centerhole.n_points(), 8u);
531  CPPUNIT_ASSERT_EQUAL(centerhole.area(), Real(1));
532  Point inside = centerhole.inside();
533  CPPUNIT_ASSERT_EQUAL(inside(0), Real(20));
534  CPPUNIT_ASSERT_EQUAL(inside(1), Real(20));
535 
536  const std::vector<TriangulatorInterface::Hole*> holes { &centerhole };
537  triangulator.attach_hole_list(&holes);
538 
539  triangulator.triangulate();
540 
541  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(12));
542 
543  // Center coordinates for all the elements we expect
544  std::vector <Point> expected_centers
545  { {-0.5, Real(2)/3}, {0, Real(5)/6},
546  {0.5, Real(2)/3},
547  {Real(2)/3, -0.5}, {Real(5)/6, 0},
548  {Real(2)/3, 0.5},
549  {-0.5, -Real(2)/3}, {0, -Real(5)/6},
550  {0.5, -Real(2)/3},
551  {-Real(2)/3, -0.5}, {-Real(5)/6, 0},
552  {-Real(2)/3, 0.5} };
553 
554  // With the same offset
555  for (auto & p : expected_centers)
556  p += Point(20,20);
557 
558  testFoundCenters(mesh, expected_centers);
559  }
560 
561 
562 #ifdef LIBMESH_ENABLE_AMR
564  TriangulatorInterface & triangulator)
565  {
566  // A square quad; we'll put a round hole in the middle.
567  mesh.add_point(Point(19,19), 0);
568  mesh.add_point(Point(21,19), 1);
569  mesh.add_point(Point(21,21), 2);
570  mesh.add_point(Point(19,21), 3);
571 
572  commonSettings(triangulator);
573 
574  // Add a square meshed hole in the center
575  const Real radius = 0.5;
576  const Point center{20,20};
577 
578  Mesh centermesh { mesh.comm() };
580  MeshTools::Modification::translate(centermesh, center(0), center(1));
581 
582  TriangulatorInterface::MeshedHole centerhole { centermesh };
583 
584  CPPUNIT_ASSERT_EQUAL(centerhole.n_points(), 8u);
585  Point inside = centerhole.inside();
586  CPPUNIT_ASSERT_EQUAL(inside(0), Real(20));
587  CPPUNIT_ASSERT_EQUAL(inside(1), Real(20));
588 
589  const std::vector<TriangulatorInterface::Hole*> holes { &centerhole };
590  triangulator.attach_hole_list(&holes);
591  triangulator.elem_type() = TRI6;
592 
593  triangulator.triangulate();
594 
595  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(12));
596  CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), static_cast<dof_id_type>(36));
597 
598  // Make sure we didn't screw up the outer sides. We should
599  // have exact values for the outer vertices, so we can use
600  // those for a map.
601  std::map<std::pair<Real, Real>, Point> outer_midpoints
602  {{{19, 19}, {20, 19}},
603  {{21, 19}, {21, 20}},
604  {{21, 21}, {20, 21}},
605  {{19, 21}, {19, 20}}
606  };
607 
608  for (const auto & elem : mesh.active_local_element_ptr_range())
609  for (const auto n : make_range(elem->n_sides()))
610  {
611  if (elem->neighbor_ptr(n))
612  continue;
613 
614  auto it =
615  outer_midpoints.find(std::make_pair(elem->point(n)(0),
616  elem->point(n)(1)));
617  if (it != outer_midpoints.end())
618  {
619  const Point error = it->second - elem->point(n+3);
620  CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
621  error.norm_sq());
622  }
623  else
624  {
625  const Point radius1 = elem->point(n) - center;
626  CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
627  std::abs(radius1.norm()-radius));
628 
629  const Point radius2 = elem->point((n+1)%3) - center;
630  CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
631  std::abs(radius2.norm()-radius));
632 
633  const Point radius3 = elem->point(n+3) - center;
634  CPPUNIT_ASSERT_LESS(TOLERANCE*TOLERANCE,
635  std::abs(radius3.norm()-radius));
636  }
637  }
638  }
639 #endif // LIBMESH_ENABLE_AMR
640 
641 
643  TriangulatorInterface & triangulator)
644  {
645  // The same quad as testTriangulator, but out of order
646  mesh.add_point(Point(0,0), 0);
647  mesh.add_point(Point(1,2), 1);
648  mesh.add_point(Point(1,0), 2);
649  mesh.add_point(Point(0,1), 3);
650 
651  // Segments to put them in order
652  triangulator.segments = {{0,2},{2,1},{1,3},{3,0}};
653 
654  this->testTriangulatorBase(mesh, triangulator);
655  }
656 
657 
658  void testEdgesMesh(MeshBase & mesh, ElemType elem_type)
659  {
660  // The same quad as testTriangulator, but out of order
661  auto node0 = mesh.add_point(Point(0,0), 0);
662  auto node1 = mesh.add_point(Point(1,2), 1);
663  auto node2 = mesh.add_point(Point(1,0), 2);
664  auto node3 = mesh.add_point(Point(0,1), 3);
665 
666  // Edges, also out of order, but enough to put them in order
667  auto edge13 = mesh.add_elem(Elem::build(elem_type));
668  edge13->set_node(0, node1);
669  edge13->set_node(1, node3);
670  auto edge02 = mesh.add_elem(Elem::build(elem_type));
671  edge02->set_node(0, node0);
672  edge02->set_node(1, node2);
673  auto edge30 = mesh.add_elem(Elem::build(elem_type));
674  edge30->set_node(0, node3);
675  edge30->set_node(1, node0);
676  auto edge21 = mesh.add_elem(Elem::build(elem_type));
677  edge21->set_node(0, node2);
678  edge21->set_node(1, node1);
679 
680  // Add mid-edge nodes if asked to
681  if (elem_type == EDGE3)
682  {
683  auto node4 = mesh.add_point(Point(.5,1.5), 4);
684  edge13->set_node(2, node4);
685  auto node5 = mesh.add_point(Point(.5,0), 5);
686  edge02->set_node(2, node5);
687  auto node6 = mesh.add_point(Point(0,.5), 6);
688  edge30->set_node(2, node6);
689  auto node7 = mesh.add_point(Point(1,1), 7);
690  edge21->set_node(2, node7);
691  }
692  else
693  libmesh_assert(elem_type == EDGE2);
694 
696  }
697 
699  {
700  const Real hsq2 = std::sqrt(2) / 2.0;
701  auto node0 = mesh.add_point(Point(1,0), 0);
702  auto node1 = mesh.add_point(Point(hsq2,hsq2), 1);
703  auto node2 = mesh.add_point(Point(0,1), 2);
704  auto node3 = mesh.add_point(Point(-hsq2,hsq2), 3);
705  auto node4 = mesh.add_point(Point(-1,0), 4);
706  auto node5 = mesh.add_point(Point(-hsq2,-hsq2), 5);
707  auto node6 = mesh.add_point(Point(0,-1), 6);
708  auto node7 = mesh.add_point(Point(hsq2,-hsq2), 7);
709 
710  auto edge021 = mesh.add_elem(Elem::build(EDGE3));
711  edge021->set_node(0, node0);
712  edge021->set_node(1, node2);
713  edge021->set_node(2, node1);
714  auto edge243 = mesh.add_elem(Elem::build(EDGE3));
715  edge243->set_node(0, node2);
716  edge243->set_node(1, node4);
717  edge243->set_node(2, node3);
718  auto edge465 = mesh.add_elem(Elem::build(EDGE3));
719  edge465->set_node(0, node4);
720  edge465->set_node(1, node6);
721  edge465->set_node(2, node5);
722  auto edge607 = mesh.add_elem(Elem::build(EDGE3));
723  edge607->set_node(0, node6);
724  edge607->set_node(1, node0);
725  edge607->set_node(2, node7);
726 
728  }
729 
730 
732  TriangulatorInterface & triangulator,
733  ElemType elem_type)
734  {
735  // We might have 4 or 8 nodes on the outer boundary, depending on
736  // whether it has mid-edge nodes
737  const dof_id_type off = (elem_type == EDGE3)*4;
738 
739  // A pentagon we'll avoid via subdomain ids
740  auto node4 = mesh.add_point(Point(2,0), 4+off);
741  auto node5 = mesh.add_point(Point(3,0), 5+off);
742  auto node6 = mesh.add_point(Point(3,2), 6+off);
743  auto node7 = mesh.add_point(Point(2,2), 7+off);
744  auto node8 = mesh.add_point(Point(2,1), 8+off);
745 
746  auto edge45 = mesh.add_elem(Elem::build(elem_type));
747  edge45->set_node(0, node4);
748  edge45->set_node(1, node5);
749  edge45->subdomain_id() = 1;
750  auto edge56 = mesh.add_elem(Elem::build(elem_type));
751  edge56->set_node(0, node5);
752  edge56->set_node(1, node6);
753  edge56->subdomain_id() = 1;
754  auto edge67 = mesh.add_elem(Elem::build(elem_type));
755  edge67->set_node(0, node6);
756  edge67->set_node(1, node7);
757  edge67->subdomain_id() = 1;
758  auto edge78 = mesh.add_elem(Elem::build(elem_type));
759  edge78->set_node(0, node7);
760  edge78->set_node(1, node8);
761  edge78->subdomain_id() = 1;
762  auto edge84 = mesh.add_elem(Elem::build(elem_type));
763  edge84->set_node(0, node8);
764  edge84->set_node(1, node4);
765  edge84->subdomain_id() = 1;
766 
767  if (elem_type == EDGE3)
768  {
769  auto node9 = mesh.add_point(Point(2.5,0), 9+off);
770  edge45->set_node(2, node9);
771  auto node10 = mesh.add_point(Point(3,1), 10+off);
772  edge56->set_node(2, node10);
773  auto node11 = mesh.add_point(Point(2.5,2), 11+off);
774  edge67->set_node(2, node11);
775  auto node12 = mesh.add_point(Point(2,1.5), 12+off);
776  edge78->set_node(2, node12);
777  auto node13 = mesh.add_point(Point(2,.5), 13+off);
778  edge84->set_node(2, node13);
779  }
780  else
781  libmesh_assert(elem_type == EDGE2);
782 
783  testEdgesMesh(mesh, elem_type);
784 
785  std::set<std::size_t> bdy_ids {0};
786  triangulator.set_outer_boundary_ids(bdy_ids);
787 
788  this->testTriangulatorBase(mesh, triangulator);
789  }
790 
791 
792 #ifdef LIBMESH_HAVE_TRIANGLE
794  {
795  LOG_UNIT_TEST;
796 
798  TriangleInterface triangle(mesh);
799  testTriangulator(mesh, triangle);
800  }
801 
802 
804  {
805  LOG_UNIT_TEST;
806 
808  TriangleInterface triangle(mesh);
809  testHalfDomain(mesh, triangle, EDGE2);
810  }
811 
812 
814  {
815  LOG_UNIT_TEST;
816 
818  TriangleInterface triangle(mesh);
819  testTriangulatorInterp(mesh, triangle, 1, 6);
820  }
821 
822 
824  {
825  LOG_UNIT_TEST;
826 
828  TriangleInterface triangle(mesh);
829  testTriangulatorInterp(mesh, triangle, 2, 10);
830  }
831 
832 
834  {
835  LOG_UNIT_TEST;
836 
838  TriangleInterface triangle(mesh);
839  testTriangulatorHoles(mesh, triangle);
840  }
841 
842 
844  {
845  LOG_UNIT_TEST;
846 
848  TriangleInterface triangle(mesh);
849  testTriangulatorMeshedHoles(mesh, triangle);
850  }
851 
852 
853 #ifdef LIBMESH_ENABLE_AMR
855  {
856  LOG_UNIT_TEST;
857 
859  TriangleInterface triangle(mesh);
860  testTriangulatorRoundHole(mesh, triangle);
861  }
862 #endif // LIBMESH_ENABLE_AMR
863 
864 
866  {
867  LOG_UNIT_TEST;
868 
870  TriangleInterface triangulator(mesh);
871  testEdgesMesh(mesh, EDGE2);
872 
873  this->testTriangulatorBase(mesh, triangulator);
874  }
875 
877  {
878  LOG_UNIT_TEST;
879 
881  TriangleInterface triangle(mesh);
882  testTriangulatorSegments(mesh, triangle);
883  }
884 
886  {
887  LOG_UNIT_TEST;
888 
890  TriangleInterface triangulator(mesh);
891 
892  testEdge3Mesh(mesh);
893 
894  this->testEdge3ToTri6Base(mesh, triangulator);
895  }
896 
897 #endif // LIBMESH_HAVE_TRIANGLE
898 
899 
900 #ifdef LIBMESH_HAVE_POLY2TRI
902  {
903  LOG_UNIT_TEST;
904 
906  Poly2TriTriangulator p2t_tri(mesh);
907  testTriangulator(mesh, p2t_tri);
908  }
909 
910 
912  {
913  LOG_UNIT_TEST;
914 
916  Poly2TriTriangulator p2t_tri(mesh);
917  testHalfDomain(mesh, p2t_tri, EDGE2);
918  }
919 
920 
922  {
923  LOG_UNIT_TEST;
924 
926  Poly2TriTriangulator p2t_tri(mesh);
927  testHalfDomain(mesh, p2t_tri, EDGE3);
928  }
929 
930 
932  {
933  LOG_UNIT_TEST;
934 
936  Poly2TriTriangulator p2t_tri(mesh);
937  testTriangulatorInterp(mesh, p2t_tri, 1, 6);
938  }
939 
940 
942  {
943  LOG_UNIT_TEST;
944 
946  Poly2TriTriangulator p2t_tri(mesh);
947  testTriangulatorInterp(mesh, p2t_tri, 2, 10);
948  }
949 
950 
952  {
953  LOG_UNIT_TEST;
954 
956  Poly2TriTriangulator p2t_tri(mesh);
957  testTriangulatorHoles(mesh, p2t_tri);
958  }
959 
960 
962  {
963  LOG_UNIT_TEST;
964 
966  Poly2TriTriangulator p2t_tri(mesh);
967  testTriangulatorMeshedHoles(mesh, p2t_tri);
968  }
969 
971  {
972  LOG_UNIT_TEST;
973 
975  Poly2TriTriangulator triangulator(mesh);
976 
977  testEdge3Mesh(mesh);
978 
979  this->testEdge3ToTri6Base(mesh, triangulator);
980  }
981 
982 
983 #ifdef LIBMESH_ENABLE_AMR
985  {
986  LOG_UNIT_TEST;
987 
989  Poly2TriTriangulator p2t_tri(mesh);
990  testTriangulatorRoundHole(mesh, p2t_tri);
991  }
992 #endif // LIBMESH_ENABLE_AMR
993 
994 
996  {
997  LOG_UNIT_TEST;
998 
1000  Poly2TriTriangulator triangulator(mesh);
1001  testEdgesMesh(mesh, EDGE2);
1002 
1003  this->testTriangulatorBase(mesh, triangulator);
1004  }
1005 
1006 
1008  {
1009  LOG_UNIT_TEST;
1010 
1012  Poly2TriTriangulator triangulator(mesh);
1013  testEdgesMesh(mesh, EDGE3);
1014 
1015  this->testTriangulatorBase(mesh, triangulator);
1016  }
1017 
1018 
1020  {
1021  LOG_UNIT_TEST;
1022 
1024  Poly2TriTriangulator triangulator(mesh);
1025 
1026  // The same quad as testTriangulator, but out of order
1027  auto node0 = mesh.add_point(Point(0,0), 0);
1028  auto node1 = mesh.add_point(Point(1,2), 1);
1029  auto node2 = mesh.add_point(Point(1,0), 2);
1030  auto node3 = mesh.add_point(Point(0,1), 3);
1031 
1032  // Edges, but not enough to complete the quad
1033  auto edge13 = mesh.add_elem(Elem::build(EDGE2));
1034  edge13->set_node(0, node1);
1035  edge13->set_node(1, node3);
1036  auto edge02 = mesh.add_elem(Elem::build(EDGE2));
1037  edge02->set_node(0, node0);
1038  edge02->set_node(1, node2);
1039  auto edge30 = mesh.add_elem(Elem::build(EDGE2));
1040  edge30->set_node(0, node3);
1041  edge30->set_node(1, node0);
1042 
1044 
1045  testExceptionBase(mesh, triangulator, "Bad edge topology");
1046  }
1047 
1048 
1050  {
1051  LOG_UNIT_TEST;
1052 
1054  Poly2TriTriangulator triangulator(mesh);
1055 
1056  // Two separate triangles
1057  auto node0 = mesh.add_point(Point(0,0), 0);
1058  auto node1 = mesh.add_point(Point(0,1), 1);
1059  auto node2 = mesh.add_point(Point(1,0), 2);
1060 
1061  auto node3 = mesh.add_point(Point(2,0), 3);
1062  auto node4 = mesh.add_point(Point(2,1), 4);
1063  auto node5 = mesh.add_point(Point(3,0), 5);
1064 
1065  auto edge01 = mesh.add_elem(Elem::build(EDGE2));
1066  edge01->set_node(0, node0);
1067  edge01->set_node(1, node1);
1068  auto edge12 = mesh.add_elem(Elem::build(EDGE2));
1069  edge12->set_node(0, node1);
1070  edge12->set_node(1, node2);
1071  auto edge20 = mesh.add_elem(Elem::build(EDGE2));
1072  edge20->set_node(0, node2);
1073  edge20->set_node(1, node0);
1074 
1075  auto edge34 = mesh.add_elem(Elem::build(EDGE2));
1076  edge34->set_node(0, node3);
1077  edge34->set_node(1, node4);
1078  auto edge45 = mesh.add_elem(Elem::build(EDGE2));
1079  edge45->set_node(0, node4);
1080  edge45->set_node(1, node5);
1081  auto edge53 = mesh.add_elem(Elem::build(EDGE2));
1082  edge53->set_node(0, node5);
1083  edge53->set_node(1, node3);
1084 
1086 
1087  testExceptionBase(mesh, triangulator, "multiple loops of Edge");
1088  }
1089 
1091  {
1092  LOG_UNIT_TEST;
1093 
1095  Poly2TriTriangulator triangulator(mesh);
1096 
1097  // Two separate triangles
1098  auto node0 = mesh.add_point(Point(0,0), 0);
1099  auto node1 = mesh.add_point(Point(1,0), 1);
1100  auto node2 = mesh.add_point(Point(0,1), 2);
1101 
1102  auto node3 = mesh.add_point(Point(2,0), 3);
1103  auto node4 = mesh.add_point(Point(3,0), 4);
1104  auto node5 = mesh.add_point(Point(2,1), 5);
1105 
1106  auto tri012 = mesh.add_elem(Elem::build(TRI3));
1107  tri012->set_node(0, node0);
1108  tri012->set_node(1, node1);
1109  tri012->set_node(2, node2);
1110  auto tri345 = mesh.add_elem(Elem::build(TRI3));
1111  tri345->set_node(0, node3);
1112  tri345->set_node(1, node4);
1113  tri345->set_node(2, node5);
1114 
1116 
1117  testExceptionBase(mesh, triangulator, "cannot choose one");
1118  }
1119 
1120 
1122  {
1123  LOG_UNIT_TEST;
1124 
1126  Poly2TriTriangulator triangulator(mesh);
1127  testEdgesMesh(mesh, EDGE2);
1128  testPoly2TriRefinementBase(mesh, nullptr, 1.5, 14);
1129  }
1130 
1131 
1133  {
1134  LOG_UNIT_TEST;
1135 
1137  Poly2TriTriangulator p2t_tri(mesh);
1138  testTriangulatorSegments(mesh, p2t_tri);
1139  }
1140 
1141  void testPoly2TriRefinementBase
1143  const std::vector<TriangulatorInterface::Hole*> * holes,
1144  Real expected_total_area,
1145  dof_id_type n_original_elem,
1146  Real desired_area = 0.1,
1147  FunctionBase<Real> * area_func = nullptr)
1148  {
1149  Poly2TriTriangulator triangulator(mesh);
1150 
1151  commonSettings(triangulator);
1152 
1153  if (holes)
1154  triangulator.attach_hole_list(holes);
1155 
1156  // Try to insert points!
1157  triangulator.desired_area() = desired_area;
1158  triangulator.set_desired_area_function(area_func);
1159 
1160  triangulator.triangulate();
1161 
1162  // If refinement should have increased our element count, check it
1163  if (desired_area || area_func)
1164  CPPUNIT_ASSERT_GREATER(n_original_elem, mesh.n_elem()); // n_elem+++
1165  else
1166  CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_original_elem);
1167 
1168  Real area = 0;
1169  for (const auto & elem : mesh.active_local_element_ptr_range())
1170  {
1171  CPPUNIT_ASSERT_EQUAL(elem->level(), 0u);
1172  CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
1173 
1174  const Real my_area = elem->volume();
1175 
1176  // my_area <= desired_area, wow this macro ordering hurts
1177  if (desired_area != 0)
1178  CPPUNIT_ASSERT_LESSEQUAL(desired_area, my_area);
1179 
1180  if (area_func != nullptr)
1181  for (auto v : make_range(elem->n_vertices()))
1182  {
1183  const Real local_desired_area =
1184  (*area_func)(elem->point(v));
1185  CPPUNIT_ASSERT_LESSEQUAL(local_desired_area, my_area);
1186  }
1187 
1188  area += my_area;
1189  }
1190 
1191  mesh.comm().sum(area);
1192 
1193  LIBMESH_ASSERT_FP_EQUAL(area, expected_total_area, TOLERANCE*TOLERANCE);
1194  }
1195 
1197  {
1198  LOG_UNIT_TEST;
1199 
1201  testTriangulatorTrapMesh(mesh);
1202  testPoly2TriRefinementBase(mesh, nullptr, 1.5, 15);
1203  }
1204 
1206  {
1207  LOG_UNIT_TEST;
1208 
1210  testTriangulatorTrapMesh(mesh);
1211  // Make sure we see 0 as "don't refine", not "infinitely refine"
1212  testPoly2TriRefinementBase(mesh, nullptr, 1.5, 2, 0);
1213  }
1214 
1215 
1217  {
1218  LOG_UNIT_TEST;
1219 
1221  testTriangulatorTrapMesh(mesh);
1222  testPoly2TriRefinementBase(mesh, nullptr, 1.5, 150, 0.01);
1223  }
1224 
1226  {
1227 #ifdef LIBMESH_HAVE_FPARSER
1228  ParsedFunction<Real> var_area {"0.002*(1+2*x)*(1+2*y)"};
1230  testTriangulatorTrapMesh(mesh);
1231  testPoly2TriRefinementBase(mesh, nullptr, 1.5, 150, 0, &var_area);
1232 #endif // LIBMESH_HAVE_FPARSER
1233  }
1234 
1236  {
1237  LOG_UNIT_TEST;
1238 
1239  // Add a diamond hole
1240  TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1241  const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1242 
1244  testTriangulatorTrapMesh(mesh);
1245  testPoly2TriRefinementBase(mesh, &holes, 1.25, 13);
1246  }
1247 
1249  {
1250  LOG_UNIT_TEST;
1251 
1253  testTriangulatorTrapMesh(mesh);
1254  Poly2TriTriangulator p2t_tri(mesh);
1255 
1256  Real total_area = 1.5;
1257 
1258  // Try a narrower point
1259  mesh.node_ref(2)(1) = 3;
1260  total_area += 0.5;
1261 
1262  // Add a bunch of tiny diamond holes
1263  const int N=3, M=3;
1264  TriangulatorInterface::PolygonHole diamond(Point(), std::sqrt(2_R)/20, 4);
1265 
1266  std::vector<TriangulatorInterface::AffineHole> hole_data;
1267  // Reserve so we don't invalidate pointers
1268  hole_data.reserve(M*N);
1269  std::vector<TriangulatorInterface::Hole*> holes;
1270  for (int i : make_range(M))
1271  for (int j : make_range(N))
1272  {
1273  Point shift(Real(i+1)/(M+1),Real(j+1)/(N+1));
1274  hole_data.emplace_back(diamond, 0, shift);
1275  holes.push_back(&hole_data.back());
1276  total_area -= hole_data.back().area();
1277  }
1278 
1279  p2t_tri.attach_hole_list(&holes);
1280 
1281  p2t_tri.set_refine_boundary_allowed(false);
1282 
1283  testTriangulatorInterp(mesh, p2t_tri, 4, 0, total_area, 0.03_R);
1284  }
1285 
1286  void testPoly2TriHolesInteriorRefinedBase
1287  (dof_id_type n_original_elem,
1288  Real desired_area)
1289  {
1290  // Add a diamond hole, disallowing refinement of it
1291  TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1292 
1293  CPPUNIT_ASSERT_EQUAL(diamond.refine_boundary_allowed(), true);
1294 
1295  diamond.set_refine_boundary_allowed(false);
1296 
1297  const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1298 
1299  // Doing extra refinement here to ensure that we had the
1300  // *opportunity* to refine the hole boundaries.
1302  testTriangulatorTrapMesh(mesh);
1303  testPoly2TriRefinementBase(mesh, &holes, 1.25, n_original_elem, desired_area);
1304 
1305  // Checking that we have more outer boundary sides than we started
1306  // with, and exactly the 4 hole boundary sides we started with.
1307  auto side_bcs = mesh.get_boundary_info().build_side_list();
1308 
1309  int n_outer_sides = std::count_if(side_bcs.begin(), side_bcs.end(),
1310  [](auto t){return std::get<2>(t) == 0;});
1311  CPPUNIT_ASSERT_GREATER(4, n_outer_sides); // n_outer_sides > 4
1312  int n_hole_sides = std::count_if(side_bcs.begin(), side_bcs.end(),
1313  [](auto t){return std::get<2>(t) == 1;});
1314  CPPUNIT_ASSERT_EQUAL(n_hole_sides, 4);
1315  }
1316 
1317 
1319  {
1320  LOG_UNIT_TEST;
1321  testPoly2TriHolesInteriorRefinedBase(25, 0.05);
1322  }
1323 
1324 
1326  {
1327  LOG_UNIT_TEST;
1328  // 0.01 creates slivers triggering a poly2tri exception for me - RHS
1329  testPoly2TriHolesInteriorRefinedBase(60, 0.02);
1330  }
1331 
1332 
1334  {
1335  // Add a diamond hole
1336  TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1337  const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1338 
1340  testTriangulatorTrapMesh(mesh);
1341  testPoly2TriRefinementBase(mesh, &holes, 1.25, 125, 0.01);
1342  }
1343 
1345  {
1346 #ifdef LIBMESH_HAVE_FPARSER
1347  // Add a diamond hole
1348  TriangulatorInterface::PolygonHole diamond(Point(0.5,0.5), std::sqrt(2_R)/4, 4);
1349  const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
1350 
1351  ParsedFunction<Real> var_area {"0.002*(0.25+2*x)*(0.25+2*y)"};
1353  testTriangulatorTrapMesh(mesh);
1354  testPoly2TriRefinementBase(mesh, &holes, 1.25, 150, 0, &var_area);
1355 #endif // LIBMESH_HAVE_FPARSER
1356  }
1357 
1358 
1359 #endif // LIBMESH_HAVE_POLY2TRI
1360 
1361 };
1362 
1363 
void testTriangulatorTrapMesh(UnstructuredMesh &mesh)
ElemType
Defines an enum for geometric element types.
void build_sphere(UnstructuredMesh &mesh, const Real rad=1, const unsigned int nr=2, const ElemType type=INVALID_ELEM, const unsigned int n_smooth=2, const bool flat=true)
Meshes a spherical or mapped-spherical domain.
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
const Real radius
virtual void triangulate()=0
This is the main public interface for this function.
auto norm() const -> decltype(std::norm(T()))
Definition: type_vector.h:907
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...
A Function generated (via FParser) by parsing a mathematical expression.
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
static constexpr Real TOLERANCE
An abstract class for defining a 2-dimensional hole.
void set_verify_hole_boundaries(bool v)
Verifying that hole boundaries don&#39;t cross the outer boundary or each other is something like O(N_bdy...
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly ecreated (or read) mesh for use.
Definition: mesh_base.C:759
void sum(T &r) const
void testTriangulatorSegments(MeshBase &mesh, TriangulatorInterface &triangulator)
Real area() const
Return the area of the hole.
MeshBase & mesh
const Parallel::Communicator & comm() const
void set_interpolate_boundary_points(int n_points)
Complicated setter, for compatibility with insert_extra_points()
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.
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...
The libMesh namespace provides an interface to certain functionality in the library.
ElemType & elem_type()
Sets and/or gets the desired element type.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:165
void testTriangulatorHoles(MeshBase &mesh, TriangulatorInterface &triangulator)
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.
void testHalfDomain(MeshBase &mesh, TriangulatorInterface &triangulator, ElemType elem_type)
This is the MeshBase class.
Definition: mesh_base.h:75
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
Creates a list of element numbers, sides, and ids for those sides.
TriangulationType & triangulation_type()
Sets and/or gets the desired triangulation type.
void testEdge3ToTri6Base(MeshBase &mesh, TriangulatorInterface &triangulator)
void commonSettings(TriangulatorInterface &triangulator)
auto norm_sq() const -> decltype(std::norm(T()))
Definition: type_vector.h:926
void attach_hole_list(const std::vector< Hole *> *holes)
Attaches a vector of Hole* pointers which will be meshed around.
virtual void set_desired_area_function(FunctionBase< Real > *desired) override
Set a function giving desired triangle area as a function of position.
The UnstructuredMesh class is derived from the MeshBase class.
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:444
Real & minimum_angle()
Sets and/or gets the minimum desired angle.
libmesh_assert(ctx)
CPPUNIT_TEST_SUITE_REGISTRATION(MeshTriangulationTest)
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)
bool absolute_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:972
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
Translates the mesh.
virtual bool refine_boundary_allowed() const
Get whether or not the triangulation is allowed to refine the mesh boundary when refining the interio...
virtual void triangulate() override
Internally, this calls the poly2tri triangulation code in a loop, inserting our owner Steiner points ...
Real & desired_area()
Sets and/or gets the desired triangle area.
Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of...
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
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 ...
void testTriangulatorRoundHole(MeshBase &mesh, TriangulatorInterface &triangulator)
void max(const T &r, T &o, Request &req) const
A class to represent the internal "this should never happen" errors, to be thrown by "libmesh_error()...
Another concrete instantiation of the hole, as general as ArbitraryHole, but based on an existing 1D ...
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...
A concrete instantiation of the Hole class that describes polygonal (triangular, square, pentagonal, ...) holes.
void testTriangulatorBase(MeshBase &mesh, TriangulatorInterface &triangulator)
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:140
void testTriangulator(MeshBase &mesh, TriangulatorInterface &triangulator)
void testEdgesMesh(MeshBase &mesh, ElemType elem_type)
Another concrete instantiation of the hole, this one should be sufficiently general for most non-poly...
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
A C++ interface between LibMesh and the poly2tri library, with custom code for Steiner point insertio...
void testEdge3Mesh(MeshBase &mesh)
virtual dof_id_type n_elem() const =0
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50
A C++ interface between LibMesh and the Triangle library written by J.R.
void testTriangulatorMeshedHoles(MeshBase &mesh, TriangulatorInterface &triangulator)
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
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:117
virtual dof_id_type n_nodes() const =0
void testExceptionBase(MeshBase &mesh, TriangulatorInterface &triangulator, const char *re)
bool & smooth_after_generating()
Sets/gets flag which tells whether to do two steps of Laplace mesh smoothing after generating the gri...
void testFoundCenters(const MeshBase &mesh, const std::vector< Point > &expected_centers)
uint8_t dof_id_type
Definition: id_types.h:67
const Real pi
.
Definition: libmesh.h:299