Read VTK_POLYHEDRON cells in the VTU reader

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 <cstring>
 #include <fstream>
 
 
@@ -183,6 +189,72 @@ std::map<ElemMappingType, VTKIO::ElementMaps> 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<Elem>
+add_vtk_polyhedron(vtkUnstructuredGrid & vtk_grid,
+                   MeshBase & mesh,
+                   const vtkIdType cell_id,
+                   vtkIntArray * node_id,
+                   const std::vector<dof_id_type> & 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<vtkIdList> face_stream = vtkSmartPointer<vtkIdList>::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<std::shared_ptr<Polygon>> sides(cast_int<std::size_t>(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<C0Polygon>(cast_int<unsigned 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<dof_id_type>(vtk_point_id);
+          side->set_node(cast_int<unsigned 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<Node> mid_elem_node;
+  auto elem = std::make_unique<C0Polyhedron>(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<vtkXMLPUnstructuredGridReader> 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 <VTKFile type="..."/>
+  // 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<vtkXMLFileReadTester> tester =
+      vtkSmartPointer<vtkXMLFileReadTester>::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<vtkXMLPUnstructuredGridReader> reader =
+        vtkSmartPointer<vtkXMLPUnstructuredGridReader>::New();
+      reader->SetFileName(name.c_str());
+      reader->Update();
+      _vtk_grid = reader->GetOutput();
+    }
+  else if (file_type == "UnstructuredGrid")
+    {
+      vtkSmartPointer<vtkXMLUnstructuredGridReader> reader =
+        vtkSmartPointer<vtkXMLUnstructuredGridReader>::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<MeshBase>::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> 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<vtkIdType>(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<dof_id_type> 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<unsigned int>(conn.size());
-           j != n_conn; ++j)
-        elem->set_node(j, mesh.node_ptr(conn[j]));
+          // then get the connectivity
+          std::vector<dof_id_type> 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<unsigned int>(conn.size());
+               j != n_conn; ++j)
+            elem->set_node(j, mesh.node_ptr(conn[j]));
+        }
 
       if (elem_id)
         {

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 \

tests/mesh/mesh_input.C

@@ -9,6 +9,7 @@
 #include <libmesh/replicated_mesh.h>
 #include <libmesh/enum_norm_type.h>
 #include <libmesh/enum_to_string.h>
+#include <libmesh/cell_c0polyhedron.h>
 
 #include <libmesh/abaqus_io.h>
 #include <libmesh/dyna_io.h>
@@ -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<std::set<dof_id_type>> faces;
+    for (auto s : elem->side_index_range())
+      {
+        std::set<dof_id_type> face;
+        for (const auto n : elem->nodes_on_side(s))
+          face.insert(elem->node_id(n));
+        faces.insert(std::move(face));
+      }
+
+    const std::set<std::set<dof_id_type>> 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
 
 

tests/meshes/hex_prism_polyhedron.pvtu

@@ -0,0 +1,9 @@
+<?xml version="1.0"?>
+<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian">
+  <PUnstructuredGrid GhostLevel="0">
+    <PPoints>
+      <PDataArray type="Float64" Name="Points" NumberOfComponents="3"/>
+    </PPoints>
+    <Piece Source="hex_prism_polyhedron_0.vtu"/>
+  </PUnstructuredGrid>
+</VTKFile>