From d4e0372fc1018875a98140ae335d03a77ab1cc31 Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Tue, 10 Jan 2023 18:42:08 -0700
Subject: [PATCH] Centralize translation between deal.II and VTK numbering for
 pyramids.

---
 include/deal.II/grid/reference_cell.h |  6 ++++
 source/base/data_out_base.cc          | 14 +++-----
 source/grid/reference_cell.cc         | 52 +++++++++++++++++++++++++++
 3 files changed, 62 insertions(+), 10 deletions(-)

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index ffa3244207..562a032b72 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -734,6 +734,12 @@ public:
     const std::array<unsigned, dim> &nodes_per_direction,
     const bool                       legacy_format) const;
 
+  /**
+   * Map a VTK vertex number to a deal.II vertex number.
+   */
+  unsigned int
+  vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const;
+
   /**
    * Return the GMSH element type code that corresponds to the reference cell.
    */
diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index 5813b8befa..d661106371 100644
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -2642,7 +2642,7 @@ namespace DataOutBase
                                   patch.reference_cell);
             first_vertex_of_patch += patch.data.n_cols();
           }
-        else
+        else // hypercube cell
           {
             const unsigned int n_subdivisions = patch.n_subdivisions;
             const unsigned int n              = n_subdivisions + 1;
@@ -5900,8 +5900,6 @@ namespace DataOutBase
               Assert(patch.n_subdivisions == 1, ExcNotImplemented());
 
               const unsigned int n_points = patch.data.n_cols();
-              static const std::array<unsigned int, 5>
-                pyramid_index_translation_table = {{0, 1, 3, 2, 4}};
 
               if (deal_ii_with_zlib &&
                   (flags.compression_level !=
@@ -5910,18 +5908,14 @@ namespace DataOutBase
                   for (unsigned int i = 0; i < n_points; ++i)
                     cells.push_back(
                       first_vertex_of_patch +
-                      (patch.reference_cell == ReferenceCells::Pyramid ?
-                         pyramid_index_translation_table[i] :
-                         i));
+                      patch.reference_cell.vtk_vertex_to_deal_vertex(i));
                 }
               else
                 {
                   for (unsigned int i = 0; i < n_points; ++i)
                     o << '\t'
-                      << first_vertex_of_patch +
-                           (patch.reference_cell == ReferenceCells::Pyramid ?
-                              pyramid_index_translation_table[i] :
-                              i);
+                      << (first_vertex_of_patch +
+                          patch.reference_cell.vtk_vertex_to_deal_vertex(i));
                   o << '\n';
                 }
 
diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc
index 2f69862f54..3ae63c4556 100644
--- a/source/grid/reference_cell.cc
+++ b/source/grid/reference_cell.cc
@@ -655,6 +655,58 @@ ReferenceCell::vtk_lexicographic_to_node_index<3>(
 
 
 
+unsigned int
+ReferenceCell::vtk_vertex_to_deal_vertex(const unsigned int vertex_index) const
+{
+  AssertIndexRange(vertex_index, n_vertices());
+
+  // For some of the following, deal.II uses the same ordering as VTK
+  // and in that case, we only need to return 'vertex_index' (i.e.,
+  // use the identity mapping). For some others, we need to translate.
+  //
+  // For the ordering, see the VTK manual (for example at
+  // http://www.princeton.edu/~efeibush/viscourse/vtk.pdf, page 9).
+  if (*this == ReferenceCells::Vertex)
+    return vertex_index;
+  else if (*this == ReferenceCells::Line)
+    return vertex_index;
+  else if (*this == ReferenceCells::Triangle)
+    return vertex_index;
+  else if (*this == ReferenceCells::Quadrilateral)
+    {
+      static constexpr std::array<unsigned int, 4> index_translation_table = {
+        {0, 1, 3, 2}};
+      return index_translation_table[vertex_index];
+    }
+  else if (*this == ReferenceCells::Tetrahedron)
+    return vertex_index;
+  else if (*this == ReferenceCells::Pyramid)
+    {
+      static constexpr std::array<unsigned int, 5> index_translation_table = {
+        {0, 1, 3, 2, 4}};
+      return index_translation_table[vertex_index];
+    }
+  else if (*this == ReferenceCells::Wedge)
+    return vertex_index;
+  else if (*this == ReferenceCells::Hexahedron)
+    {
+      static constexpr std::array<unsigned int, 8> index_translation_table = {
+        {0, 1, 3, 2, 4, 5, 7, 6}};
+      return index_translation_table[vertex_index];
+    }
+  else if (*this == ReferenceCells::Invalid)
+    {
+      Assert(false, ExcNotImplemented());
+      return numbers::invalid_unsigned_int;
+    }
+
+  Assert(false, ExcNotImplemented());
+
+  return numbers::invalid_unsigned_int;
+}
+
+
+
 unsigned int
 ReferenceCell::gmsh_element_type() const
 {
-- 
2.39.5