diff --git a/src/mesh/vtk_io.C b/src/mesh/vtk_io.C index d479c6e4c5..85e684aa69 100644 --- a/src/mesh/vtk_io.C +++ b/src/mesh/vtk_io.C @@ -25,6 +25,9 @@ #include "libmesh/system.h" #include "libmesh/node.h" #include "libmesh/elem.h" +#include "libmesh/cell_c0polyhedron.h" +#include "libmesh/face_c0polygon.h" +#include "libmesh/face_polygon.h" #include "libmesh/enum_io_package.h" #include "libmesh/utility.h" @@ -38,12 +41,14 @@ #include "vtkXMLPUnstructuredGridReader.h" #include "vtkXMLUnstructuredGridWriter.h" #include "vtkXMLPUnstructuredGridWriter.h" +#include "vtkXMLFileReadTester.h" #include "vtkUnstructuredGrid.h" #include "vtkIntArray.h" #include "vtkCellArray.h" #include "vtkCellData.h" #include "vtkDoubleArray.h" #include "vtkGenericCell.h" +#include "vtkIdList.h" #include "vtkPointData.h" #include "vtkPoints.h" #include "vtkSmartPointer.h" @@ -57,6 +62,7 @@ #include "libmesh/restore_warnings.h" // C++ includes +#include #include @@ -183,6 +189,72 @@ std::map VTKIO::build_element_maps() +namespace { + +// VTK describes a polyhedron by a "face stream" (a list of polygonal +// faces, each a list of global point ids) rather than by the flat, +// ordered node list used for every other supported element type. This +// helper builds an equivalent libMesh C0Polyhedron from that face +// stream, assigning its node pointers and adding to \p mesh any interior +// "mid-element" node the C0Polyhedron construction requires. +std::unique_ptr +add_vtk_polyhedron(vtkUnstructuredGrid & vtk_grid, + MeshBase & mesh, + const vtkIdType cell_id, + vtkIntArray * node_id, + const std::vector & vtk_node_to_libmesh) +{ + // GetFaceStream() fills face_stream with + // [ n_faces, + // n_face0_pts, id0, id1, ..., + // n_face1_pts, id0, id1, ..., ] + // where all the point ids are global VTK point ids. + vtkSmartPointer face_stream = vtkSmartPointer::New(); + vtk_grid.GetFaceStream(cell_id, face_stream); + + vtkIdType pos = 0; + const vtkIdType n_faces = face_stream->GetId(pos++); + + libmesh_error_msg_if + (n_faces < 4, + "Error: VTK polyhedron cell " << cell_id << " has only " << n_faces << + " faces, but a polyhedron requires at least 4."); + + std::vector> sides(cast_int(n_faces)); + for (std::size_t s = 0; s != sides.size(); ++s) + { + const vtkIdType n_face_pts = face_stream->GetId(pos++); + auto side = std::make_shared(cast_int(n_face_pts)); + + for (vtkIdType n = 0; n != n_face_pts; ++n) + { + const vtkIdType vtk_point_id = face_stream->GetId(pos++); + const dof_id_type libmesh_node_id = node_id ? + vtk_node_to_libmesh[vtk_point_id] : + cast_int(vtk_point_id); + side->set_node(cast_int(n), + mesh.node_ptr(libmesh_node_id)); + } + + sides[s] = std::move(side); + } + + // Constructing the C0Polyhedron may create an interior "mid-element" + // node for its default tetrahedralization; if so, it is our job to + // add it to the mesh. The polyhedron's vertex node pointers were + // already assigned above via its polygonal sides. + std::unique_ptr mid_elem_node; + auto elem = std::make_unique(sides, mid_elem_node); + if (mid_elem_node) + mesh.add_node(std::move(mid_elem_node)); + + return elem; +} + +} // anonymous namespace + + + void VTKIO::read (const std::string & name) { // This is a serial-only process for now; @@ -194,18 +266,53 @@ void VTKIO::read (const std::string & name) elems_of_dimension.clear(); elems_of_dimension.resize(4, false); - // Use a typedef, because these names are just crazy - typedef vtkSmartPointer MyReader; - MyReader reader = MyReader::New(); - - // Pass the filename along to the reader - reader->SetFileName(name.c_str()); - - // Force reading - reader->Update(); + // VTK XML unstructured grid files come in two root-element forms: a + // plain serial "UnstructuredGrid" file (conventionally .vtu), and a + // "PUnstructuredGrid" descriptor (conventionally .pvtu) referencing one + // or more serial piece files for parallel use. Each needs a different + // reader class. Sniff the file's actual root + // attribute rather than assuming one from the extension (or, worse, + // always assuming one -- vtkXMLPUnstructuredGridReader silently yields + // an empty grid, rather than any diagnostic, when handed the other + // format). + std::string file_type; + { + vtkSmartPointer tester = + vtkSmartPointer::New(); + tester->SetFileName(name.c_str()); + + libmesh_error_msg_if + (!tester->TestReadFile(), + "Error: " << name << " does not appear to be a VTK XML file."); + + const char * file_type_cstr = tester->GetFileDataType(); + libmesh_error_msg_if + (!file_type_cstr, + "Error: " << name << " does not specify a VTK XML data type."); + file_type = file_type_cstr; + } - // read in the grid - _vtk_grid = reader->GetOutput(); + if (file_type == "PUnstructuredGrid") + { + vtkSmartPointer reader = + vtkSmartPointer::New(); + reader->SetFileName(name.c_str()); + reader->Update(); + _vtk_grid = reader->GetOutput(); + } + else if (file_type == "UnstructuredGrid") + { + vtkSmartPointer reader = + vtkSmartPointer::New(); + reader->SetFileName(name.c_str()); + reader->Update(); + _vtk_grid = reader->GetOutput(); + } + else + libmesh_error_msg + ("Error: " << name << " is a VTK XML file of type \"" << file_type + << "\", but libMesh's VTKIO reader only supports \"UnstructuredGrid\" " + "(.vtu) and \"PUnstructuredGrid\" (.pvtu) files."); // Get a reference to the mesh MeshBase & mesh = MeshInput::mesh(); @@ -286,30 +393,43 @@ void VTKIO::read (const std::string & name) { _vtk_grid->GetCell(i, cell); - // Get the libMesh element type corresponding to this VTK element type. - ElemType libmesh_elem_type = element_map.find(cell->GetCellType()); - auto elem = Elem::build(libmesh_elem_type); + std::unique_ptr elem; - // get the straightforward numbering from the VTK cells - for (auto j : elem->node_index_range()) + // VTK polyhedra are described by a stream of polygonal faces + // rather than by an ordered node list, so they are built + // separately as C0Polyhedron elements with their nodes already + // assigned via their sides. + if (cell->GetCellType() == VTK_POLYHEDRON) + elem = add_vtk_polyhedron(*_vtk_grid, mesh, + static_cast(i), + node_id, vtk_node_to_libmesh); + else { - const auto vtk_point_id = cell->GetPointId(j); - const dof_id_type libmesh_node_id = node_id ? - vtk_node_to_libmesh[vtk_point_id] : vtk_point_id; + // Get the libMesh element type corresponding to this VTK element type. + ElemType libmesh_elem_type = element_map.find(cell->GetCellType()); + elem = Elem::build(libmesh_elem_type); - elem->set_node(j, mesh.node_ptr(libmesh_node_id)); - } + // get the straightforward numbering from the VTK cells + for (auto j : elem->node_index_range()) + { + const auto vtk_point_id = cell->GetPointId(j); + const dof_id_type libmesh_node_id = node_id ? + vtk_node_to_libmesh[vtk_point_id] : vtk_point_id; - // then get the connectivity - std::vector conn; - elem->connectivity(0, VTK, conn); + elem->set_node(j, mesh.node_ptr(libmesh_node_id)); + } - // then reshuffle the nodes according to the connectivity, this - // two-time-assign would evade the definition of the vtk_mapping - for (unsigned int j=0, - n_conn = cast_int(conn.size()); - j != n_conn; ++j) - elem->set_node(j, mesh.node_ptr(conn[j])); + // then get the connectivity + std::vector conn; + elem->connectivity(0, VTK, conn); + + // then reshuffle the nodes according to the connectivity, this + // two-time-assign would evade the definition of the vtk_mapping + for (unsigned int j=0, + n_conn = cast_int(conn.size()); + j != n_conn; ++j) + elem->set_node(j, mesh.node_ptr(conn[j])); + } if (elem_id) { diff --git a/tests/Makefile.am b/tests/Makefile.am index 4abebf8114..d37b168fa6 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -230,6 +230,8 @@ data = matrices/geom_1_extraction_op.h5 \ meshes/quad4_tri3_smoothed.xda.gz \ meshes/quad4_tri3_hourglass.xda.gz \ meshes/hex8_prism6_smoothed.xda.gz \ + meshes/hex_prism_polyhedron.pvtu \ + meshes/hex_prism_polyhedron_0.vtu \ solutions/lagrange_vec_solution_mesh.xda \ solutions/lagrange_vec_solution.xda \ solutions/nedelec_one_solution_mesh.xda \ diff --git a/tests/mesh/mesh_input.C b/tests/mesh/mesh_input.C index 8dd2f8cbe0..59d8f6ccf2 100644 --- a/tests/mesh/mesh_input.C +++ b/tests/mesh/mesh_input.C @@ -9,6 +9,7 @@ #include #include #include +#include #include #include @@ -104,6 +105,9 @@ public: #ifdef LIBMESH_HAVE_VTK CPPUNIT_TEST( testVTKPreserveElemIds ); CPPUNIT_TEST( testVTKPreserveSubdomainIds ); + CPPUNIT_TEST( testVTKReadPolyhedra ); + CPPUNIT_TEST( testVTKReadPolyhedraSerial ); + CPPUNIT_TEST( testVTKReadNotAnXMLFile ); #endif #ifdef LIBMESH_HAVE_EXODUS_API @@ -331,6 +335,125 @@ public: } } } + + // Shared checks for a mesh that should contain exactly the hexagonal + // prism polyhedron (12 vertices, 8 faces -- two hexagons and six + // quads) written to meshes/hex_prism_polyhedron.{pvtu,_0.vtu}. + void checkHexPrismPolyhedron (MeshBase & mesh) + { + CPPUNIT_ASSERT_EQUAL(dof_id_type(1), mesh.n_elem()); + + // The element may live on a single processor when the mesh is + // distributed, so guard the checks accordingly. + const Elem * elem = mesh.query_elem_ptr(0); + bool found_elem = elem; + mesh.comm().max(found_elem); + CPPUNIT_ASSERT(found_elem); + + if (!elem) + return; + + CPPUNIT_ASSERT_EQUAL(C0POLYHEDRON, elem->type()); + CPPUNIT_ASSERT_EQUAL(12u, elem->n_vertices()); + CPPUNIT_ASSERT_EQUAL(8u, elem->n_sides()); + + // The file has no libmesh_node_id array, so libMesh node ids match + // the VTK point ordering, i.e. the coordinates we wrote out. + LIBMESH_ASSERT_FP_EQUAL(6.0, elem->volume(), TOLERANCE); + + // Collect the faces (as node-id sets) and compare against the faces + // that were written to the file, independent of any reordering the + // C0Polyhedron may do to its sides. + std::set> faces; + for (auto s : elem->side_index_range()) + { + std::set face; + for (const auto n : elem->nodes_on_side(s)) + face.insert(elem->node_id(n)); + faces.insert(std::move(face)); + } + + const std::set> expected_faces = + { {0, 1, 2, 3, 4, 5}, + {0, 1, 7, 6}, + {1, 2, 8, 7}, + {2, 3, 9, 8}, + {3, 4, 10, 9}, + {4, 5, 11, 10}, + {5, 0, 6, 11}, + {6, 7, 8, 9, 10, 11} }; + + CPPUNIT_ASSERT(faces == expected_faces); + } + + void testVTKReadPolyhedra () + { + LOG_UNIT_TEST; + + // This .pvtu (+ piece) file contains a single VTK_POLYHEDRON cell. + // VTKIO::read() sniffs the file's actual root element to tell a + // parallel PUnstructuredGrid (.pvtu) descriptor from a plain serial + // UnstructuredGrid (.vtu) file and uses the matching reader class -- + // see testVTKReadPolyhedraSerial below for the latter. + Mesh mesh(*TestCommWorld); + // Without this, prepare_for_use()'s default renumbering can reassign + // node ids (by local element-traversal order rather than preserving + // the ids we read), which would invalidate the "ids match VTK point + // ordering" assumption checkHexPrismPolyhedron() relies on. + mesh.allow_renumbering(false); + mesh.read("meshes/hex_prism_polyhedron.pvtu"); + mesh.prepare_for_use(); + + checkHexPrismPolyhedron(mesh); + } + + // The same polyhedron, read directly from the underlying serial + // UnstructuredGrid piece file rather than its .pvtu descriptor. This + // is a perfectly valid, common thing to hand a mesh reader -- many + // external tools write only this form, with no parallel wrapper -- but + // VTKIO::read() briefly lost the ability to read it at all: a Jan 2024 + // change to add .pvtu support replaced the serial reader class instead + // of adding the parallel one alongside it, silently regressing plain + // .vtu files to fail (which is what caused this CI failure in the + // first place). This is the regression test for that. + void testVTKReadPolyhedraSerial () + { + LOG_UNIT_TEST; + + Mesh mesh(*TestCommWorld); + mesh.allow_renumbering(false); + mesh.read("meshes/hex_prism_polyhedron_0.vtu"); + mesh.prepare_for_use(); + + checkHexPrismPolyhedron(mesh); + } + + // A file that isn't a VTK XML file at all (wrong format entirely, or a + // corrupted/truncated download) should be rejected clearly rather than + // being handed to vtkXMLFileReadTester/the VTK XML readers and failing + // in some less legible way. We reuse an existing, unrelated mesh + // fixture here rather than adding a new one; any non-VTK-XML file + // demonstrates this code path. + // + // We construct VTKIO directly here rather than going through + // mesh.read(), which would dispatch via NameBasedIO's rank-0-reads/ + // then-broadcast pattern: throwing out of that on rank 0 before the + // broadcast would leave every other rank blocked on a broadcast that + // never arrives. Every rank hits this same, rank-independent error + // when calling VTKIO::read() directly, so there's no such hazard here. + void testVTKReadNotAnXMLFile () + { + LOG_UNIT_TEST; + +#ifdef LIBMESH_ENABLE_EXCEPTIONS + Mesh mesh(*TestCommWorld); + VTKIO vtk(mesh); + CPPUNIT_ASSERT_THROW_MESSAGE + ("Non-VTK-XML file not detected", + vtk.read("meshes/circle.msh"), + libMesh::LogicError); +#endif + } #endif // LIBMESH_HAVE_VTK diff --git a/tests/meshes/hex_prism_polyhedron.pvtu b/tests/meshes/hex_prism_polyhedron.pvtu new file mode 100644 index 0000000000..cc88c97c1a --- /dev/null +++ b/tests/meshes/hex_prism_polyhedron.pvtu @@ -0,0 +1,9 @@ + + + + + + + + + diff --git a/tests/meshes/hex_prism_polyhedron_0.vtu b/tests/meshes/hex_prism_polyhedron_0.vtu new file mode 100644 index 0000000000..a3d82e21d2 --- /dev/null +++ b/tests/meshes/hex_prism_polyhedron_0.vtu @@ -0,0 +1,38 @@ + + + + + + + 0 -2 0 -1 -1 0 -1 1 0 0 2 0 1 1 0 1 -1 0 + 0 -2 1 -1 -1 1 -1 1 1 0 2 1 1 1 1 1 -1 1 + + + + + 0 1 2 3 4 5 6 7 8 9 10 11 + + + 12 + + + 42 + + + 8 + 6 0 1 2 3 4 5 + 4 0 1 7 6 + 4 1 2 8 7 + 4 2 3 9 8 + 4 3 4 10 9 + 4 4 5 11 10 + 4 5 0 6 11 + 6 6 7 8 9 10 11 + + + 45 + + + + +