More Poly2TriTriangulator robustness fixes

Our interpretation of number of interpolation points didn't match our
documentation or intentions, so some unit tests here had to be changed
along with the code.

Other than that, this just refactors us to be more robust to use cases
like incoming meshes with weird node numberings.

Poly2Tri is now passing all the tests I can throw at this in both MOOSE
and libMesh.

Author: roystgnr

include/mesh/triangulator_interface.h

@@ -234,11 +234,17 @@ class TriangulatorInterface
 protected:
   /**
    * Helper function to create PSLG segments from our other
-   * boundary-defining options (nodal ordering, 1D mesh edges, 2D mesh
+   * boundary-defining options (1D mesh edges, 2D mesh
    * boundary sides), if no segments already exist.
    */
   void elems_to_segments();
 
+  /**
+   * Helper function to create PSLG segments from our node ordering,
+   * up to the maximum node id, if no segments already exist.
+   */
+  void nodes_to_segments(dof_id_type max_node_id);
+
   /**
    * Helper function to add extra points (midpoints of initial
    * segments) to a PSLG triangulation

src/mesh/mesh_triangle_interface.C

@@ -60,12 +60,16 @@ void TriangleInterface::triangulate()
   // Will the triangulation have holes?
   const bool have_holes = ((_holes != nullptr) && (!_holes->empty()));
 
+  unsigned int n_hole_points = this->total_hole_points();
+
+  // If we're doing PSLG without segments, construct them from all our
+  // mesh nodes
+  this->nodes_to_segments(_mesh.max_node_id());
+
   // Insert additional new points in between existing boundary points,
   // if that is requested and reasonable
   this->insert_any_extra_boundary_points();
 
-  unsigned int n_hole_points = this->total_hole_points();
-
   // Regardless of whether we added additional points, the set of points to
   // triangulate is now sitting in the mesh.
 
@@ -339,6 +343,8 @@ void TriangleInterface::triangulate()
   }
 
 
+  _mesh.set_mesh_dimension(2);
+
 
   // To the naked eye, a few smoothing iterations usually looks better,
   // so we do this by default unless the user says not to.

src/mesh/poly2tri_triangulator.C

@@ -189,6 +189,14 @@ void Poly2TriTriangulator::triangulate()
       _elem_type != TRI7)
     libmesh_not_implemented();
 
+  // If we have no explicit segments defined, we may get them from
+  // mesh elements
+  this->elems_to_segments();
+
+  // If we *still* have no explicit segments defined, we get them from
+  // the order of nodes.
+  this->nodes_to_segments(_n_boundary_nodes);
+
   // Insert additional new points in between existing boundary points,
   // if that is requested and reasonable
   this->insert_any_extra_boundary_points();
@@ -204,6 +212,8 @@ void Poly2TriTriangulator::triangulate()
   if (_markers)
     libmesh_not_implemented();
 
+  _mesh.set_mesh_dimension(2);
+
   // To the naked eye, a few smoothing iterations usually looks better,
   // so we do this by default unless the user says not to.
   if (this->_smooth_after_generating)
@@ -307,10 +317,6 @@ void Poly2TriTriangulator::triangulate_current_points()
   // Prepare poly2tri points for our nodes, sorted into outer boundary
   // points and interior Steiner points.
 
-  // If we have no explicit segments defined, we may get them from
-  // mesh elements
-  this->elems_to_segments();
-
   if (this->segments.empty())
     {
       // If we have no segments even after taking elems into account,
@@ -586,6 +592,7 @@ bool Poly2TriTriangulator::insert_refinement_points()
   // In cases where we've been working with contiguous node id ranges;
   // let's keep it that way.
   dof_id_type nn = _mesh.max_node_id();
+  dof_id_type ne = _mesh.max_elem_id();
 
   for (auto & elem : mesh.element_ptr_range())
     {
@@ -855,7 +862,7 @@ bool Poly2TriTriangulator::insert_refinement_points()
                       continue;
                     }
 
-                  auto new_elem = Elem::build(TRI3);
+                  auto new_elem = Elem::build_with_id(TRI3, ne++);
                   new_elem->set_node(0) = new_node;
                   new_elem->set_node(1) = node_CW;
                   new_elem->set_node(2) = node_CCW;

src/mesh/triangulator_interface.C

@@ -128,6 +128,35 @@ void TriangulatorInterface::elems_to_segments()
 
 
 
+void TriangulatorInterface::nodes_to_segments(dof_id_type max_node_id)
+{
+  // Don't try to override manually specified segments, or try to add
+  // segments if we're doing a convex hull
+  if (!this->segments.empty() || _triangulation_type != PSLG)
+    return;
+
+  for (auto node_it = _mesh.nodes_begin(),
+       node_end = _mesh.nodes_end();
+       node_it != node_end;)
+    {
+      Node * node = *node_it;
+
+      // If we're out of boundary nodes, the rest are going to be
+      // Steiner points or hole points
+      if (node->id() >= max_node_id)
+        break;
+
+      ++node_it;
+
+      Node * next_node = (node_it == node_end) ?
+        *_mesh.nodes_begin() : *node_it;
+
+      this->segments.emplace_back(node->id(), next_node->id());
+    }
+}
+
+
+
 void TriangulatorInterface::insert_any_extra_boundary_points()
 {
   // If the initial PSLG is really simple, e.g. an L-shaped domain or
@@ -141,40 +170,38 @@ void TriangulatorInterface::insert_any_extra_boundary_points()
   const int n_interpolated = this->get_interpolate_boundary_points();
   if ((_triangulation_type==PSLG) && n_interpolated)
     {
-      // Make a copy of the original points from the Mesh
-      std::vector<Point> original_points;
-      original_points.reserve (_mesh.n_nodes());
-      for (auto & node : _mesh.node_ptr_range())
-        original_points.push_back(*node);
+      std::vector<std::pair<unsigned int, unsigned int>> old_segments =
+        std::move(this->segments);
 
-      // Clear out the mesh
-      _mesh.clear();
+      // We expect to have converted any elems and/or nodes into
+      // segments by now.
+      libmesh_assert(!old_segments.empty());
 
-      // Make sure the new Mesh will be 2D
-      _mesh.set_mesh_dimension(2);
+      this->segments.clear();
 
-      // Insert a new point on each PSLG at evenly spaced locations
+      // Insert a new point on each segment at evenly spaced locations
       // between existing boundary points.
-      // np=index into new points vector, also an explicit node index
-      // (to avoid the stride in automatic DistributedMesh indexing)
+      // np=index into new points vector
       // n =index into original points vector
-      for (std::size_t np=0, n=0, tops=n_interpolated*original_points.size(); np<tops; ++np)
+      for (auto old_segment : old_segments)
         {
-          int offset = np % n_interpolated;
-
-          // Some entries are the original points
-          if (!offset)
-            _mesh.add_point(original_points[n++], np);
-
-          else // the odd entries are equispaced along the original PSLG segments
+          Node * begin_node = _mesh.node_ptr(old_segment.first);
+          Node * end_node = _mesh.node_ptr(old_segment.second);
+          dof_id_type current_id = begin_node->id();
+          for (auto i : make_range(n_interpolated))
             {
-              libmesh_assert_less(n-1, original_points.size());
-              const std::size_t next_n = n % original_points.size();
+              // new points are equispaced along the original segments
               const Point new_point =
-                (offset*original_points[next_n] +
-                 (n_interpolated-offset)*original_points[n-1])/n_interpolated;
-              _mesh.add_point(new_point, np);
+                ((n_interpolated-i) * *(Point *)(begin_node) +
+                 (i+1) * *(Point *)(end_node)) /
+                (n_interpolated + 1);
+              Node * next_node = _mesh.add_point(new_point);
+              this->segments.emplace_back(current_id,
+                                          next_node->id());
+              current_id = next_node->id();
             }
+          this->segments.emplace_back(current_id,
+                                      end_node->id());
         }
     }
 }

tests/mesh/mesh_triangulation.C

@@ -32,7 +32,7 @@ public:
 #ifdef LIBMESH_HAVE_POLY2TRI
   CPPUNIT_TEST( testPoly2Tri );
   CPPUNIT_TEST( testPoly2TriInterp );
-  CPPUNIT_TEST( testPoly2TriInterp3 );
+  CPPUNIT_TEST( testPoly2TriInterp2 );
   CPPUNIT_TEST( testPoly2TriHoles );
   CPPUNIT_TEST( testPoly2TriMeshedHoles );
   CPPUNIT_TEST( testPoly2TriEdges );
@@ -54,7 +54,7 @@ public:
 #ifdef LIBMESH_HAVE_TRIANGLE
   CPPUNIT_TEST( testTriangle );
   CPPUNIT_TEST( testTriangleInterp );
-  CPPUNIT_TEST( testTriangleInterp3 );
+  CPPUNIT_TEST( testTriangleInterp2 );
   CPPUNIT_TEST( testTriangleHoles );
   CPPUNIT_TEST( testTriangleMeshedHoles );
   CPPUNIT_TEST( testTriangleSegments );
@@ -397,17 +397,17 @@ public:
 
     Mesh mesh(*TestCommWorld);
     TriangleInterface triangle(mesh);
-    testTriangulatorInterp(mesh, triangle, 2, 6);
+    testTriangulatorInterp(mesh, triangle, 1, 6);
   }
 
 
-  void testTriangleInterp3()
+  void testTriangleInterp2()
   {
     LOG_UNIT_TEST;
 
     Mesh mesh(*TestCommWorld);
     TriangleInterface triangle(mesh);
-    testTriangulatorInterp(mesh, triangle, 3, 10);
+    testTriangulatorInterp(mesh, triangle, 2, 10);
   }
 
 
@@ -459,17 +459,17 @@ public:
 
     Mesh mesh(*TestCommWorld);
     Poly2TriTriangulator p2t_tri(mesh);
-    testTriangulatorInterp(mesh, p2t_tri, 2, 6);
+    testTriangulatorInterp(mesh, p2t_tri, 1, 6);
   }
 
 
-  void testPoly2TriInterp3()
+  void testPoly2TriInterp2()
   {
     LOG_UNIT_TEST;
 
     Mesh mesh(*TestCommWorld);
     Poly2TriTriangulator p2t_tri(mesh);
-    testTriangulatorInterp(mesh, p2t_tri, 3, 10);
+    testTriangulatorInterp(mesh, p2t_tri, 2, 10);
   }