Only write faces and faceoffsets fields in pvtu when needed

Actually, it seems that reading of general polyhedral cells is broken
in early versions of Paraview
- is it a change in the format (poorly documented...)
- or just a Paraview's bug?

src/output/VTKWriter.cpp

@@ -368,80 +368,7 @@ VTKWriter::_write(const MeshType& mesh,
 
   createDirectoryIfNeeded(m_base_filename);
 
-  if (parallel::rank() == 0) {   // write PVTK file
-    std::ofstream fout(_getFilenamePVTU());
-    if (not fout) {
-      std::ostringstream error_msg;
-      error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenamePVTU() << rang::fg::reset << '"';
-      throw NormalError(error_msg.str());
-    }
-
-    fout << "<?xml version=\"1.0\"?>\n";
-    fout << _getDateAndVersionComment();
-    fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
-    fout << "<PUnstructuredGrid GhostLevel=\"0\">\n";
-
-    fout << "<PPoints>\n";
-    fout << "<PDataArray Name=\"Positions\" NumberOfComponents=\"3\" "
-            "type=\"Float64\"/>\n";
-    fout << "</PPoints>\n";
-
-    fout << "<PCells>\n";
-    fout << "<PDataArray type=\"Int32\" Name=\"connectivity\" "
-            "NumberOfComponents=\"1\"/>\n";
-    fout << "<PDataArray type=\"UInt32\" Name=\"offsets\" "
-            "NumberOfComponents=\"1\"/>\n";
-    fout << "<PDataArray type=\"Int8\" Name=\"types\" "
-            "NumberOfComponents=\"1\"/>\n";
-    if constexpr (MeshType::Dimension == 3) {
-      fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
-      fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
-    }
-    fout << "</PCells>\n";
-
-    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit(
-        [&, name = name](auto&& item_value) {
-          using ItemValueType = std::decay_t<decltype(item_value)>;
-          if constexpr (ItemValueType::item_t == ItemType::edge) {
-            std::ostringstream error_msg;
-            error_msg << "VTK format does not support edge data, cannot save variable \"" << rang::fgB::yellow << name
-                      << rang::fg::reset << '"';
-            throw NormalError(error_msg.str());
-          }
-          if constexpr (ItemValueType::item_t == ItemType::face) {
-            std::ostringstream error_msg;
-            error_msg << "VTK format does not support face data, cannot save variable \"" << rang::fgB::yellow << name
-                      << rang::fg::reset << '"';
-            throw NormalError(error_msg.str());
-          }
-        },
-        item_value_variant);
-    }
-
-    fout << "<PPointData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
-                 item_value_variant);
-    }
-    fout << "</PPointData>\n";
-
-    fout << "<PCellData>\n";
-    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
-      std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
-                 item_value_variant);
-    }
-    if (parallel::size() > 1) {
-      fout << "<PDataArray type=\"UInt8\" Name=\"vtkGhostType\" NumberOfComponents=\"1\"/>\n";
-    }
-    fout << "</PCellData>\n";
-
-    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
-      fout << "<Piece Source=" << std::filesystem::path{_getFilenameVTU(i_rank)}.filename() << "/>\n";
-    }
-    fout << "</PUnstructuredGrid>\n";
-    fout << "</VTKFile>\n";
-  }
+  bool has_general_polyhedron = false;
 
   {
     SerializedDataList serialize_data_list;
@@ -584,7 +511,7 @@ VTKWriter::_write(const MeshType& mesh,
         });
       _write_array(fout, "types", types, serialize_data_list);
       if constexpr (MeshType::Dimension == 3) {
-        const bool has_general_polyhedron = [&] {
+        has_general_polyhedron = [&] {
           for (size_t i = 0; i < types.size(); ++i) {
             if (types[i] == 41)
               return true;
@@ -660,6 +587,83 @@ VTKWriter::_write(const MeshType& mesh,
     fout << "</VTKFile>\n";
   }
 
+  if (parallel::rank() == 0) {   // write PVTK file
+    std::ofstream fout(_getFilenamePVTU());
+    if (not fout) {
+      std::ostringstream error_msg;
+      error_msg << "cannot create file \"" << rang::fgB::yellow << _getFilenamePVTU() << rang::fg::reset << '"';
+      throw NormalError(error_msg.str());
+    }
+
+    fout << "<?xml version=\"1.0\"?>\n";
+    fout << _getDateAndVersionComment();
+    fout << "<VTKFile type=\"PUnstructuredGrid\">\n";
+    fout << "<PUnstructuredGrid GhostLevel=\"0\">\n";
+
+    fout << "<PPoints>\n";
+    fout << "<PDataArray Name=\"Positions\" NumberOfComponents=\"3\" "
+            "type=\"Float64\"/>\n";
+    fout << "</PPoints>\n";
+
+    fout << "<PCells>\n";
+    fout << "<PDataArray type=\"Int32\" Name=\"connectivity\" "
+            "NumberOfComponents=\"1\"/>\n";
+    fout << "<PDataArray type=\"UInt32\" Name=\"offsets\" "
+            "NumberOfComponents=\"1\"/>\n";
+    fout << "<PDataArray type=\"Int8\" Name=\"types\" "
+            "NumberOfComponents=\"1\"/>\n";
+    if constexpr (MeshType::Dimension == 3) {
+      if (has_general_polyhedron) {
+        fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faces\"/>\n";
+        fout << "<PDataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\"/>\n";
+      }
+    }
+    fout << "</PCells>\n";
+
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
+      std::visit(
+        [&, name = name](auto&& item_value) {
+          using ItemValueType = std::decay_t<decltype(item_value)>;
+          if constexpr (ItemValueType::item_t == ItemType::edge) {
+            std::ostringstream error_msg;
+            error_msg << "VTK format does not support edge data, cannot save variable \"" << rang::fgB::yellow << name
+                      << rang::fg::reset << '"';
+            throw NormalError(error_msg.str());
+          }
+          if constexpr (ItemValueType::item_t == ItemType::face) {
+            std::ostringstream error_msg;
+            error_msg << "VTK format does not support face data, cannot save variable \"" << rang::fgB::yellow << name
+                      << rang::fg::reset << '"';
+            throw NormalError(error_msg.str());
+          }
+        },
+        item_value_variant);
+    }
+
+    fout << "<PPointData>\n";
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
+      std::visit([&, name = name](auto&& item_value) { return this->_write_node_pvtu(fout, name, item_value); },
+                 item_value_variant);
+    }
+    fout << "</PPointData>\n";
+
+    fout << "<PCellData>\n";
+    for (const auto& [name, item_value_variant] : output_named_item_data_set) {
+      std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
+                 item_value_variant);
+    }
+    if (parallel::size() > 1) {
+      fout << "<PDataArray type=\"UInt8\" Name=\"vtkGhostType\" NumberOfComponents=\"1\"/>\n";
+    }
+    fout << "</PCellData>\n";
+
+    for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+      fout << "<Piece Source=" << std::filesystem::path{_getFilenameVTU(i_rank)}.filename() << "/>\n";
+    }
+    fout << "</PUnstructuredGrid>\n";
+    fout << "</VTKFile>\n";
+  }
+
   if (parallel::rank() == 0) {   // write PVD file
     const std::string pvd_filename = m_base_filename + ".pvd";
     std::ofstream fout(pvd_filename);