Make DataOutFaces work with simplices.

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 b48c7f7019834340c75035a152a03c6729c5933d..ea788c8023755344bc6e9e0ba38cba566af5d79a 100644 (file)
--- a/include/deal.II/numerics/data_out_dof_data.templates.h
+++ b/include/deal.II/numerics/data_out_dof_data.templates.h
@@ -282,6 +282,7 @@ namespace internal
                                                ReferenceCells::Pyramid);
            }));
 
+      // Decide whether we want to work on cell or face FE(Face)Values objects:
       if (use_face_values == false)
         {
           std::unique_ptr<dealii::Quadrature<dim>> quadrature_simplex;
@@ -325,21 +326,23 @@ namespace internal
           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 one of the previous finite elements is equal to the
+              // present one. If so, re-use the FEValues object.
               for (unsigned int j = 0; j < i; ++j)
                 if (finite_elements[i].get() == finite_elements[j].get())
                   {
                     x_fe_values[i] = x_fe_values[j];
                     break;
                   }
+
+              // If none was found, create an FEValues object:
               if (x_fe_values[i].get() == nullptr)
                 {
                   dealii::hp::QCollection<dim> quadrature;
 
                   for (unsigned int j = 0; j < finite_elements[i]->size(); ++j)
                     {
-                      const auto reference_cell =
+                      const dealii::ReferenceCell reference_cell =
                         (*finite_elements[i])[j].reference_cell();
 
                       if (reference_cell.is_hyper_cube())
@@ -379,31 +382,84 @@ namespace internal
         }
       else // build FEFaceValues objects instead
         {
-          dealii::hp::QCollection<dim - 1> quadrature(
-            QIterated<dim - 1>(QTrapezoid<1>(), n_subdivisions));
+          std::unique_ptr<dealii::Quadrature<dim - 1>> quadrature_simplex;
+          std::unique_ptr<dealii::Quadrature<dim - 1>> quadrature_hypercube;
+
+          // See whether we need simplex or hypercube quadrature formulas.
+          // This is only an issue in 3d (in 2d every face integral is just
+          // a line integral, so we can deal with that via the usual hypercube
+          // quadrature rule).
+          if ((dim == 3) &&
+              (needs_simplex_setup || needs_pyramid_setup || needs_wedge_setup))
+            {
+              quadrature_simplex = std::make_unique<Quadrature<dim - 1>>(
+                generate_simplex_evaluation_points<dim - 1>(n_subdivisions));
+            }
+
+          if ((dim < 3) || (needs_hypercube_setup || needs_pyramid_setup ||
+                            needs_wedge_setup))
+            {
+              quadrature_hypercube =
+                std::make_unique<QIterated<dim - 1>>(QTrapezoid<1>(),
+                                                     n_subdivisions);
+            }
 
           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 one of the previous finite elements is equal to the
+              // present one. If so, re-use the FEValues object.
               for (unsigned int j = 0; j < i; ++j)
                 if (finite_elements[i].get() == finite_elements[j].get())
                   {
                     x_fe_face_values[i] = x_fe_face_values[j];
                     break;
                   }
+
+              // If none was found, create an FEFaceValues object:
               if (x_fe_face_values[i].get() == nullptr)
-                x_fe_face_values[i] =
-                  std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
-                    mapping_collection,
-                    *finite_elements[i],
-                    quadrature,
-                    update_flags);
+                {
+                  dealii::hp::QCollection<dim - 1> quadrature;
+
+                  for (unsigned int j = 0; j < finite_elements[i]->size(); ++j)
+                    {
+                      const dealii::ReferenceCell reference_cell =
+                        (*finite_elements[i])[j].reference_cell();
+
+                      // In 1d/2d and for hypercube/wedge/pyramid elements, we
+                      // need hypercube quadratures.
+                      if ((dim < 3) ||
+                          (reference_cell.is_hyper_cube() ||
+                           (reference_cell == dealii::ReferenceCells::Wedge) ||
+                           (reference_cell == dealii::ReferenceCells::Pyramid)))
+                        quadrature.push_back(*quadrature_hypercube);
+
+                      // In 3d, if the element is for simplex/wedge/pyramid
+                      // cells, then we also need simplex quadratures
+                      if ((dim == 3) &&
+                          (reference_cell.is_simplex() ||
+                           (reference_cell == dealii::ReferenceCells::Wedge) ||
+                           (reference_cell == dealii::ReferenceCells::Pyramid)))
+                        quadrature.push_back(*quadrature_simplex);
+                    }
+
+                  x_fe_face_values[i] =
+                    std::make_shared<dealii::hp::FEFaceValues<dim, spacedim>>(
+                      mapping_collection,
+                      *finite_elements[i],
+                      quadrature,
+                      update_flags);
+                }
             }
 
           // Return maximal number of evaluation points:
-          return quadrature[0].size();
+          return std::max(
+            {(dim == 3) && (needs_simplex_setup || needs_pyramid_setup ||
+                            needs_wedge_setup) ?
+               quadrature_simplex->size() :
+               0,
+             (dim < 3) || needs_hypercube_setup ? quadrature_hypercube->size() :
+                                                  0});
         }
     }
 
@@ -496,11 +552,14 @@ namespace internal
             {
               if (cell->is_active())
                 {
-                  typename DoFHandler<dim, spacedim>::active_cell_iterator
+                  const typename DoFHandler<dim, spacedim>::active_cell_iterator
                     dh_cell(&cell->get_triangulation(),
                             cell->level(),
                             cell->index(),
                             dof_data[dataset]->dof_handler);
+
+                  // Check whether we need cell or face FEValues objects by
+                  // testing which of the two arrays actually has any content.
                   if (x_fe_values.empty())
                     {
                       AssertIndexRange(face, GeometryInfo<dim>::faces_per_cell);
@@ -513,6 +572,11 @@ namespace internal
                 x_fe_values[dataset]->reinit(cell);
             }
         }
+
+      // If there is are no DoF-associated data (just cell-associated ones),
+      // then the loop above will not execute any iterations. In that case,
+      // do the initialization for the first FE(Face)Values object by
+      // hand, using only the (triangulation) cell without a DoFHandler.
       if (dof_data.empty())
         {
           if (x_fe_values.empty())
@@ -533,6 +597,9 @@ namespace internal
       const unsigned int dataset) const
     {
       AssertIndexRange(dataset, finite_elements.size());
+
+      // Check whether we need cell or face FEValues objects by testing
+      // which of the two arrays actually has any content.
       if (x_fe_values.empty())
         return x_fe_face_values[dataset]->get_present_fe_values();
       else

diff --git a/tests/fail/step-04-data_out_faces.output b/tests/fail/step-04-data_out_faces.output
index 484c320883fff536050ede820ba5f69be77b2add..a795ce86268bac26171a4626381ba0e2bc3072e5 100644 (file)
--- a/tests/fail/step-04-data_out_faces.output
+++ b/tests/fail/step-04-data_out_faces.output
@@ -1,4 +1,4 @@
-
+JobId strange.math.colostate.edu Tue Oct 26 17:19:44 2021
 DEAL::Solving problem in 2 space dimensions.
 DEAL::   Number of active cells: 8
 DEAL::   Total number of cells: 8
@@ -7,7 +7,7 @@ DEAL:cg::Starting value 4.18267
 DEAL:cg::Convergence step 1 value 0.00000
 DEAL::   1 CG iterations needed to obtain convergence.
 # vtk DataFile Version 3.0
-#This file was generated 
+#This file was generated by the deal.II library on 2021/10/26 at 17:19:44
 ASCII
 DATASET UNSTRUCTURED_GRID
 
@@ -53,137 +53,113 @@ DEAL:cg::Starting value 0.00000
 DEAL:cg::Convergence step 0 value 0.00000
 DEAL::   0 CG iterations needed to obtain convergence.
 # vtk DataFile Version 3.0
-#This file was generated 
+#This file was generated by the deal.II library on 2021/10/26 at 17:19:44
 ASCII
 DATASET UNSTRUCTURED_GRID
 
-POINTS 96 double
+POINTS 72 double
 -1.00000 -1.00000 -1.00000
 1.00000 -1.00000 -1.00000
 0.00000 0.00000 -1.00000
-2.00000 0.00000 -1.00000
 1.00000 -1.00000 -1.00000
 -1.00000 -1.00000 -1.00000
 0.00000 -1.00000 0.00000
--2.00000 -1.00000 0.00000
 -1.00000 1.00000 -1.00000
 1.00000 1.00000 -1.00000
 0.00000 1.00000 0.00000
-2.00000 1.00000 0.00000
 1.00000 1.00000 -1.00000
 -1.00000 1.00000 -1.00000
 0.00000 0.00000 -1.00000
--2.00000 0.00000 -1.00000
 1.00000 1.00000 1.00000
 -1.00000 1.00000 1.00000
 0.00000 1.00000 0.00000
--2.00000 1.00000 0.00000
 -1.00000 1.00000 1.00000
 1.00000 1.00000 1.00000
 0.00000 0.00000 1.00000
-2.00000 0.00000 1.00000
 1.00000 -1.00000 1.00000
 -1.00000 -1.00000 1.00000
 0.00000 0.00000 1.00000
--2.00000 0.00000 1.00000
 -1.00000 -1.00000 1.00000
 1.00000 -1.00000 1.00000
 0.00000 -1.00000 0.00000
-2.00000 -1.00000 0.00000
 -1.00000 -1.00000 -1.00000
 -1.00000 1.00000 -1.00000
 -1.00000 0.00000 0.00000
--1.00000 2.00000 0.00000
 -1.00000 1.00000 -1.00000
 -1.00000 -1.00000 -1.00000
 0.00000 0.00000 -1.00000
-0.00000 -2.00000 -1.00000
 -1.00000 -1.00000 1.00000
 -1.00000 1.00000 1.00000
 0.00000 0.00000 1.00000
-0.00000 2.00000 1.00000
 -1.00000 1.00000 1.00000
 -1.00000 -1.00000 1.00000
 -1.00000 0.00000 0.00000
--1.00000 -2.00000 0.00000
 1.00000 -1.00000 1.00000
 0.00000 0.00000 1.00000
 1.00000 1.00000 1.00000
-0.00000 2.00000 1.00000
 1.00000 -1.00000 1.00000
 1.00000 1.00000 1.00000
 1.00000 0.00000 0.00000
-1.00000 2.00000 0.00000
 1.00000 -1.00000 -1.00000
 1.00000 0.00000 0.00000
 1.00000 1.00000 -1.00000
-1.00000 2.00000 0.00000
 1.00000 -1.00000 -1.00000
 1.00000 1.00000 -1.00000
 0.00000 0.00000 -1.00000
-0.00000 2.00000 -1.00000
 -1.00000 -1.00000 -1.00000
 -1.00000 0.00000 0.00000
 -1.00000 -1.00000 1.00000
--1.00000 0.00000 2.00000
 -1.00000 -1.00000 -1.00000
 -1.00000 -1.00000 1.00000
 0.00000 -1.00000 0.00000
-0.00000 -1.00000 2.00000
 1.00000 -1.00000 -1.00000
 1.00000 -1.00000 1.00000
 1.00000 0.00000 0.00000
-1.00000 0.00000 2.00000
 1.00000 -1.00000 1.00000
 1.00000 -1.00000 -1.00000
 0.00000 -1.00000 0.00000
-0.00000 -1.00000 -2.00000
 1.00000 1.00000 -1.00000
 1.00000 1.00000 1.00000
 0.00000 1.00000 0.00000
-0.00000 1.00000 2.00000
 1.00000 1.00000 1.00000
 1.00000 1.00000 -1.00000
 1.00000 0.00000 0.00000
-1.00000 0.00000 -2.00000
 -1.00000 1.00000 -1.00000
 -1.00000 1.00000 1.00000
 -1.00000 0.00000 0.00000
--1.00000 0.00000 2.00000
 -1.00000 1.00000 1.00000
 -1.00000 1.00000 -1.00000
 0.00000 1.00000 0.00000
-0.00000 1.00000 -2.00000
 
-CELLS 24 120
-4      0       1       3       2
-4      4       5       7       6
-4      8       9       11      10
-4      12      13      15      14
-4      16      17      19      18
-4      20      21      23      22
-4      24      25      27      26
-4      28      29      31      30
-4      32      33      35      34
-4      36      37      39      38
-4      40      41      43      42
-4      44      45      47      46
-4      48      49      51      50
-4      52      53      55      54
-4      56      57      59      58
-4      60      61      63      62
-4      64      65      67      66
-4      68      69      71      70
-4      72      73      75      74
-4      76      77      79      78
-4      80      81      83      82
-4      84      85      87      86
-4      88      89      91      90
-4      92      93      95      94
+CELLS 24 96
+       3       0       1       2
+       3       3       4       5
+       3       6       7       8
+       3       9       10      11
+       3       12      13      14
+       3       15      16      17
+       3       18      19      20
+       3       21      22      23
+       3       24      25      26
+       3       27      28      29
+       3       30      31      32
+       3       33      34      35
+       3       36      37      38
+       3       39      40      41
+       3       42      43      44
+       3       45      46      47
+       3       48      49      50
+       3       51      52      53
+       3       54      55      56
+       3       57      58      59
+       3       60      61      62
+       3       63      64      65
+       3       66      67      68
+       3       69      70      71
 
 CELL_TYPES 24
- 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9
-POINT_DATA 96
+ 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
+POINT_DATA 72
 SCALARS solution double 1
 LOOKUP_TABLE default
-3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 3.00000 3.00000 1.00000 1.00000 
+3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 1.00000 3.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000 3.00000 3.00000 1.00000