From efc2d9073ccdb09123ce49a1299d2470f49532b7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Tue, 25 Mar 2025 19:34:00 +0100
Subject: [PATCH] Add VTK output for polynomial meshes of degree > 2

---
 src/output/VTKWriter.cpp | 113 +++++++++++++++++++++++++++++++++------
 1 file changed, 98 insertions(+), 15 deletions(-)

diff --git a/src/output/VTKWriter.cpp b/src/output/VTKWriter.cpp
index 19a613610..88922a65d 100644
--- a/src/output/VTKWriter.cpp
+++ b/src/output/VTKWriter.cpp
@@ -477,7 +477,8 @@ VTKWriter::_write(const MeshType& mesh,
     if constexpr (is_polygonal_mesh_v<MeshType>) {
       fout << mesh.numberOfNodes();
     } else if constexpr (is_polynomial_mesh_v<MeshType>) {
-      fout << mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces();
+      fout << mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
+                (mesh.degree() > 2) * mesh.numberOfCells();
     } else {
       throw NotImplementedError("unexpected mesh type");
     }
@@ -521,7 +522,26 @@ VTKWriter::_write(const MeshType& mesh,
         item_value_variant);
     }
     if constexpr (is_polynomial_mesh_v<MeshType>) {
-      Array<int> node_number(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces());
+      size_t number_of_inner_nodes = 0;
+      if (mesh.degree() > 2) {
+        const auto cell_type = mesh.connectivity().cellType();
+        for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
+          switch (cell_type[cell_id]) {
+          case CellType::Triangle: {
+            number_of_inner_nodes += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
+            break;
+          }
+          case CellType::Quadrangle: {
+            number_of_inner_nodes += (mesh.degree() - 2) * (mesh.degree() - 2);
+            break;
+          }
+          default: {
+            throw NormalError("invalid cell type for polynomal mesh of degree > 2");
+          }
+          }
+        }
+      }
+      Array<int> node_number(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() + number_of_inner_nodes);
       node_number.fill(0);
       _write_array(fout, "node_number", node_number, serialize_data_list);
     }
@@ -545,7 +565,8 @@ VTKWriter::_write(const MeshType& mesh,
     } else if constexpr (is_polynomial_mesh_v<MeshType>) {
       const NodeValue<const Rd>& xr = mesh.xr();
       const FaceArray<const Rd>& xl = mesh.xl();
-      Array<TinyVector<3>> positions(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces());
+      Array<TinyVector<3>> positions(mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces() +
+                                     (mesh.degree() > 2) * mesh.numberOfCells());
       parallel_for(
         mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
           for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
@@ -570,6 +591,25 @@ VTKWriter::_write(const MeshType& mesh,
             }
           }
         });
+
+      if (mesh.degree() > 2) {
+        const size_t cell_begining = mesh.numberOfNodes() + (mesh.degree() - 1) * mesh.numberOfFaces();
+
+        auto xj = MeshDataManager::instance().getMeshData(mesh).xj();
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
+            for (size_t i_node = 0; i_node < number_of_node_per_face; ++i_node) {
+              const size_t index = cell_begining + cell_id;
+              for (unsigned short i = 0; i < MeshType::Dimension; ++i) {
+                positions[index][i] = xj[cell_id][i];
+              }
+              for (unsigned short i = MeshType::Dimension; i < 3; ++i) {
+                positions[index][i] = 0;
+              }
+            }
+          });
+      }
+
       _write_array(fout, "Positions", positions, serialize_data_list);
     } else {
       throw NotImplementedError("unexpected mesh type");
@@ -596,25 +636,31 @@ VTKWriter::_write(const MeshType& mesh,
         _write_array(fout, "offsets", offsets, serialize_data_list);
       }
     } else if constexpr (is_polynomial_mesh_v<MeshType>) {
-      if (mesh.degree() > 2) {
-        throw NotImplementedError("only parabolic meshes are supported");
-      }
       const size_t degree = mesh.degree();
 
       const size_t number_of_vertices = mesh.numberOfNodes();
       const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
+      const auto& cell_type = mesh.connectivity().cellType();
+
       Array<unsigned int> offsets(mesh.numberOfCells());
       unsigned int offset = 0;
       for (CellId cell_id = 0; cell_id < mesh.numberOfCells(); ++cell_id) {
         const auto& cell_nodes = cell_to_node_matrix[cell_id];
         offset += degree * cell_nodes.size();
+        if (degree > 2) {
+          if (cell_type[cell_id] == CellType::Triangle) {
+            offset += (mesh.degree() - 1) * (mesh.degree() - 2) / 2;
+          } else if (cell_type[cell_id] == CellType::Quadrangle) {
+            offset += (mesh.degree() - 1) * (mesh.degree() - 1);
+          }
+        }
         offsets[cell_id] = offset;
       }
       _write_array(fout, "offsets", offsets, serialize_data_list);
 
-      const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
-      const auto& cell_type           = mesh.connectivity().cellType();
+      const auto& cell_to_face_matrix          = mesh.connectivity().cellToFaceMatrix();
+      const auto& cell_face_is_reversed_matrix = mesh.connectivity().cellFaceIsReversed();
 
       Array<NodeId::base_type> nodes{offset};
       size_t index = 0;
@@ -624,14 +670,51 @@ VTKWriter::_write(const MeshType& mesh,
         for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
           nodes[index++] = cell_nodes[i_node];
         }
-        if (degree == 2) {
+        if (degree > 1) {
+          auto cell_face_is_reversed = cell_face_is_reversed_matrix[cell_id];
           if (cell_type[cell_id] == CellType::Triangle) {
             for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-              nodes[index++] = number_of_vertices + cell_faces[(i_face + 2) % 3];
+              const size_t i_face_local_number = (i_face + 2) % 3;
+              if (cell_face_is_reversed[i_face_local_number]) {
+                for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
+                  nodes[index++] =
+                    number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
+                }
+              } else {
+                for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
+                  nodes[index++] =
+                    number_of_vertices + (degree - 1) * cell_faces[i_face_local_number] + i_face_inner_node;
+                }
+              }
+            }
+            if (degree > 2) {
+              for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 2) / 2; ++i) {
+                nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
+              }
             }
           } else {
             for (size_t i_face = 0; i_face < cell_faces.size(); ++i_face) {
-              nodes[index++] = number_of_vertices + cell_faces[i_face];
+              if ((i_face >= 2) xor cell_face_is_reversed[i_face]) {
+                for (ssize_t i_face_inner_node = degree - 2; i_face_inner_node >= 0; --i_face_inner_node) {
+                  nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
+                }
+              } else {
+                for (size_t i_face_inner_node = 0; i_face_inner_node < degree - 1; ++i_face_inner_node) {
+                  nodes[index++] = number_of_vertices + (degree - 1) * cell_faces[i_face] + i_face_inner_node;
+                }
+              }
+            }
+            if (degree > 2) {
+              if (cell_type[cell_id] != CellType::Quadrangle) {
+                std::ostringstream error_msg;
+                error_msg << "cell type: " << name(cell_type[cell_id])
+                          << " is not supported for VTK output of polynomial meshes of degree " << mesh.degree()
+                          << '\n';
+                throw NormalError(error_msg.str());
+              }
+              for (size_t i = 0; i < (mesh.degree() - 1) * (mesh.degree() - 1); ++i) {
+                nodes[index++] = number_of_vertices + (degree - 1) * mesh.numberOfFaces() + cell_id;
+              }
             }
           }
         }
@@ -673,7 +756,7 @@ VTKWriter::_write(const MeshType& mesh,
               } else if (element_degree == 2) {
                 types[j] = 22;
               } else {
-                throw NotImplementedError("unexpected degree for mesh writer");
+                types[j] = 69;
               }
             }
             break;
@@ -687,7 +770,7 @@ VTKWriter::_write(const MeshType& mesh,
               } else if (element_degree == 2) {
                 types[j] = 23;
               } else {
-                throw NotImplementedError("unexpected degree for mesh writer");
+                types[j] = 70;
               }
             }
             break;
@@ -696,12 +779,12 @@ VTKWriter::_write(const MeshType& mesh,
             if constexpr (is_polygonal_mesh_v<MeshType>) {
               types[j] = 7;
             } else if constexpr (is_polynomial_mesh_v<MeshType>) {
-              if (element_degree == 2) {
+              if (element_degree == 1) {
                 types[j] = 7;
               } else if (element_degree == 2) {
                 types[j] = 36;
               } else {
-                throw NotImplementedError("unexpected degree for mesh writer");
+                throw NormalError("output for polynomial cells of degree > 2 is not available for polygonal cells");
               }
             }
             break;
-- 
GitLab