DataOut for pyramids: reorder at the end

dealii · Oct 19, 2021 · c083c45 · c083c45
1 parent 6e182af
commit c083c45
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 16 deletions.
diff --git 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
@@ -781,21 +781,28 @@ namespace
   {
     Point<spacedim> node;
 
+    unsigned int point_no_actual = point_no;
+
+    if (patch.reference_cell == ReferenceCells::Pyramid)
+      {
+        AssertDimension(patch.n_subdivisions, 1);
+
+        static 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;