From e82b008b3a02804efebee3f98a1d254c9224b22d Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Sat, 5 Dec 2020 12:51:53 +0100
Subject: [PATCH] Enable DataOut for Pyramdis

---
 .../numerics/data_out_dof_data.templates.h    | 36 ++++++++++++++++---
 source/base/data_out_base.cc                  |  6 ++++
 2 files changed, 37 insertions(+), 5 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
index e9f7ee2b61..efa7446471 100644
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -110,6 +110,7 @@ namespace internal
           bool needs_hypercube_setup = false;
           bool needs_simplex_setup   = false;
           bool needs_wedge_setup     = false;
+          bool needs_pyramid_setup   = false;
 
           for (const auto &fe : finite_elements)
             for (unsigned int i = 0; i < fe->size(); ++i)
@@ -128,6 +129,9 @@ namespace internal
                   case ReferenceCell::Type::Wedge:
                     needs_wedge_setup |= true;
                     break;
+                  case ReferenceCell::Type::Pyramid:
+                    needs_pyramid_setup |= true;
+                    break;
                   default:
                     Assert(false, ExcNotImplemented());
                 }
@@ -135,6 +139,7 @@ namespace internal
           std::unique_ptr<dealii::Quadrature<dim>> quadrature_simplex;
           std::unique_ptr<dealii::Quadrature<dim>> quadrature_hypercube;
           std::unique_ptr<dealii::Quadrature<dim>> quadrature_wedge;
+          std::unique_ptr<dealii::Quadrature<dim>> quadrature_pyramid;
 
           if (needs_simplex_setup)
             {
@@ -158,10 +163,27 @@ namespace internal
                   .get_unit_support_points());
             }
 
-          n_q_points = std::max(
-            {needs_wedge_setup ? quadrature_wedge->size() : 0,
-             needs_simplex_setup ? quadrature_simplex->size() : 0,
-             needs_hypercube_setup ? quadrature_hypercube->size() : 0});
+          if (needs_pyramid_setup)
+            {
+              Assert(1 <= n_subdivisions && n_subdivisions <= 2,
+                     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);
+            }
+
+          n_q_points =
+            std::max({needs_wedge_setup ? quadrature_wedge->size() : 0,
+                      needs_simplex_setup ? quadrature_simplex->size() : 0,
+                      needs_hypercube_setup ? quadrature_hypercube->size() : 0,
+                      needs_pyramid_setup ? quadrature_pyramid->size() : 0});
 
           x_fe_values.resize(finite_elements.size());
           for (unsigned int i = 0; i < finite_elements.size(); ++i)
@@ -194,6 +216,9 @@ namespace internal
                         case ReferenceCell::Type::Wedge:
                           quadrature.push_back(*quadrature_wedge);
                           break;
+                        case ReferenceCell::Type::Pyramid:
+                          quadrature.push_back(*quadrature_pyramid);
+                          break;
                         default:
                           Assert(false, ExcNotImplemented());
                       }
@@ -206,7 +231,8 @@ namespace internal
                       // certain reference-cell kinds take by default the
                       // curved code path in DataOut::build_one_patch, for
                       // which the quadrature points need to be enabled
-                      (needs_simplex_setup || needs_wedge_setup) ?
+                      (needs_simplex_setup || needs_wedge_setup ||
+                       needs_pyramid_setup) ?
                         (update_flags | update_quadrature_points) :
                         update_flags);
                 }
diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index 860afc305b..1f2e3c1387 100644
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -640,6 +640,12 @@ namespace
         vtk_cell_id[0] = 13;
         vtk_cell_id[1] = 1;
       }
+    else if (patch.reference_cell_type == ReferenceCell::Type::Pyramid &&
+             patch.data.n_cols() == 5)
+      {
+        vtk_cell_id[0] = 14;
+        vtk_cell_id[1] = 1;
+      }
     else if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
       {
         vtk_cell_id[0] = vtk_cell_type[dim];
-- 
2.39.5