Use a set of enums for tet hull integrity

An integer error code was kind of lazy of us, and using a set lets us
cleanly fix errors rather than just report them.

Author: roystgnr

include/mesh/mesh_tet_interface.h

@@ -103,30 +103,40 @@ class MeshTetInterface
    */
   static BoundingBox volume_to_surface_mesh (UnstructuredMesh & mesh);
 
+  /**
+   * Enumeration of possible surface mesh integrity issues.
+   */
+  enum SurfaceIntegrity {
+    NON_TRI3 = 1,           // a non-TRI3 element is found
+    MISSING_NEIGHBOR = 2,   // an element with a nullptr-neighbor is found
+    EMPTY_MESH = 3,         // the mesh is empty
+    MISSING_BACKLINK = 4,   // an element neighbor isn't linked back to it
+    BAD_NEIGHBOR_NODES = 5, // an element neighbor isn't linked to expected nodes
+    NON_ORIENTED = 6,       // an element neighbor has inconsistent orientation
+    BAD_NEIGHBOR_LINKS = 7  // an element neighbor has other inconsistent links
+  };
+
   /**
    * This function checks the integrity of the current set of
    * elements in the Mesh to see if they comprise a topological
    * manifold that (if it's also geometrically valid) would define
-   * valid hull for a tetrahedralized volume.
-   * That is:
-   * - If they are all TRI3 elements
-   * - They all have non-nullptr neighbors
+   * valid boundary for a tetrahedralized volume.
+   *
+   * Named \p check_hull_integrity() for backward compatibility, but
+   * now accepts non-convex manifolds.
    *
-   * \returns
-   * - 0 if the mesh forms a topologically valid hull
-   * - 1 if a non-TRI3 element is found
-   * - 2 if an element with a nullptr-neighbor is found
-   * - 3 if the mesh is empty
+   * \returns a set of \p  enums describing problems
+   * found, or an empty set if no problems are found.
    */
-  [[nodiscard]] unsigned int check_hull_integrity();
+  [[nodiscard]] std::set<SurfaceIntegrity> check_hull_integrity() const;
 
   /**
    * This function prints an informative message and throws an
    * exception based on the output of the check_hull_integrity()
    * function.  It is a separate function so that you can check hull
    * integrity without exiting or catching an exception if desired.
    */
-  void process_hull_integrity_result(unsigned int result);
+  void process_hull_integrity_result(const std::set<SurfaceIntegrity> & result) const;
 
   /**
    * Delete original convex hull elements from the Mesh

src/mesh/mesh_netgen_interface.C

@@ -109,6 +109,11 @@ void NetGenMeshInterface::triangulate ()
            /* serial_only_needed_on_proc_0 */ true);
       }
 
+  // This should probably only be done on rank 0, but the API is
+  // designed with the hope that we'll parallelize it eventually
+  auto integrity = this->improve_hull_integrity();
+  this->process_hull_integrity_result(integrity);
+
   // If we're not rank 0, we're just going to wait for rank 0 to call
   // Netgen, then receive its data afterward, we're not going to hope
   // that Netgen does the exact same thing on every processor.
@@ -127,9 +132,6 @@ void NetGenMeshInterface::triangulate ()
       return;
     }
 
-  auto integrity = this->check_hull_integrity();
-  this->process_hull_integrity_result(integrity);
-
   Ng_Meshing_Parameters params;
 
   // Override any default parameters we might need to, to avoid

src/mesh/mesh_tet_interface.C

@@ -335,65 +335,127 @@ BoundingBox MeshTetInterface::volume_to_surface_mesh(UnstructuredMesh & mesh)
 }
 
 
-unsigned int MeshTetInterface::check_hull_integrity()
+std::set<MeshTetInterface::SurfaceIntegrity> MeshTetInterface::check_hull_integrity() const
 {
   // Check for easy return: if the Mesh is empty (i.e. if
   // somebody called triangulate_conformingDelaunayMesh on
   // a Mesh with no elements, then hull integrity check must
   // fail...
   if (_mesh.n_elem() == 0)
-    return 3;
+    return {EMPTY_MESH};
+
+  std::set<MeshTetInterface::SurfaceIntegrity> returnval;
 
   for (auto & elem : this->_mesh.element_ptr_range())
     {
       // Check for proper element type
       if (elem->type() != TRI3)
-        {
-          //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!");
-          return 1;
-        }
+        returnval.insert(NON_TRI3);
 
-      for (auto neigh : elem->neighbor_ptr_range())
+      for (auto s : elem->side_index_range())
         {
+          const Elem * const neigh = elem->neighbor_ptr(s);
+
           if (neigh == nullptr)
             {
-              // libmesh_error_msg("ERROR: Non-convex hull, cannot be tetrahedralized.");
-              return 2;
+              returnval.insert(MISSING_NEIGHBOR);
+              continue;
+            }
+
+          // Make sure our neighbor points back to us
+          const unsigned int nn = neigh->which_neighbor_am_i(elem);
+
+          if (nn >= 3)
+            {
+              returnval.insert(MISSING_BACKLINK);
+              continue;
+            }
+
+          // Our neighbor should have the same the edge nodes we do on
+          // the neighboring edgei
+          const Node * const n1 = elem->node_ptr(s);
+          const Node * const n2 = elem->node_ptr((s+1)%3);
+
+          const unsigned int i1 = neigh->local_node(n1->id());
+          const unsigned int i2 = neigh->local_node(n2->id());
+          if (i1 >= 3 || i2 >= 3)
+            {
+              returnval.insert(BAD_NEIGHBOR_NODES);
+              continue;
+            }
+
+          // It should have those edge nodes in the opposite order
+          // (because they have the same orientation we do)
+          if ((i2 + 1)%3 != i1)
+            {
+              returnval.insert(NON_ORIENTED);
+              continue;
+            }
+
+          // And it should have those edge nodes in the expected
+          // places relative to its neighbor link
+          if (i2 != nn)
+            {
+              returnval.insert(BAD_NEIGHBOR_LINKS);
+              continue;
             }
         }
     }
 
-  // If we made it here, return success!
-  return 0;
+  // Return anything and everything we found
+  return returnval;
 }
 
 
 
-void MeshTetInterface::process_hull_integrity_result(unsigned result)
-{
-  std::ostringstream err_msg;
 
-  if (result != 0)
-    {
-      err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
 
-      if (result==1)
-        {
-          err_msg << "Non-TRI3 elements were found in the input Mesh.  ";
-          err_msg << "A constrained Delaunay tetrahedralization requires a convex hull of TRI3 elements." << std::endl;
-        }
+void MeshTetInterface::process_hull_integrity_result
+  (const std::set<MeshTetInterface::SurfaceIntegrity> & result) const
+{
+  std::ostringstream err_msg;
 
-      if (result==2)
-        {
-          err_msg << "At least one triangle without three neighbors was found in the input Mesh.  ";
-          err_msg << "A constrained Delaunay tetrahedralization must be a triangular manifold without boundary." << std::endl;
-        }
+  if (result.empty()) // success
+    return;
 
-      if (result==3)
-        err_msg << "The input Mesh was empty!" << std::endl;
+  err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl;
 
-      libmesh_error_msg(err_msg.str());
+  if (result.count(NON_TRI3))
+    {
+      err_msg << "At least one non-Tri3 element was found in the input boundary mesh.  ";
+      err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
+    }
+  if (result.count(MISSING_NEIGHBOR))
+    {
+      err_msg << "At least one triangle without three neighbors was found in the input boundary mesh.  ";
+      err_msg << "A constrained Delaunay tetrahedralization boundary must be a triangular manifold without boundary." << std::endl;
+    }
+  if (result.count(EMPTY_MESH))
+    {
+      err_msg << "The input boundary mesh was empty!" << std::endl;
+      err_msg << "Our constrained Delaunay tetrahedralization boundary must be a triangulation of Tri3 elements." << std::endl;
+    }
+  if (result.count(MISSING_BACKLINK))
+    {
+      err_msg << "At least one triangle neighbor without a return neighbor link was found in the input boundary mesh.  ";
+      err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
+    }
+  if (result.count(BAD_NEIGHBOR_NODES))
+    {
+      err_msg << "At least one triangle neighbor without expected node links was found in the input boundary mesh.  ";
+      err_msg << "A constrained Delaunay tetrahedralization boundary must be a conforming and non-adaptively-refined mesh." << std::endl;
     }
+  if (result.count(NON_ORIENTED))
+    {
+      err_msg << "At least one triangle neighbor with an inconsistent orientation was found in the input boundary mesh.  ";
+      err_msg << "A constrained Delaunay tetrahedralization boundary must be an oriented Tri3 mesh." << std::endl;
+    }
+  if (result.count(BAD_NEIGHBOR_LINKS))
+    {
+      err_msg << "At least one triangle neighbor with inconsistent node and neighbor links was found in the input boundary mesh." << std::endl;
+    }
+
+  libmesh_error_msg(err_msg.str());
 }
 
 

src/mesh/mesh_tetgen_interface.C

@@ -202,7 +202,7 @@ void TetGenMeshInterface::triangulate_conformingDelaunayMesh_carvehole  (const s
 {
   // Before calling this function, the Mesh must contain a convex hull
   // of TRI3 elements which define the boundary.
-  unsigned hull_integrity_check = check_hull_integrity();
+  auto hull_integrity_check = this->improve_hull_integrity();
 
   // Possibly die if hull integrity check failed
   this->process_hull_integrity_result(hull_integrity_check);