From a8ae801e52ea10e513e50d1ce4f03ca2a59c888d Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Tue, 19 Oct 2021 09:52:11 +0200
Subject: [PATCH] DataOut for pyramids: reorder at the end

---
 include/deal.II/base/data_out_base.h          |  7 +-
 .../numerics/data_out_dof_data.templates.h    | 14 +--
 source/base/data_out_base.cc                  | 99 +++++++++++++------
 tests/simplex/data_out_write_vtk_04.output    | 50 +++++-----
 4 files changed, 100 insertions(+), 70 deletions(-)

diff --git a/include/deal.II/base/data_out_base.h b/include/deal.II/base/data_out_base.h
index 528f7fc190..4a049fdbae 100644
--- a/include/deal.II/base/data_out_base.h
+++ b/include/deal.II/base/data_out_base.h
@@ -1393,9 +1393,10 @@ namespace DataOutBase
      * internal reordered format.
      */
     void
-    write_cell_single(const unsigned int index,
-                      const unsigned int start,
-                      const unsigned int n_points);
+    write_cell_single(const unsigned int   index,
+                      const unsigned int   start,
+                      const unsigned int   n_points,
+                      const ReferenceCell &reference_cell);
 
     /**
      * Filter and record a data set. If there are multiple values at a given
diff --git a/include/deal.II/numerics/data_out_dof_data.templates.h b/include/deal.II/numerics/data_out_dof_data.templates.h
index 5a3fdaa67a..763083b353 100644
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -320,16 +320,10 @@ namespace internal
           if (needs_pyramid_setup)
             {
               Assert(n_subdivisions == 1, ExcNotImplemented());
-
-              std::vector<Point<dim>> points;
-
-              points.emplace_back(-1.0, -1.0, 0.0);
-              points.emplace_back(+1.0, -1.0, 0.0);
-              points.emplace_back(+1.0, +1.0, 0.0);
-              points.emplace_back(-1.0, +1.0, 0.0);
-              points.emplace_back(+0.0, +0.0, 1.0);
-
-              quadrature_pyramid = std::make_unique<Quadrature<dim>>(points);
+              quadrature_pyramid = std::make_unique<Quadrature<dim>>(
+                FE_PyramidP<dim, spacedim>(
+                  1 /*note: vtk only supports linear wedges*/)
+                  .get_unit_support_points());
             }
 
           n_q_points =
diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index 3121729503..4a66fadcd3 100644
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -496,18 +496,22 @@ namespace DataOutBase
 
 
   void
-  DataOutFilter::write_cell_single(const unsigned int index,
-                                   const unsigned int start,
-                                   const unsigned int n_points)
+  DataOutFilter::write_cell_single(const unsigned int   index,
+                                   const unsigned int   start,
+                                   const unsigned int   n_points,
+                                   const ReferenceCell &reference_cell)
   {
     ++num_cells;
 
     const unsigned int base_entry = index * n_points;
 
+    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
+
     for (unsigned int i = 0; i < n_points; ++i)
-      {
-        internal_add_cell(base_entry + i, start + i);
-      }
+      internal_add_cell(base_entry + i,
+                        start + (reference_cell == ReferenceCells::Pyramid ?
+                                   table[i] :
+                                   i));
   }
 
 
@@ -781,21 +785,28 @@ namespace
   {
     Point<spacedim> node;
 
+    unsigned int point_no_actual = point_no;
+
+    if (patch.reference_cell == ReferenceCells::Pyramid)
+      {
+        AssertDimension(patch.n_subdivisions, 1);
+
+        static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
+        point_no_actual                                = table[point_no];
+      }
+
     if (patch.points_are_available)
       {
         for (unsigned int d = 0; d < spacedim; ++d)
-          node[d] = patch.data(patch.data.size(0) - spacedim + d, point_no);
+          node[d] =
+            patch.data(patch.data.size(0) - spacedim + d, point_no_actual);
         return node;
       }
     else
       {
         AssertDimension(patch.n_subdivisions, 1);
-        Assert(
-          patch.reference_cell != ReferenceCells::Pyramid,
-          ExcMessage(
-            "Pyramids need different ordering of the vertices, which is not implemented yet here."));
 
-        node = patch.vertices[point_no];
+        node = patch.vertices[point_no_actual];
       }
 
     return node;
@@ -1072,13 +1083,15 @@ namespace
      * @note All inheriting classes should implement this function.
      */
     void
-    write_cell_single(const unsigned int index,
-                      const unsigned int start,
-                      const unsigned int n_points)
+    write_cell_single(const unsigned int   index,
+                      const unsigned int   start,
+                      const unsigned int   n_points,
+                      const ReferenceCell &reference_cell)
     {
       (void)index;
       (void)start;
       (void)n_points;
+      (void)reference_cell;
 
       Assert(false,
              ExcMessage("The derived class you are using needs to "
@@ -1299,9 +1312,10 @@ namespace
      * Print vertices [start, start+n_points[
      */
     void
-    write_cell_single(const unsigned int index,
-                      const unsigned int start,
-                      const unsigned int n_points);
+    write_cell_single(const unsigned int   index,
+                      const unsigned int   start,
+                      const unsigned int   n_points,
+                      const ReferenceCell &reference_cell);
 
     /**
      * Write a high-order cell type, i.e., a Lagrange cell
@@ -1350,9 +1364,10 @@ namespace
      * Print vertices [start, start+n_points[
      */
     void
-    write_cell_single(const unsigned int index,
-                      const unsigned int start,
-                      const unsigned int n_points);
+    write_cell_single(const unsigned int   index,
+                      const unsigned int   start,
+                      const unsigned int   n_points,
+                      const ReferenceCell &reference_cell);
 
     /**
      * Write a high-order cell type, i.e., a Lagrange cell
@@ -1749,15 +1764,20 @@ namespace
   }
 
   void
-  VtkStream::write_cell_single(const unsigned int index,
-                               const unsigned int start,
-                               const unsigned int n_points)
+  VtkStream::write_cell_single(const unsigned int   index,
+                               const unsigned int   start,
+                               const unsigned int   n_points,
+                               const ReferenceCell &reference_cell)
   {
     (void)index;
 
+    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
+
     stream << '\t' << n_points;
     for (unsigned int i = 0; i < n_points; ++i)
-      stream << '\t' << start + i;
+      stream << '\t'
+             << start +
+                  (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
     stream << '\n';
   }
 
@@ -1857,19 +1877,25 @@ namespace
   }
 
   void
-  VtuStream::write_cell_single(const unsigned int index,
-                               const unsigned int start,
-                               const unsigned int n_points)
+  VtuStream::write_cell_single(const unsigned int   index,
+                               const unsigned int   start,
+                               const unsigned int   n_points,
+                               const ReferenceCell &reference_cell)
   {
     (void)index;
 
+    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
+
 #if !defined(DEAL_II_WITH_ZLIB)
     for (unsigned int i = 0; i < n_points; ++i)
-      stream << '\t' << start + i;
+      stream << '\t'
+             << start +
+                  (reference_cell == ReferenceCells::Pyramid ? table[i] : i);
     stream << '\n';
 #else
     for (unsigned int i = 0; i < n_points; ++i)
-      cells.push_back(start + i);
+      cells.push_back(
+        start + (reference_cell == ReferenceCells::Pyramid ? table[i] : i));
 #endif
   }
 
@@ -2701,6 +2727,8 @@ namespace DataOutBase
     Assert(dim <= 3, ExcNotImplemented());
     unsigned int count = 0;
 
+    static const std::array<unsigned int, 5> table = {{0, 1, 3, 2, 4}};
+
     for (const auto &patch : patches)
       {
         // special treatment of non-hypercube cells
@@ -2708,7 +2736,13 @@ namespace DataOutBase
           {
             for (unsigned int point_no = 0; point_no < patch.data.n_cols();
                  ++point_no)
-              out.write_point(count++, compute_arbitrary_node(patch, point_no));
+              out.write_point(
+                count++,
+                compute_arbitrary_node(patch,
+                                       (patch.reference_cell ==
+                                            ReferenceCells::Pyramid ?
+                                          table[point_no] :
+                                          point_no)));
           }
         else
           {
@@ -2745,7 +2779,8 @@ namespace DataOutBase
           {
             out.write_cell_single(count++,
                                   first_vertex_of_patch,
-                                  patch.data.n_cols());
+                                  patch.data.n_cols(),
+                                  patch.reference_cell);
             first_vertex_of_patch += patch.data.n_cols();
           }
         else
diff --git a/tests/simplex/data_out_write_vtk_04.output b/tests/simplex/data_out_write_vtk_04.output
index f65b314b34..b10ff6ab57 100644
--- a/tests/simplex/data_out_write_vtk_04.output
+++ b/tests/simplex/data_out_write_vtk_04.output
@@ -62,49 +62,49 @@ DATASET UNSTRUCTURED_GRID
 POINTS 30 double
 0.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
-0.00000 1.00000 1.00000
 0.00000 0.00000 1.00000
+0.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-1.00000 1.00000 1.00000
 0.00000 1.00000 1.00000
+1.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 1.00000 1.00000
 1.00000 1.00000 1.00000
-1.00000 0.00000 1.00000
 0.00000 0.00000 1.00000
+1.00000 0.00000 1.00000
 0.500000 0.500000 0.500000
 1.00000 0.00000 0.00000
 1.00000 0.00000 1.00000
-1.00000 1.00000 1.00000
 1.00000 1.00000 0.00000
+1.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 0.00000 0.00000
 1.00000 0.00000 0.00000
-1.00000 1.00000 0.00000
 0.00000 1.00000 0.00000
+1.00000 1.00000 0.00000
 0.500000 0.500000 0.500000
 0.00000 0.00000 0.00000
 0.00000 0.00000 1.00000
-1.00000 0.00000 1.00000
 1.00000 0.00000 0.00000
+1.00000 0.00000 1.00000
 0.500000 0.500000 0.500000
 
 CELLS 6 36
-	5	0	1	2	3	4
-	5	5	6	7	8	9
-	5	10	11	12	13	14
-	5	15	16	17	18	19
-	5	20	21	22	23	24
-	5	25	26	27	28	29
+	5	0	1	3	2	4
+	5	5	6	8	7	9
+	5	10	11	13	12	14
+	5	15	16	18	17	19
+	5	20	21	23	22	24
+	5	25	26	28	27	29
 
 CELL_TYPES 6
  14 14 14 14 14 14
 POINT_DATA 30
 SCALARS solution double 1
 LOOKUP_TABLE default
-0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 1.00000 1.00000 0.00000 0.250000 0.00000 1.00000 1.00000 0.00000 0.250000 1.00000 1.00000 1.00000 1.00000 0.250000 0.00000 1.00000 1.00000 0.00000 0.250000 0.00000 0.00000 1.00000 1.00000 0.250000 
+0.00000 0.00000 0.00000 0.00000 0.250000 0.00000 1.00000 0.00000 1.00000 0.250000 0.00000 1.00000 0.00000 1.00000 0.250000 1.00000 1.00000 1.00000 1.00000 0.250000 0.00000 1.00000 0.00000 1.00000 0.250000 0.00000 0.00000 1.00000 1.00000 0.250000 
 # vtk DataFile Version 3.0
 #This file was generated 
 ASCII
@@ -113,42 +113,42 @@ DATASET UNSTRUCTURED_GRID
 POINTS 30 double
 0.00000 0.00000 0.00000
 0.00000 1.00000 0.00000
-0.00000 1.00000 1.00000
 0.00000 0.00000 1.00000
+0.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 1.00000 0.00000
 1.00000 1.00000 0.00000
-1.00000 1.00000 1.00000
 0.00000 1.00000 1.00000
+1.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 1.00000 1.00000
 1.00000 1.00000 1.00000
-1.00000 0.00000 1.00000
 0.00000 0.00000 1.00000
+1.00000 0.00000 1.00000
 0.500000 0.500000 0.500000
 1.00000 0.00000 0.00000
 1.00000 0.00000 1.00000
-1.00000 1.00000 1.00000
 1.00000 1.00000 0.00000
+1.00000 1.00000 1.00000
 0.500000 0.500000 0.500000
 0.00000 0.00000 0.00000
 1.00000 0.00000 0.00000
-1.00000 1.00000 0.00000
 0.00000 1.00000 0.00000
+1.00000 1.00000 0.00000
 0.500000 0.500000 0.500000
 0.00000 0.00000 0.00000
 0.00000 0.00000 1.00000
-1.00000 0.00000 1.00000
 1.00000 0.00000 0.00000
+1.00000 0.00000 1.00000
 0.500000 0.500000 0.500000
 
 CELLS 6 36
-	5	0	1	2	3	4
-	5	5	6	7	8	9
-	5	10	11	12	13	14
-	5	15	16	17	18	19
-	5	20	21	22	23	24
-	5	25	26	27	28	29
+	5	0	1	3	2	4
+	5	5	6	8	7	9
+	5	10	11	13	12	14
+	5	15	16	18	17	19
+	5	20	21	23	22	24
+	5	25	26	28	27	29
 
 CELL_TYPES 6
  14 14 14 14 14 14
-- 
2.39.5