diff --git a/include/mesh/mesh_tet_interface.h b/include/mesh/mesh_tet_interface.h index 013b920c30..588afd43f7 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 + */ + unsigned int _verbosity; + /** * Remove volume elements from the given mesh, after converting * their outer boundary faces to surface elements. @@ -103,22 +118,44 @@ 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 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 improve_hull_integrity(); /** * This function prints an informative message and throws an @@ -126,7 +163,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 @@ -163,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 diff --git a/src/mesh/mesh_netgen_interface.C b/src/mesh/mesh_netgen_interface.C index 742d289c04..c8b3affd71 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. @@ -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(); diff --git a/src/mesh/mesh_tet_interface.C b/src/mesh/mesh_tet_interface.C index 0d72d8aaa5..5a7fcb29a1 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) { } @@ -335,65 +335,305 @@ 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; + + 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 if (elem->type() != TRI3) { - //libmesh_error_msg("ERROR: Some of the elements in the original mesh were not TRI3!"); - return 1; + if (this->_verbosity >= 50) + std::cerr << "Non-Tri3: " << elem->get_info() << std::endl; + 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 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; + if (this->_verbosity >= 50) + std::cerr << "Element missing neighbor " << s << ": " << elem->get_info() << std::endl; + 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) + { + if (this->_verbosity >= 50) + std::cerr << "Element missing backlink " << s << ": " << elem->get_info() << std::endl; + 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) + { + if (this->_verbosity >= 50) + std::cerr << "Element with bad neighbor " << s << " nodes: " << elem->get_info() << std::endl; + 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) + { + if (this->_verbosity >= 50) + std::cerr << "Element orientation mismatch with neighbor " << s << ": " << elem->get_info() << std::endl; + returnval.insert(NON_ORIENTED); + continue; + } + + // And it should have those edge nodes in the expected + // 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; } } } - // 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::set MeshTetInterface::improve_hull_integrity() { - std::ostringstream err_msg; + // 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); + + // 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(); - if (result != 0) + 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()) { - err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl; + 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()); + } + } - if (result==1) + 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) { - 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; + 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); } + } - if (result==2) + // Now flood-fill from that element to get a consistent orientation + // for the others. + std::unordered_set frontier_elements{best_elem->id()}, + finished_elements{}; + + 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()) { - 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; + 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); } - if (result==3) - err_msg << "The input Mesh was empty!" << std::endl; + finished_elements.insert(elem_id); + frontier_elements.erase(elem_id); + } + + this->_mesh.find_neighbors(); + + libmesh_assert(this->check_hull_integrity().empty()); + + return {}; +} + - libmesh_error_msg(err_msg.str()); +void MeshTetInterface::process_hull_integrity_result + (const std::set & result) const +{ + std::ostringstream err_msg; + + if (result.empty()) // success + return; + + err_msg << "Error! Conforming Delaunay mesh tetrahedralization requires a convex hull." << std::endl; + + 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 a629ff8323..76cce90cd4 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); diff --git a/tests/mesh/mesh_tet_test.C b/tests/mesh/mesh_tet_test.C index 197b7d6326..4472defdfc 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); } @@ -357,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); } @@ -394,6 +416,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;