Merge pull request #3396 from roystgnr/hole_ordering_fix

More MeshedHole fixes

Author: roystgnr

src/mesh/mesh_triangle_holes.C

@@ -201,7 +201,9 @@ RealGradient TriangulatorInterface::Hole::areavec() const
   // the division by 2 and the norm for the end.
   //
   // Your hole points had best be coplanar, but this should work
-  // regardless of which plane they're in.
+  // regardless of which plane they're in.  If you're in the XY plane,
+  // then the standard counter-clockwise hole point ordering gives you
+  // a positive areavec(2);
 
   RealGradient areavec = 0;
 
@@ -210,7 +212,7 @@ RealGradient TriangulatorInterface::Hole::areavec() const
       const Point e_0im = this->point(i-1) - p0,
                   e_0i  = this->point(i) - p0;
 
-      areavec += e_0im.cross(e_0i);
+      areavec += e_0i.cross(e_0im);
     }
 
   return areavec;
@@ -648,13 +650,13 @@ TriangulatorInterface::MeshedHole::MeshedHole(const MeshBase & mesh,
         const Point e_0im = *hole_points[i-1] - p0,
                     e_0i  = *hole_points[i] - p0;
 
-        twice_this_area += e_0im.cross(e_0i)(2);
+        twice_this_area += e_0i.cross(e_0im)(2);
       }
 
     auto abs_twice_this_area = std::abs(twice_this_area);
 
-    if (((abs_twice_this_area == twice_this_area) && edge_type == 2) ||
-        ((abs_twice_this_area != twice_this_area) && edge_type == 1))
+    if (((twice_this_area > 0) && edge_type == 2) ||
+        ((twice_this_area < 0) && edge_type == 1))
       ++n_positive_areas;
     else
       ++n_negative_areas;
@@ -700,7 +702,7 @@ TriangulatorInterface::MeshedHole::MeshedHole(const MeshBase & mesh,
 #endif
 
   if (((twice_outer_area > 0) && outer_edge_type == 2) ||
-      outer_edge_type == 1)
+      ((twice_outer_area < 0) && outer_edge_type == 1))
     {
       if (n_positive_areas > 1)
         {

src/mesh/poly2tri_triangulator.C

@@ -1107,6 +1107,20 @@ bool Poly2TriTriangulator::insert_refinement_points()
               // "to the left as viewed from the domain interior".  We
               // need to build the updated hole ordering "backwards".
 
+              // We should never see duplicate points when we add one
+              // to a hole; if we do then we did something wrong.
+              auto push_back_new_point = [&new_points](const Point & p) {
+                // O(1) assert in devel
+                libmesh_assert(new_points.empty() ||
+                               new_points.back() != p);
+#ifdef DEBUG
+                // O(N) asserts in dbg
+                for (auto old_p : new_points)
+                  libmesh_assert_not_equal_to(old_p, p);
+#endif
+                new_points.push_back(p);
+              };
+
               for (auto point_it = arb.get_points().rbegin(),
                    point_end = arb.get_points().rend();
                    point_it != point_end;)
@@ -1117,7 +1131,7 @@ bool Poly2TriTriangulator::insert_refinement_points()
                   if (new_points.empty() ||
                       (point != new_points.back() &&
                        point != new_points.front()))
-                    new_points.push_back(point);
+                    push_back_new_point(point);
 
                   auto it = next_boundary_node.find(point);
                   while (it != next_boundary_node.end())
@@ -1128,7 +1142,7 @@ bool Poly2TriTriangulator::insert_refinement_points()
                       if (point_it != point_end &&
                           point == *point_it)
                         ++point_it;
-                      new_points.push_back(point);
+                      push_back_new_point(point);
                       it = next_boundary_node.find(point);
                     }
                 }

src/mesh/triangulator_interface.C

@@ -112,8 +112,7 @@ void TriangulatorInterface::elems_to_segments()
         }
 
       // We'll steal the ordering calculation from
-      // the MeshedHole code, but reverse the ordering since it's
-      // to be used as an outer rather than an inner boundary.
+      // the MeshedHole code
       const TriangulatorInterface::MeshedHole mh { _mesh, this->_bdy_ids };
 
       // And now we're done with elements.  Delete them lest they have
@@ -125,7 +124,6 @@ void TriangulatorInterface::elems_to_segments()
       // inside that boundary.  Triangulator code doesn't like Steiner
       // points that aren't inside the triangulation domain, so we
       // need to get rid of them.
-
       std::unordered_set<Node *> nodes_to_delete;
       if (!this->_bdy_ids.empty())
         {
@@ -139,10 +137,10 @@ void TriangulatorInterface::elems_to_segments()
       const std::size_t np = mh.n_points();
       for (auto i : make_range(np))
         {
-          const Point pt = mh.point(np-i-1);
+          const Point pt = mh.point(i);
           const dof_id_type id0 = libmesh_map_find(point_id_map, pt);
           nodes_to_delete.erase(_mesh.node_ptr(id0));
-          const Point next_pt = mh.point((2*np-i-2)%np);
+          const Point next_pt = mh.point((np+i+1)%np);
           const dof_id_type id1 = libmesh_map_find(point_id_map, next_pt);
           this->segments.emplace_back(id0, id1);
         }

tests/mesh/mesh_triangulation.C

@@ -904,7 +904,7 @@ public:
 
     // If refinement should have increased our element count, check it
     if (desired_area || area_func)
-      CPPUNIT_ASSERT_LESS(mesh.n_elem(), n_original_elem); // n_elem+++
+      CPPUNIT_ASSERT_GREATER(n_original_elem, mesh.n_elem()); // n_elem+++
     else
       CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_original_elem);