Begin VTK output for general polyhedral meshes

src/output/VTKWriter.hpp

@@ -8,12 +8,14 @@
 #include <output/OutputNamedItemValueSet.hpp>
 #include <utils/Exceptions.hpp>
 #include <utils/Messenger.hpp>
+#include <utils/PugsAssert.hpp>
 
 #include <fstream>
 #include <iomanip>
 #include <iostream>
 #include <sstream>
 #include <string>
+#include <unordered_map>
 
 class VTKWriter
 {
@@ -234,6 +236,10 @@ class VTKWriter
         std::visit([&, name = name](auto&& item_value) { return this->_write_cell_pvtu(fout, name, item_value); },
                    item_value_variant);
       }
+      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";
 
       fout << "<PPointData>\n";
@@ -347,6 +353,10 @@ class VTKWriter
               types[j] = 12;
               break;
             }
+            case CellType::Diamond: {
+              types[j] = 42;
+              break;
+            }
             default: {
               std::ostringstream os;
               os << __FILE__ << ':' << __LINE__ << ": unknown cell type";
@@ -356,7 +366,62 @@ class VTKWriter
           });
         _write_array(fout, "types", types);
       }
+      if constexpr (MeshType::Dimension == 3) {
+        const auto& cell_to_face_matrix   = mesh->connectivity().cellToFaceMatrix();
+        const auto& cell_to_node_matrix   = mesh->connectivity().cellToNodeMatrix();
+        const auto& face_to_node_matrix   = mesh->connectivity().faceToNodeMatrix();
+        const auto& cell_face_is_reversed = mesh->connectivity().cellFaceIsReversed();
+
+        Array<size_t> faces_offsets(mesh->numberOfCells());
+        size_t next_offset = 0;
+        fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faces\">\n";
+        for (CellId cell_id = 0; cell_id < mesh->numberOfCells(); ++cell_id) {
+          const auto& cell_nodes = cell_to_node_matrix[cell_id];
+          std::unordered_map<size_t, size_t> node_id_to_local_number_in_cell;
+          for (size_t i_cell_node = 0; i_cell_node < cell_nodes.size(); ++i_cell_node) {
+            node_id_to_local_number_in_cell[cell_nodes[i_cell_node]] = i_cell_node;
+          }
+
+          const auto& cell_faces = cell_to_face_matrix[cell_id];
+          fout << cell_faces.size() << '\n';
+          next_offset++;
+          for (size_t i_cell_face = 0; i_cell_face < cell_faces.size(); ++i_cell_face) {
+            const FaceId& face_id  = cell_faces[i_cell_face];
+            const auto& face_nodes = face_to_node_matrix[face_id];
+            fout << face_nodes.size();
+            next_offset++;
+            Array<size_t> face_node_in_cell(face_nodes.size());
+
+            for (size_t i_face_node = 0; i_face_node < face_nodes.size(); ++i_face_node) {
+              const NodeId& node_id                  = face_nodes[i_face_node];
+              auto i_node_id_to_local_number_in_cell = node_id_to_local_number_in_cell.find(node_id);
+              Assert(i_node_id_to_local_number_in_cell != node_id_to_local_number_in_cell.end());
+              face_node_in_cell[i_face_node] = i_node_id_to_local_number_in_cell->second;
+            }
+
+            if (cell_face_is_reversed(cell_id, i_cell_face)) {
+              for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+                fout << ' ' << face_node_in_cell[face_node_in_cell.size() - 1 - i];
+              }
+            } else {
+              for (size_t i = 0; i < face_node_in_cell.size(); ++i) {
+                fout << ' ' << face_node_in_cell[i];
+              }
+            }
+
+            next_offset += face_nodes.size();
 
+            fout << '\n';
+          }
+          faces_offsets[cell_id] = next_offset;
+        }
+        fout << "</DataArray>\n";
+        fout << "<DataArray type=\"Int64\" IdType=\"1\" Name=\"faceoffsets\">\n";
+        for (size_t i_face_offsets = 0; i_face_offsets < faces_offsets.size(); ++i_face_offsets) {
+          fout << faces_offsets[i_face_offsets] << '\n';
+        }
+        fout << "</DataArray>\n";
+      }
       fout << "</Cells>\n";
       fout << "</Piece>\n";
       fout << "</UnstructuredGrid>\n";