Add more integrity checks for tetrahedralization

include/mesh/mesh_tet_interface.h

@@ -92,8 +92,23 @@ class MeshTetInterface
    */
   virtual void triangulate () = 0;
 
+  /**
+   * Sets a verbosity level, defaulting to 0 (print nothing), to be
+   * set as high as 100 (print everything).
+   *
+   * For verbosity >= 50, print all detected surface mesh integrity
+   * issues as they're found.  Subclasses may add other output at
+   * other verbosity levels.
+   */
+  void set_verbosity(unsigned int v);
+
 protected:
 
+  /**
+   * verbosity setting
+   */
+  unsigned int _verbosity;
+
   /**
    * Remove volume elements from the given mesh, after converting
    * their outer boundary faces to surface elements.
@@ -103,30 +118,52 @@ 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
+    DEGENERATE_ELEMENT = 8, // an element has zero area
+    DEGENERATE_MESH = 9     // the mesh clearly bounds zero volume
+  };
+
   /**
    * 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.
    *
-   * \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
+   * Named \p check_hull_integrity() for backward compatibility, but
+   * now accepts non-convex manifolds.
+   *
+   * \returns a set of \p  enums describing problems
+   * found, or an empty set if no problems are found.
+   */
+  [[nodiscard]] std::set<SurfaceIntegrity> check_hull_integrity() const;
+
+  /**
+   * This function checks the integrity of the current set of
+   * elements in the Mesh, and corrects what it can.
+   *
+   * \returns A set of SurfaceIntegrity codes from \p
+   * check_hull_integrity() if there are problems it can't fix, or an
+   * empty set otherwise.
    */
-  [[nodiscard]] unsigned int check_hull_integrity();
+  [[nodiscard]] std::set<SurfaceIntegrity> improve_hull_integrity();
 
   /**
    * 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
@@ -163,6 +200,17 @@ class MeshTetInterface
   std::unique_ptr<std::vector<std::unique_ptr<UnstructuredMesh>>> _holes;
 };
 
+
+// ------------------------------------------------------------
+// MeshTetInterface class member functions
+inline
+void
+MeshTetInterface::set_verbosity(unsigned int v)
+{
+  this->_verbosity = v;
+}
+
+
 } // namespace libMesh
 
 #endif // LIBMESH_MESH_TET_INTERFACE_H

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.
@@ -123,13 +128,16 @@ void NetGenMeshInterface::triangulate ()
       // Receive the mesh data rank 0 will send later, then fix it up
       // together
       MeshCommunication().broadcast(this->_mesh);
+
+      // If we got an empty mesh here then our tetrahedralization
+      // failed.
+      libmesh_error_msg_if (!this->_mesh.n_elem(),
+                            "NetGen failed to generate any tetrahedra");
+
       this->_mesh.prepare_for_use();
       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
@@ -324,8 +332,17 @@ void NetGenMeshInterface::triangulate ()
 
   const int n_elem = Ng_GetNE(ngmesh);
 
-  libmesh_error_msg_if (n_elem <= 0,
-                        "NetGen failed to generate any tetrahedra");
+  // If Netgen fails us, we're likely to get n_elem <= 0.  This is a
+  // common enough failure from bad setups that I want to make sure
+  // it's thrown in parallel so as to not desynchronize any unit tests
+  // that trigger it.  So we'll broadcast the empty mesh to indicate
+  // the problem and enable throwing exceptions in parallel.
+  if (n_elem <= 0)
+    {
+      this->_mesh.clear();
+      MeshCommunication().broadcast(this->_mesh);
+      libmesh_error_msg ("NetGen failed to generate any tetrahedra");
+    }
 
   const dof_id_type n_points = Ng_GetNP(ngmesh);
   const dof_id_type old_nodes = this->_mesh.n_nodes();