From 019725723e2a06e9ec5c71606480d9bcff27134b Mon Sep 17 00:00:00 2001
From: Yinbin Miao <ymiao@anl.gov>
Date: Wed, 21 Jun 2023 14:02:59 -0500
Subject: [PATCH 1/3] add random shifts to find_intersection()
closes #3496
---
src/mesh/mesh_triangle_holes.C | 5 ++++-
1 file changed, 4 insertions(+), 1 deletion(-)
diff --git a/src/mesh/mesh_triangle_holes.C b/src/mesh/mesh_triangle_holes.C
index 0250c1ea29d..5b45e593d56 100644
--- a/src/mesh/mesh_triangle_holes.C
+++ b/src/mesh/mesh_triangle_holes.C
@@ -105,11 +105,14 @@ namespace
//
// If the intersection is a "glancing" one at a corner, return -1.
Real find_intersection(const Point & source,
- const Point & ray_target,
+ const Point & ray_target_0,
const Point & edge_pt0,
const Point & edge_pt1,
const Point & edge_pt2)
{
+ // Add a random small shift to the source and three points
+ const Point ray_target = ray_target_0 + Point((Real)(rand() % 5000 + 5000) / 5000.0 * libMesh::TOLERANCE * libMesh::TOLERANCE, (Real)(rand() % 5000 + 5000) / 5000.0 * libMesh::TOLERANCE * libMesh::TOLERANCE, 0.0);
+
// Quick and more numerically stable check
if (!is_intersection(source, ray_target, edge_pt0, edge_pt1))
return -1;
From 57bab4ba8c145381ebf9eff588f8f0247f34711c Mon Sep 17 00:00:00 2001
From: Roy Stogner <Roy.Stogner@inl.gov>
Date: Mon, 20 Mar 2023 16:57:10 -0500
Subject: [PATCH 2/3] New unit test
This fails for me on master, and works for me with the new patch, so
hopefully it'll be a reliable way of making sure that any improved code
doesn't revert the fix.
---
tests/mesh/mesh_triangulation.C | 70 ++++++++++++++++++++++++++++++++-
1 file changed, 69 insertions(+), 1 deletion(-)
diff --git a/tests/mesh/mesh_triangulation.C b/tests/mesh/mesh_triangulation.C
index f28e817a4ab..440aefe51f9 100644
--- a/tests/mesh/mesh_triangulation.C
+++ b/tests/mesh/mesh_triangulation.C
@@ -54,6 +54,9 @@ public:
// This covers an old poly2tri collinearity-tolerance bug
CPPUNIT_TEST( testPoly2TriHolesExtraRefined );
+ // This covers a more recent tolerance bug when verifying holes
+ CPPUNIT_TEST( testPoly2TriHolePerturbed );
+
CPPUNIT_TEST( testPoly2TriNonUniformRefined );
CPPUNIT_TEST( testPoly2TriHolesNonUniformRefined );
#endif
@@ -64,6 +67,7 @@ public:
CPPUNIT_TEST( testTriangleInterp );
CPPUNIT_TEST( testTriangleInterp2 );
CPPUNIT_TEST( testTriangleHoles );
+ CPPUNIT_TEST( testTriangleHolePerturbed );
CPPUNIT_TEST( testTriangleMeshedHoles );
CPPUNIT_TEST( testTriangleEdges );
CPPUNIT_TEST( testTriangleSegments );
@@ -84,7 +88,7 @@ public:
triangulator.triangulation_type() = TriangulatorInterface::PSLG;
// Don't try to insert points unless we're requested to later
- triangulator.desired_area() = 1000;
+ triangulator.desired_area() = 1e16;
triangulator.minimum_angle() = 0;
triangulator.smooth_after_generating() = false;
triangulator.set_verify_hole_boundaries(true);
@@ -469,6 +473,50 @@ public:
}
+ void testTriangulatorHolePerturbed(MeshBase & mesh,
+ TriangulatorInterface & triangulator)
+ {
+ // Points based on a simplification of a hole verification failure
+ // case
+ mesh.add_point(Point(100,0), 0);
+ mesh.add_point(Point(100,100), 1);
+ mesh.add_point(Point(0,100), 2);
+ mesh.add_point(Point(-100,100), 3);
+ mesh.add_point(Point(-100,0), 4);
+ mesh.add_point(Point(-100,-100), 5);
+ mesh.add_point(Point(0,-100), 6);
+ mesh.add_point(Point(100,-100), 7);
+
+ commonSettings(triangulator);
+
+ // Add a diamond hole in *almost* the center
+ TriangulatorInterface::PolygonHole
+ diamond(Point(0,4.e-16),
+ std::sqrt(2)/2, 4);
+ const std::vector<TriangulatorInterface::Hole*> holes { &diamond };
+ triangulator.attach_hole_list(&holes);
+
+ triangulator.triangulate();
+
+ CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(12));
+
+ // Center coordinates for all the elements we expect
+ const Real r2p200o6 = (std::sqrt(Real(2))+200)/6,
+ r2p400o6 = (std::sqrt(Real(2))+400)/6;
+
+ std::vector <Point> expected_centers
+ { {r2p400o6,100./3}, {-r2p400o6,100./3},
+ {r2p400o6,-100./3}, {-r2p400o6,-100./3},
+ {100./3,r2p400o6}, {-100./3,r2p400o6},
+ {100./3,-r2p400o6}, {-100./3,-r2p400o6},
+ {r2p200o6,r2p200o6}, {-r2p200o6,r2p200o6},
+ {r2p200o6,-r2p200o6}, {-r2p200o6,-r2p200o6},
+ };
+
+ testFoundCenters(mesh, expected_centers);
+ }
+
+
void testTriangulatorMeshedHoles(MeshBase & mesh,
TriangulatorInterface & triangulator)
{
@@ -652,6 +700,16 @@ public:
}
+ void testTriangleHolePerturbed()
+ {
+ LOG_UNIT_TEST;
+
+ Mesh mesh(*TestCommWorld);
+ TriangleInterface triangle(mesh);
+ testTriangulatorHolePerturbed(mesh, triangle);
+ }
+
+
void testTriangleMeshedHoles()
{
LOG_UNIT_TEST;
@@ -736,6 +794,16 @@ public:
}
+ void testPoly2TriHolePerturbed()
+ {
+ LOG_UNIT_TEST;
+
+ Mesh mesh(*TestCommWorld);
+ Poly2TriTriangulator p2t_tri(mesh);
+ testTriangulatorHolePerturbed(mesh, p2t_tri);
+ }
+
+
void testPoly2TriMeshedHoles()
{
LOG_UNIT_TEST;
From 256691480b92c34383107b420c1e77fab2cc4df8 Mon Sep 17 00:00:00 2001
From: Roy Stogner <Roy.Stogner@inl.gov>
Date: Tue, 21 Mar 2023 17:20:08 -0500
Subject: [PATCH 3/3] Unit tests for "tangent" ray/hole intersection
---
tests/mesh/mesh_triangulation.C | 83 ++++++++++++++++++++++++++++++++-
1 file changed, 82 insertions(+), 1 deletion(-)
diff --git a/tests/mesh/mesh_triangulation.C b/tests/mesh/mesh_triangulation.C
index 440aefe51f9..86e3a54fa6a 100644
--- a/tests/mesh/mesh_triangulation.C
+++ b/tests/mesh/mesh_triangulation.C
@@ -54,8 +54,9 @@ public:
// This covers an old poly2tri collinearity-tolerance bug
CPPUNIT_TEST( testPoly2TriHolesExtraRefined );
- // This covers a more recent tolerance bug when verifying holes
+ // These cover more recent tolerance issues when verifying holes
CPPUNIT_TEST( testPoly2TriHolePerturbed );
+ CPPUNIT_TEST( testPoly2TriHoleTangentPerturbed );
CPPUNIT_TEST( testPoly2TriNonUniformRefined );
CPPUNIT_TEST( testPoly2TriHolesNonUniformRefined );
@@ -68,6 +69,7 @@ public:
CPPUNIT_TEST( testTriangleInterp2 );
CPPUNIT_TEST( testTriangleHoles );
CPPUNIT_TEST( testTriangleHolePerturbed );
+ CPPUNIT_TEST( testTriangleHoleTangentPerturbed );
CPPUNIT_TEST( testTriangleMeshedHoles );
CPPUNIT_TEST( testTriangleEdges );
CPPUNIT_TEST( testTriangleSegments );
@@ -517,6 +519,65 @@ public:
}
+ void testTriangulatorHoleTangentPerturbed(MeshBase & mesh,
+ TriangulatorInterface & triangulator)
+ {
+ // Points based on a simplification of a hole verification failure
+ // case
+ mesh.add_point(Point(200,0), 0);
+ mesh.add_point(Point(200,100), 1);
+ mesh.add_point(Point(100,100), 2);
+ mesh.add_point(Point(0,100), 3);
+ mesh.add_point(Point(-100,100), 4);
+ mesh.add_point(Point(-200,100), 5);
+ mesh.add_point(Point(-200,0), 6);
+ mesh.add_point(Point(-200,-100), 7);
+ mesh.add_point(Point(-100,-100), 8);
+ mesh.add_point(Point(0,-100), 9);
+ mesh.add_point(Point(100,-100), 10);
+ mesh.add_point(Point(200,-100), 11);
+
+ commonSettings(triangulator);
+
+ // Two diamond holes, in *almost* the center of each half of that
+ // rectangle
+ TriangulatorInterface::PolygonHole
+ left_diamond(Point(-100,4.e-16),
+ std::sqrt(2)/2, 4),
+ right_diamond(Point(100,-4.e-16),
+ std::sqrt(2)/2, 4);
+ const std::vector<TriangulatorInterface::Hole*> holes
+ { &left_diamond, &right_diamond };
+ triangulator.attach_hole_list(&holes);
+
+ triangulator.triangulate();
+
+ CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), dof_id_type(22));
+
+ // Center coordinates for all the elements we expect
+ const Real r2p200o6 = (std::sqrt(Real(2))+200)/6,
+ r2p400o6 = (std::sqrt(Real(2))+400)/6;
+
+ std::vector <Point> expected_centers
+ { {50+r2p400o6,100./3},
+ {50+r2p400o6,-100./3},
+ {50+100./3,r2p400o6}, {50-100./3,r2p400o6},
+ {50+100./3,-r2p400o6}, {50-100./3,-r2p400o6},
+ {50+r2p200o6,r2p200o6}, {50-r2p200o6,r2p200o6},
+ {50+r2p200o6,-r2p200o6}, {50-r2p200o6,-r2p200o6},
+ {50-r2p400o6,100./3},
+ {50-r2p400o6,-100./3},
+ {-50+100./3,r2p400o6}, {50-100./3,r2p400o6},
+ {-50+100./3,-r2p400o6}, {50-100./3,-r2p400o6},
+ {-50+r2p200o6,r2p200o6}, {50-r2p200o6,r2p200o6},
+ {-50+r2p200o6,-r2p200o6}, {50-r2p200o6,-r2p200o6},
+ {0,100./3}, {0,-100./3}
+ };
+
+ testFoundCenters(mesh, expected_centers);
+ }
+
+
void testTriangulatorMeshedHoles(MeshBase & mesh,
TriangulatorInterface & triangulator)
{
@@ -710,6 +771,16 @@ public:
}
+ void testTriangleHoleTangentPerturbed()
+ {
+ LOG_UNIT_TEST;
+
+ Mesh mesh(*TestCommWorld);
+ TriangleInterface triangle(mesh);
+ testTriangulatorHoleTangentPerturbed(mesh, triangle);
+ }
+
+
void testTriangleMeshedHoles()
{
LOG_UNIT_TEST;
@@ -804,6 +875,16 @@ public:
}
+ void testPoly2TriHoleTangentPerturbed()
+ {
+ LOG_UNIT_TEST;
+
+ Mesh mesh(*TestCommWorld);
+ Poly2TriTriangulator p2t_tri(mesh);
+ testTriangulatorHoleTangentPerturbed(mesh, p2t_tri);
+ }
+
+
void testPoly2TriMeshedHoles()
{
LOG_UNIT_TEST;