Enable DataOut for Wedges

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 f7b9b6fbe2b9e2aaefbdac97ea9dc23a358b6f87..e9f7ee2b61082946e2a9624348e576e3e8558470 100644 (file)
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -82,78 +82,94 @@ namespace internal
     {}
 
 
+
+    /**
+     * Set up vectors of FEValues and FEFaceValues needed inside of
+     * ParallelDataBase and return the maximum number of quadrature points
+     * needed to allocate enough memory for the scratch data.
+     */
     template <int dim, int spacedim>
-    ParallelDataBase<dim, spacedim>::ParallelDataBase(
-      const unsigned int               n_datasets,
-      const unsigned int               n_subdivisions,
-      const std::vector<unsigned int> &n_postprocessor_outputs,
-      const dealii::hp::MappingCollection<dim, spacedim> &mapping,
+    unsigned int
+    setup_parallel_data_base_internal(
+      const unsigned int                                  n_subdivisions,
+      const dealii::hp::MappingCollection<dim, spacedim> &mapping_collection,
       const std::vector<
         std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
         &               finite_elements,
       const UpdateFlags update_flags,
-      const bool        use_face_values)
-      : n_datasets(n_datasets)
-      , n_subdivisions(n_subdivisions)
-      , postprocessed_values(n_postprocessor_outputs.size())
-      , mapping_collection(mapping)
-      , finite_elements(finite_elements)
-      , update_flags(update_flags)
+      const bool        use_face_values,
+      std::vector<std::shared_ptr<dealii::hp::FEValues<dim, spacedim>>>
+        &x_fe_values,
+      std::vector<std::shared_ptr<dealii::hp::FEFaceValues<dim, spacedim>>>
+        &x_fe_face_values)
     {
       unsigned int n_q_points = 0;
       if (use_face_values == false)
         {
-          const bool do_simplex =
-            std::any_of(finite_elements.begin(),
-                        finite_elements.end(),
-                        [](const auto &fe) {
-                          for (unsigned int i = 0; i < fe->size(); ++i)
-                            if ((*fe)[i].reference_cell_type() !=
-                                ReferenceCell::get_hypercube(dim))
-                              return true;
-                          return false;
-                        });
-
-          const bool do_hex =
-            std::any_of(finite_elements.begin(),
-                        finite_elements.end(),
-                        [](const auto &fe) {
-                          for (unsigned int i = 0; i < fe->size(); ++i)
-                            if ((*fe)[i].reference_cell_type() ==
-                                ReferenceCell::get_hypercube(dim))
-                              return true;
-                          return false;
-                        });
-
-          std::unique_ptr<dealii::Quadrature<dim>> quadrature1;
-          std::unique_ptr<dealii::Quadrature<dim>> quadrature2;
-
-          if (do_simplex)
+          // determine if specific quadrature rules need to set up
+          bool needs_hypercube_setup = false;
+          bool needs_simplex_setup   = false;
+          bool needs_wedge_setup     = false;
+
+          for (const auto &fe : finite_elements)
+            for (unsigned int i = 0; i < fe->size(); ++i)
+              switch ((*fe)[i].reference_cell_type())
+                {
+                  case ReferenceCell::Type::Vertex:
+                  case ReferenceCell::Type::Line:
+                  case ReferenceCell::Type::Quad:
+                  case ReferenceCell::Type::Hex:
+                    needs_hypercube_setup |= true;
+                    break;
+                  case ReferenceCell::Type::Tri:
+                  case ReferenceCell::Type::Tet:
+                    needs_simplex_setup |= true;
+                    break;
+                  case ReferenceCell::Type::Wedge:
+                    needs_wedge_setup |= true;
+                    break;
+                  default:
+                    Assert(false, ExcNotImplemented());
+                }
+
+          std::unique_ptr<dealii::Quadrature<dim>> quadrature_simplex;
+          std::unique_ptr<dealii::Quadrature<dim>> quadrature_hypercube;
+          std::unique_ptr<dealii::Quadrature<dim>> quadrature_wedge;
+
+          if (needs_simplex_setup)
             {
-              Assert(1 <= n_subdivisions && n_subdivisions <= 2,
-                     ExcNotImplemented());
-              quadrature1.reset(new Quadrature<dim>(
-                Simplex::FE_P<dim, spacedim>(n_subdivisions == 1 ? 1 : 2)
-                  .get_unit_support_points()));
+              quadrature_simplex = std::make_unique<Quadrature<dim>>(
+                Simplex::FE_P<dim, spacedim>(n_subdivisions)
+                  .get_unit_support_points());
             }
 
-          if (do_hex)
+          if (needs_hypercube_setup)
             {
-              quadrature2.reset(new Quadrature<dim>(
-                QIterated<dim>(QTrapezoid<1>(), n_subdivisions)));
+              quadrature_hypercube =
+                std::make_unique<QIterated<dim>>(QTrapezoid<1>(),
+                                                 n_subdivisions);
             }
 
-          n_q_points = std::max(do_simplex ? quadrature1->size() : 0,
-                                do_hex ? quadrature2->size() : 1);
+          if (needs_wedge_setup)
+            {
+              quadrature_wedge = std::make_unique<Quadrature<dim>>(
+                Simplex::FE_WedgeP<dim, spacedim>(
+                  1 /*note: vtk only supports linear wedges*/)
+                  .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});
 
-          x_fe_values.resize(this->finite_elements.size());
-          for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
+          x_fe_values.resize(finite_elements.size());
+          for (unsigned int i = 0; i < finite_elements.size(); ++i)
             {
-              // check if there is a finite element that is equal to the present
-              // one, then we can re-use the FEValues object
+              // check if there is a finite element that is equal to the
+              // present one, then we can re-use the FEValues object
               for (unsigned int j = 0; j < i; ++j)
-                if (this->finite_elements[i].get() ==
-                    this->finite_elements[j].get())
+                if (finite_elements[i].get() == finite_elements[j].get())
                   {
                     x_fe_values[i] = x_fe_values[j];
                     break;
@@ -162,21 +178,37 @@ namespace internal
                 {
                   dealii::hp::QCollection<dim> quadrature;
 
-                  for (unsigned int j = 0; j < this->finite_elements[i]->size();
-                       ++j)
-                    if ((*this->finite_elements[i])[j].reference_cell_type() !=
-                        ReferenceCell::get_hypercube(dim))
-                      quadrature.push_back(*quadrature1);
-                    else
-                      quadrature.push_back(*quadrature2);
+                  for (unsigned int j = 0; j < finite_elements[i]->size(); ++j)
+                    switch ((*finite_elements[i])[j].reference_cell_type())
+                      {
+                        case ReferenceCell::Type::Vertex:
+                        case ReferenceCell::Type::Line:
+                        case ReferenceCell::Type::Quad:
+                        case ReferenceCell::Type::Hex:
+                          quadrature.push_back(*quadrature_hypercube);
+                          break;
+                        case ReferenceCell::Type::Tri:
+                        case ReferenceCell::Type::Tet:
+                          quadrature.push_back(*quadrature_simplex);
+                          break;
+                        case ReferenceCell::Type::Wedge:
+                          quadrature.push_back(*quadrature_wedge);
+                          break;
+                        default:
+                          Assert(false, ExcNotImplemented());
+                      }
 
                   x_fe_values[i] =
                     std::make_shared<dealii::hp::FEValues<dim, spacedim>>(
-                      this->mapping_collection,
-                      *this->finite_elements[i],
+                      mapping_collection,
+                      *finite_elements[i],
                       quadrature,
-                      do_simplex ? (update_flags | update_quadrature_points) :
-                                   update_flags);
+                      // 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) ?
+                        (update_flags | update_quadrature_points) :
+                        update_flags);
                 }
             }
         }
@@ -185,14 +217,13 @@ namespace internal
           dealii::hp::QCollection<dim - 1> quadrature(
             QIterated<dim - 1>(QTrapezoid<1>(), n_subdivisions));
           n_q_points = quadrature[0].size();
-          x_fe_face_values.resize(this->finite_elements.size());
-          for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
+          x_fe_face_values.resize(finite_elements.size());
+          for (unsigned int i = 0; i < finite_elements.size(); ++i)
             {
-              // check if there is a finite element that is equal to the present
-              // one, then we can re-use the FEValues object
+              // check if there is a finite element that is equal to the
+              // present one, then we can re-use the FEValues object
               for (unsigned int j = 0; j < i; ++j)
-                if (this->finite_elements[i].get() ==
-                    this->finite_elements[j].get())
+                if (finite_elements[i].get() == finite_elements[j].get())
                   {
                     x_fe_face_values[i] = x_fe_face_values[j];
                     break;
@@ -200,13 +231,45 @@ namespace internal
               if (x_fe_face_values[i].get() == nullptr)
                 x_fe_face_values[i] =
                   std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
-                    this->mapping_collection,
-                    *this->finite_elements[i],
+                    mapping_collection,
+                    *finite_elements[i],
                     quadrature,
-                    this->update_flags);
+                    update_flags);
             }
         }
 
+      return n_q_points;
+    }
+
+
+
+    template <int dim, int spacedim>
+    ParallelDataBase<dim, spacedim>::ParallelDataBase(
+      const unsigned int               n_datasets,
+      const unsigned int               n_subdivisions,
+      const std::vector<unsigned int> &n_postprocessor_outputs,
+      const dealii::hp::MappingCollection<dim, spacedim> &mapping,
+      const std::vector<
+        std::shared_ptr<dealii::hp::FECollection<dim, spacedim>>>
+        &               finite_elements,
+      const UpdateFlags update_flags,
+      const bool        use_face_values)
+      : n_datasets(n_datasets)
+      , n_subdivisions(n_subdivisions)
+      , postprocessed_values(n_postprocessor_outputs.size())
+      , mapping_collection(mapping)
+      , finite_elements(finite_elements)
+      , update_flags(update_flags)
+    {
+      const unsigned int n_q_points =
+        setup_parallel_data_base_internal(n_subdivisions,
+                                          mapping,
+                                          finite_elements,
+                                          update_flags,
+                                          use_face_values,
+                                          x_fe_values,
+                                          x_fe_face_values);
+
       patch_values_scalar.solution_values.resize(n_q_points);
       patch_values_scalar.solution_gradients.resize(n_q_points);
       patch_values_scalar.solution_hessians.resize(n_q_points);
@@ -238,110 +301,13 @@ namespace internal
       , finite_elements(data.finite_elements)
       , update_flags(data.update_flags)
     {
-      if (data.x_fe_values.empty() == false)
-        {
-          Assert(data.x_fe_face_values.empty() == true, ExcInternalError());
-
-          const bool do_simplex =
-            std::any_of(finite_elements.begin(),
-                        finite_elements.end(),
-                        [](const auto &fe) {
-                          for (unsigned int i = 0; i < fe->size(); ++i)
-                            if ((*fe)[i].reference_cell_type() !=
-                                ReferenceCell::get_hypercube(dim))
-                              return true;
-                          return false;
-                        });
-
-          const bool do_hex =
-            std::any_of(finite_elements.begin(),
-                        finite_elements.end(),
-                        [](const auto &fe) {
-                          for (unsigned int i = 0; i < fe->size(); ++i)
-                            if ((*fe)[i].reference_cell_type() ==
-                                ReferenceCell::get_hypercube(dim))
-                              return true;
-                          return false;
-                        });
-
-          std::unique_ptr<dealii::Quadrature<dim>> quadrature1;
-          std::unique_ptr<dealii::Quadrature<dim>> quadrature2;
-
-          if (do_simplex)
-            {
-              Assert(1 <= n_subdivisions && n_subdivisions <= 2,
-                     ExcNotImplemented());
-              quadrature1.reset(new Quadrature<dim>(
-                Simplex::FE_P<dim, spacedim>(n_subdivisions == 1 ? 1 : 2)
-                  .get_unit_support_points()));
-            }
-
-          if (do_hex)
-            {
-              quadrature2.reset(new Quadrature<dim>(
-                QIterated<dim>(QTrapezoid<1>(), n_subdivisions)));
-            }
-
-          x_fe_values.resize(this->finite_elements.size());
-
-          for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
-            {
-              // check if there is a finite element that is equal to the present
-              // one, then we can re-use the FEValues object
-              for (unsigned int j = 0; j < i; ++j)
-                if (this->finite_elements[i].get() ==
-                    this->finite_elements[j].get())
-                  {
-                    x_fe_values[i] = x_fe_values[j];
-                    break;
-                  }
-              if (x_fe_values[i].get() == nullptr)
-                {
-                  dealii::hp::QCollection<dim> quadrature;
-
-                  for (unsigned int j = 0; j < this->finite_elements[i]->size();
-                       ++j)
-                    if ((*this->finite_elements[i])[j].reference_cell_type() !=
-                        ReferenceCell::get_hypercube(dim))
-                      quadrature.push_back(*quadrature1);
-                    else
-                      quadrature.push_back(*quadrature2);
-
-                  x_fe_values[i] =
-                    std::make_shared<dealii::hp::FEValues<dim, spacedim>>(
-                      this->mapping_collection,
-                      *this->finite_elements[i],
-                      quadrature,
-                      do_simplex ? (update_flags | update_quadrature_points) :
-                                   update_flags);
-                }
-            }
-        }
-      else
-        {
-          dealii::hp::QCollection<dim - 1> quadrature(
-            QIterated<dim - 1>(QTrapezoid<1>(), n_subdivisions));
-          x_fe_face_values.resize(this->finite_elements.size());
-          for (unsigned int i = 0; i < this->finite_elements.size(); ++i)
-            {
-              // check if there is a finite element that is equal to the present
-              // one, then we can re-use the FEValues object
-              for (unsigned int j = 0; j < i; ++j)
-                if (this->finite_elements[i].get() ==
-                    this->finite_elements[j].get())
-                  {
-                    x_fe_face_values[i] = x_fe_face_values[j];
-                    break;
-                  }
-              if (x_fe_face_values[i].get() == nullptr)
-                x_fe_face_values[i] =
-                  std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
-                    this->mapping_collection,
-                    *this->finite_elements[i],
-                    quadrature,
-                    this->update_flags);
-            }
-        }
+      setup_parallel_data_base_internal(n_subdivisions,
+                                        mapping_collection,
+                                        finite_elements,
+                                        update_flags,
+                                        data.x_fe_values.empty(),
+                                        x_fe_values,
+                                        x_fe_face_values);
     }
 
 

diff --git a/source/base/data_out_base.cc b/source/base/data_out_base.cc
index 0d91d13b26fbde9b81635621d407cd3f082ca9f2..f3a51f184b978c254abad5c536dc3f869ded683a 100644 (file)
--- a/source/base/data_out_base.cc
+++ b/source/base/data_out_base.cc
@@ -628,6 +628,12 @@ namespace
         vtk_cell_id[0] = 24;
         vtk_cell_id[1] = 1;
       }
+    else if (patch.reference_cell_type == ReferenceCell::Type::Wedge &&
+             patch.data.n_cols() == 6)
+      {
+        vtk_cell_id[0] = 13;
+        vtk_cell_id[1] = 1;
+      }
     else if (patch.reference_cell_type == ReferenceCell::get_hypercube(dim))
       {
         vtk_cell_id[0] = vtk_cell_type[dim];