From 7a0714836027037720ad97c266e09be44cbdfb0e Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 23 Feb 2026 20:59:56 -0600 Subject: [PATCH 1/9] 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. --- include/mesh/mesh_tet_interface.h | 32 +++++--- src/mesh/mesh_netgen_interface.C | 8 +- src/mesh/mesh_tet_interface.C | 122 ++++++++++++++++++++++-------- src/mesh/mesh_tetgen_interface.C | 2 +- 4 files changed, 119 insertions(+), 45 deletions(-) diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index 013b920c30b..c8297e374c2 100644 --- a/include/mesh/mesh_tet_interface.h +++ b/include/mesh/mesh_tet_interface.h @@ -103,22 +103,32 @@ 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 check_hull_integrity() const; /** * This function prints an informative message and throws an @@ -126,7 +136,7 @@ class MeshTetInterface * 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 & result) const; /** * Delete original convex hull elements from the Mesh diff --git a/src/mesh/mesh_netgen_interface.C b/src/mesh/mesh_netgen_interface.C index 742d289c042..bc26d5bbbe5 100644 --- a/src/mesh/mesh_netgen_interface.C +++ b/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 diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index 0d72d8aaa57..d7e0181a486 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/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::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 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 & 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()); } diff --git a/src/mesh/mesh_tetgen_interface.C b/src/mesh/mesh_tetgen_interface.C index a629ff83239..76cce90cd44 100644 --- a/src/mesh/mesh_tetgen_interface.C +++ b/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); From 3b36a9cb720ec0b786f71d4efcbe38296184ef9d Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 23 Feb 2026 21:00:31 -0600 Subject: [PATCH 2/9] Fail in parallel when we've confused NetGen This might be useful in letting application code recover from such errors, and it's definitely useful in debugging unit tests where I've screwed something up. --- src/mesh/mesh_netgen_interface.C | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_netgen_interface.C b/src/mesh/mesh_netgen_interface.C index bc26d5bbbe5..c8b3affd711 100644 --- a/src/mesh/mesh_netgen_interface.C +++ b/src/mesh/mesh_netgen_interface.C @@ -128,6 +128,12 @@ 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; } @@ -326,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(); From 3607fca62fcc652436fc39d562a458af2b71addf Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 23 Feb 2026 21:01:20 -0600 Subject: [PATCH 3/9] Try to fix hull problems, not just diagnose them --- include/mesh/mesh_tet_interface.h | 10 +++ src/mesh/mesh_tet_interface.C | 136 ++++++++++++++++++++++++++++++ 2 files changed, 146 insertions(+) diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index c8297e374c2..5673747b85c 100644 --- a/include/mesh/mesh_tet_interface.h +++ b/include/mesh/mesh_tet_interface.h @@ -130,6 +130,16 @@ class MeshTetInterface */ [[nodiscard]] std::set 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]] std::set improve_hull_integrity(); + /** * This function prints an informative message and throws an * exception based on the output of the check_hull_integrity() diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index d7e0181a486..67cdf4090f7 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/src/mesh/mesh_tet_interface.C @@ -408,6 +408,142 @@ std::set MeshTetInterface::check_hull_integr +std::set MeshTetInterface::improve_hull_integrity() +{ + // We don't really do anything parallel here, but we aspire to. + libmesh_parallel_only(this->_mesh.comm()); + + std::set integrityproblems = + this->check_hull_integrity(); + + // If we have no problem, or a problem we can't fix, we're done. + if (integrityproblems.empty() || + integrityproblems.count(NON_TRI3) || + integrityproblems.count(EMPTY_MESH)) + return integrityproblems; + + // Possibly the user gave us an unprepared mesh with missing or bad + // neighbor links? + if (integrityproblems.count(MISSING_NEIGHBOR) || + integrityproblems.count(MISSING_BACKLINK) || + integrityproblems.count(BAD_NEIGHBOR_LINKS)) + { + this->_mesh.find_neighbors(); + integrityproblems = this->check_hull_integrity(); + } + + // If find_neighbors() doesn't fix these, I give up. + if (integrityproblems.count(MISSING_NEIGHBOR) || + integrityproblems.count(MISSING_BACKLINK) || + integrityproblems.count(BAD_NEIGHBOR_LINKS)) + return integrityproblems; + + // find_neighbors() might have fixed everything + if (integrityproblems.empty()) + return integrityproblems; + + // A non-oriented (but orientable!) surface is the only thing we + // shouldn't have fixed or given up on by now. + libmesh_assert_equal_to(integrityproblems.size(), 1); + libmesh_assert_equal_to(integrityproblems.count(NON_ORIENTED), 1); + + // We need one known-good triangle to start from. We'll pick the + // most-negative-x normal among the triangles on the most-negative-x + // point. + + // We'll just implement this in serial for now. + MeshSerializer mesh_serializer(this->_mesh); + + const Node * lowest_point = (*this->_mesh.elements_begin())->node_ptr(0); + + // Index by ids, not pointers, for consistency in parallel + std::unordered_set attached_elements; + + for (Elem * elem : this->_mesh.element_ptr_range()) + { + for (const Node & node : elem->node_ref_range()) + { + if (node(0) < (*lowest_point)(0)) + { + lowest_point = &node; + attached_elements.clear(); + } + if (&node == lowest_point) + attached_elements.insert(elem->id()); + } + } + + Elem * best_elem = nullptr; + Real best_abs_normal_0; + + for (dof_id_type id : attached_elements) + { + Elem * elem = this->_mesh.elem_ptr(id); + const Point e01 = elem->point(1) - elem->point(0); + const Point e02 = elem->point(2) - elem->point(0); + const Point normal = e01.cross(e02).unit(); + const Real abs_normal_0 = std::abs(normal(0)); + + if (!best_elem || abs_normal_0 > best_abs_normal_0) + { + best_elem = elem; + best_abs_normal_0 = abs_normal_0; + } + } + + // Now flood-fill from that element to get a consistent orientation + // for the others. + std::unordered_set frontier_elements{best_elem->id()}, + finished_elements{}; + + BoundaryInfo & bi = this->_mesh.get_boundary_info(); + + while (!frontier_elements.empty()) + { + const dof_id_type elem_id = *frontier_elements.begin(); + Elem & elem = this->_mesh.elem_ref(elem_id); + for (auto s : elem.side_index_range()) + { + Elem * neigh = elem.neighbor_ptr(s); + libmesh_assert(neigh); + libmesh_assert_less(neigh->which_neighbor_am_i(&elem), 3); + + 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()); + libmesh_assert_less(i1, 3); + libmesh_assert_less(i2, 3); + + const dof_id_type neigh_id = neigh->id(); + + const bool frontier_neigh = frontier_elements.count(neigh_id); + const bool finished_neigh = finished_elements.count(neigh_id); + + // Are we flipped? + if ((i2 + 1)%3 != i1) + { + // Are we a Moebius strip??? We give up. + if (frontier_neigh || finished_neigh) + return integrityproblems; + + neigh->flip(&bi); + } + + if (!frontier_neigh && !finished_neigh) + frontier_elements.insert(neigh_id); + } + + finished_elements.insert(elem_id); + frontier_elements.erase(elem_id); + } + + this->_mesh.find_neighbors(); + + libmesh_assert(this->check_hull_integrity().empty()); + + return {}; +} void MeshTetInterface::process_hull_integrity_result From 55a5fefe6d0bea4ee31a295ba7eca8a8a9cbcdd2 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Mon, 23 Feb 2026 20:41:29 -0600 Subject: [PATCH 4/9] Unit test coverage of badly oriented tet boundary --- tests/mesh/mesh_tet_test.C | 31 ++++++++++++++++++++++++++++++- 1 file changed, 30 insertions(+), 1 deletion(-) diff --git a/tests/mesh/mesh_tet_test.C b/tests/mesh/mesh_tet_test.C index 197b7d63263..f116294296f 100644 --- a/tests/mesh/mesh_tet_test.C +++ b/tests/mesh/mesh_tet_test.C @@ -50,6 +50,7 @@ public: CPPUNIT_TEST( testNetGenError ); CPPUNIT_TEST( testNetGenTets ); CPPUNIT_TEST( testNetGenFlippedTris ); + CPPUNIT_TEST( testNetGenNonOriented ); CPPUNIT_TEST( testNetGenHole ); #ifdef LIBMESH_ENABLE_AMR @@ -232,13 +233,30 @@ public: void testTrisToTets(UnstructuredMesh & mesh, MeshTetInterface & triangulator, - bool flip_tris = false) + bool flip_tris = false, + bool flip_some_tris = false) { // An asymmetric octahedron, so we hopefully have an unambiguous // choice of shortest diagonal for a Delaunay algorithm to pick. const Real expected_volume = build_octahedron(mesh, flip_tris, -1, 1, -1, 1, -0.1, 0.1); + // Flip a couple tri, breaking the mesh in a way we can fix + if (flip_some_tris) + for (auto elem : mesh.element_ptr_range()) + { + Point center = elem->vertex_average(); + if ((center(0) > 0 && + center(1) > 0 && + center(2) > 0) || + (center(0) < 0 && + center(1) < 0 && + center(2) < 0)) + elem->flip(&mesh.get_boundary_info()); + + mesh.unset_is_prepared(); + } + this->testTetInterfaceBase(mesh, triangulator, /* n_elem = */ 4, /* n_nodes = */ 6, expected_volume); } @@ -394,6 +412,17 @@ public: } + void testNetGenNonOriented() + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + NetGenMeshInterface net_tet(mesh); + testTrisToTets(mesh, net_tet, true, true); + testBcids(mesh); + } + + void testNetGenHole() { LOG_UNIT_TEST; From f5c7f9a2ddabfbb301a524f15c022edbab1b3fd9 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 24 Feb 2026 10:51:55 -0600 Subject: [PATCH 5/9] Double-check that "best" element in boundary prep NetGen is only happy with oriented manifolds in one orientation --- src/mesh/mesh_tet_interface.C | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index 67cdf4090f7..c17724c00a1 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/src/mesh/mesh_tet_interface.C @@ -454,6 +454,10 @@ std::set MeshTetInterface::improve_hull_inte // We'll just implement this in serial for now. MeshSerializer mesh_serializer(this->_mesh); + // I don't see why we'd need boundary info here, but maybe we'll + // want to preserve edge/node conditions eventually? + BoundaryInfo & bi = this->_mesh.get_boundary_info(); + const Node * lowest_point = (*this->_mesh.elements_begin())->node_ptr(0); // Index by ids, not pointers, for consistency in parallel @@ -488,6 +492,11 @@ std::set MeshTetInterface::improve_hull_inte { best_elem = elem; best_abs_normal_0 = abs_normal_0; + + // Make sure that element is actually a good one, by + // flipping it if it's not. + if (abs_normal_0 == normal(0)) + elem->flip(&bi); } } @@ -496,8 +505,6 @@ std::set MeshTetInterface::improve_hull_inte std::unordered_set frontier_elements{best_elem->id()}, finished_elements{}; - BoundaryInfo & bi = this->_mesh.get_boundary_info(); - while (!frontier_elements.empty()) { const dof_id_type elem_id = *frontier_elements.begin(); From 9038a112855eefe8e8164451472f58b6b5835001 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 24 Feb 2026 16:24:56 -0600 Subject: [PATCH 6/9] Added MeshTetInterface verbosity option --- include/mesh/mesh_tet_interface.h | 15 +++++++++++++++ src/mesh/mesh_tet_interface.C | 18 ++++++++++++++++-- 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index 5673747b85c..69c01675237 100644 --- a/include/mesh/mesh_tet_interface.h +++ b/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 + */ + const unsigned int _verbosity; + /** * Remove volume elements from the given mesh, after converting * their outer boundary faces to surface elements. diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index c17724c00a1..2a24bb9e4bf 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/src/mesh/mesh_tet_interface.C @@ -112,7 +112,7 @@ namespace libMesh //---------------------------------------------------------------------- // MeshTetInterface class members MeshTetInterface::MeshTetInterface (UnstructuredMesh & mesh) : - _desired_volume(0), _smooth_after_generating(false), + _verbosity(0), _desired_volume(0), _smooth_after_generating(false), _elem_type(TET4), _mesh(mesh) { } @@ -350,7 +350,11 @@ std::set MeshTetInterface::check_hull_integr { // Check for proper element type if (elem->type() != TRI3) - returnval.insert(NON_TRI3); + { + if (this->_verbosity >= 50) + std::cerr << "Non-Tri3: " << elem->get_info() << std::endl; + returnval.insert(NON_TRI3); + } for (auto s : elem->side_index_range()) { @@ -358,6 +362,8 @@ std::set MeshTetInterface::check_hull_integr if (neigh == nullptr) { + if (this->_verbosity >= 50) + std::cerr << "Element missing neighbor " << s << ": " << elem->get_info() << std::endl; returnval.insert(MISSING_NEIGHBOR); continue; } @@ -367,6 +373,8 @@ std::set MeshTetInterface::check_hull_integr if (nn >= 3) { + if (this->_verbosity >= 50) + std::cerr << "Element missing backlink " << s << ": " << elem->get_info() << std::endl; returnval.insert(MISSING_BACKLINK); continue; } @@ -380,6 +388,8 @@ std::set MeshTetInterface::check_hull_integr const unsigned int i2 = neigh->local_node(n2->id()); if (i1 >= 3 || i2 >= 3) { + if (this->_verbosity >= 50) + std::cerr << "Element with bad neighbor " << s << " nodes: " << elem->get_info() << std::endl; returnval.insert(BAD_NEIGHBOR_NODES); continue; } @@ -388,6 +398,8 @@ std::set MeshTetInterface::check_hull_integr // (because they have the same orientation we do) if ((i2 + 1)%3 != i1) { + if (this->_verbosity >= 50) + std::cerr << "Element orientation mismatch with neighbor " << s << ": " << elem->get_info() << std::endl; returnval.insert(NON_ORIENTED); continue; } @@ -396,6 +408,8 @@ std::set MeshTetInterface::check_hull_integr // places relative to its neighbor link if (i2 != nn) { + if (this->_verbosity >= 50) + std::cerr << "Element with bad links on neighbor " << s << ": " << elem->get_info() << std::endl; returnval.insert(BAD_NEIGHBOR_LINKS); continue; } From b05519ff7702d93817253b847a23a07e9c3f9806 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Tue, 24 Feb 2026 16:25:28 -0600 Subject: [PATCH 7/9] MeshTetInterface degeneracy detection --- include/mesh/mesh_tet_interface.h | 4 +++- src/mesh/mesh_tet_interface.C | 21 +++++++++++++++++++++ 2 files changed, 24 insertions(+), 1 deletion(-) diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index 69c01675237..7a85cf9a12b 100644 --- a/include/mesh/mesh_tet_interface.h +++ b/include/mesh/mesh_tet_interface.h @@ -128,7 +128,9 @@ class MeshTetInterface 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 + 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 }; /** diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index 2a24bb9e4bf..5a7fcb29a1c 100644 --- a/src/mesh/mesh_tet_interface.C +++ b/src/mesh/mesh_tet_interface.C @@ -346,6 +346,19 @@ std::set MeshTetInterface::check_hull_integr std::set returnval; + const BoundingBox bb = MeshTools::create_bounding_box(this->_mesh); + const Point extents = bb.max() - bb.min(); + if (extents(0) == 0 || + extents(1) == 0 || + extents(2) == 0) + returnval.insert(DEGENERATE_MESH); + + // Figure a area to use for relative tolerances when detecting + // degenerate elements + const Real ref_area = std::abs(extents(0) * extents(1)) + + std::abs(extents(0) * extents(2)) + + std::abs(extents(1) * extents(2)); + for (auto & elem : this->_mesh.element_ptr_range()) { // Check for proper element type @@ -356,6 +369,14 @@ std::set MeshTetInterface::check_hull_integr returnval.insert(NON_TRI3); } + // Make sure it's a decent element. + if (elem->volume() < ref_area * TOLERANCE * TOLERANCE) + { + if (this->_verbosity >= 50) + std::cerr << "Degenerate element: " << elem->get_info() << std::endl; + returnval.insert(DEGENERATE_ELEMENT); + } + for (auto s : elem->side_index_range()) { const Elem * const neigh = elem->neighbor_ptr(s); From 4c06a3e844666db0437d877c29ed0b43a1908597 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 25 Feb 2026 10:18:42 -0600 Subject: [PATCH 8/9] Actually *define* verbosity setter --- include/mesh/mesh_tet_interface.h | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index 7a85cf9a12b..588afd43f73 100644 --- a/include/mesh/mesh_tet_interface.h +++ b/include/mesh/mesh_tet_interface.h @@ -107,7 +107,7 @@ class MeshTetInterface /** * verbosity setting */ - const unsigned int _verbosity; + unsigned int _verbosity; /** * Remove volume elements from the given mesh, after converting @@ -200,6 +200,17 @@ class MeshTetInterface std::unique_ptr>> _holes; }; + +// ------------------------------------------------------------ +// MeshTetInterface class member functions +inline +void +MeshTetInterface::set_verbosity(unsigned int v) +{ + this->_verbosity = v; +} + + } // namespace libMesh #endif // LIBMESH_MESH_TET_INTERFACE_H From 9435cce20b954147f842538b67a40599155c1dc7 Mon Sep 17 00:00:00 2001 From: Roy Stogner Date: Wed, 25 Feb 2026 10:18:55 -0600 Subject: [PATCH 9/9] Call verbosity setter You'd think that adding one line of do-nothing test coverage is pointless, but at least it would have caught that embarrassing omitted definition as a linker error. --- tests/mesh/mesh_tet_test.C | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/tests/mesh/mesh_tet_test.C b/tests/mesh/mesh_tet_test.C index f116294296f..4472defdfcb 100644 --- a/tests/mesh/mesh_tet_test.C +++ b/tests/mesh/mesh_tet_test.C @@ -375,6 +375,10 @@ public: Mesh mesh(*TestCommWorld); NetGenMeshInterface net_tet(mesh); + + // This should give no output for our correctly-set-up inputs + net_tet.set_verbosity(50); + testTrisToTets(mesh, net_tet); testBcids(mesh); }