Do not preset Patch::reference_cell.

include/deal.II/base/bounding_box_data_out.h

@@ -141,9 +141,10 @@ BoundingBoxDataOut<dim>::build_patches(
       boost::geometry::convert(getter(*value), box);
       for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; ++v)
         {
-          patches[i].vertices[v]          = box.vertex(v);
-          patches[i].patch_index          = i;
-          patches[i].n_subdivisions       = 1;
+          patches[i].vertices[v]    = box.vertex(v);
+          patches[i].patch_index    = i;
+          patches[i].n_subdivisions = 1;
+          patches[i].reference_cell = ReferenceCells::get_hypercube<dim>();
           patches[i].points_are_available = false;
         }
       ++i;

include/deal.II/lac/matrix_out.h

@@ -349,12 +349,15 @@ MatrixOut::build_patches(const Matrix &     matrix,
   for (size_type i = 0; i < gridpoints_y; ++i)
     for (size_type j = 0; j < gridpoints_x; ++j, ++index)
       {
+        patches[index].n_subdivisions = 1;
+        patches[index].reference_cell = ReferenceCells::Quadrilateral;
+
         // within each patch, order the points in such a way that if some
         // graphical output program (such as gnuplot) plots the quadrilaterals
         // as two triangles, then the diagonal of the quadrilateral which cuts
         // it into the two printed triangles is parallel to the diagonal of the
         // matrix, rather than perpendicular to it. this has the advantage that,
-        // for example, the unit matrix is plotted as a straight rim, rather
+        // for example, the unit matrix is plotted as a straight ridge, rather
         // than as a series of bumps and valleys along the diagonal
         patches[index].vertices[0](0) = j;
         patches[index].vertices[0](1) = -static_cast<signed int>(i);

source/base/data_out_base.cc

@@ -959,6 +959,12 @@ namespace
     n_cells = 0;
     for (const auto &patch : patches)
       {
+        Assert(patch.reference_cell != ReferenceCells::Invalid,
+               ExcMessage(
+                 "The reference cell for this patch is set to 'Invalid', "
+                 "but that is clearly not a valid choice. Did you forget "
+                 "to set the reference cell for the patch?"));
+
         // The following formula doesn't hold for non-tensor products.
         if (patch.reference_cell == ReferenceCells::get_hypercube<dim>())
           {
@@ -1964,7 +1970,7 @@ namespace DataOutBase
     : patch_index(no_neighbor)
     , n_subdivisions(1)
     , points_are_available(false)
-    , reference_cell(ReferenceCells::get_hypercube<dim>())
+    , reference_cell(ReferenceCells::Invalid)
   // all the other data has a constructor of its own, except for the "neighbors"
   // field, which we set to invalid values.
   {

source/numerics/data_out_rotation.cc

@@ -524,6 +524,8 @@ DataOutRotation<dim, spacedim>::build_patches(
   for (unsigned int i = 0; i < new_patches.size(); ++i)
     {
       new_patches[i].n_subdivisions = n_subdivisions;
+      new_patches[i].reference_cell = ReferenceCells::get_hypercube<dim + 1>();
+
       new_patches[i].data.reinit(
         n_datasets, Utilities::fixed_power<patch_dim>(n_subdivisions + 1));
     }

source/numerics/data_out_stack.cc

@@ -310,6 +310,7 @@ DataOutStack<dim, spacedim, void>::build_patches(
   // the time direction) points
   dealii::DataOutBase::Patch<patch_dim, patch_spacedim> default_patch;
   default_patch.n_subdivisions = n_subdivisions;
+  default_patch.reference_cell = ReferenceCells::get_hypercube<dim + 1>();
   default_patch.data.reinit(n_datasets, n_q_points * (n_subdivisions + 1));
   patches.insert(patches.end(), n_patches, default_patch);
 

tests/base/functions_04.cc

@@ -61,6 +61,7 @@ check_function(const Functions::FlowFunction<dim> &f,
     patches[0].neighbors[i] = numbers::invalid_unsigned_int;
   patches[0].patch_index          = 0;
   patches[0].n_subdivisions       = sub;
+  patches[0].reference_cell       = ReferenceCells::get_hypercube<dim>();
   patches[0].points_are_available = false;
 
   vertex_number = 1;

tests/data_out/data_out_base.cc

@@ -161,6 +161,7 @@ create_patches(std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
       DataOutBase::Patch<dim, spacedim> &p = patches[c];
       p.patch_index                        = c;
       p.n_subdivisions                     = nsub;
+      p.reference_cell = ReferenceCells::get_hypercube<dim>();
 
       for (unsigned int i = 0; i < ncells; ++i)
         for (unsigned int j = 0; j < spacedim; ++j)

tests/data_out/data_out_hdf5_02.cc

@@ -41,6 +41,7 @@ create_patches(std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
       const unsigned int nsubp = nsub + 1;
 
       patch.n_subdivisions = nsub;
+      patch.reference_cell = ReferenceCells::get_hypercube<dim>();
       for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
         for (unsigned int d = 0; d < spacedim; ++d)
           patch.vertices[v](d) =

tests/data_out/patches.h

@@ -35,6 +35,8 @@ create_patches(std::vector<DataOutBase::Patch<dim, spacedim>> &patches)
       const unsigned int nsubp = nsub + 1;
 
       patch.n_subdivisions = nsub;
+      patch.reference_cell = ReferenceCells::get_hypercube<dim>();
+
       for (const unsigned int v : GeometryInfo<dim>::vertex_indices())
         for (unsigned int d = 0; d < spacedim; ++d)
           patch.vertices[v](d) =
@@ -98,6 +100,7 @@ create_continuous_patches(std::vector<DataOutBase::Patch<dim, dim>> &patches,
         {
           DataOutBase::Patch<dim, dim> patch;
           patch.n_subdivisions = n_sub;
+          patch.reference_cell = ReferenceCells::get_hypercube<dim>();
           for (unsigned int k = 0; k < trapez.size(); ++k)
             {
               Point<dim> p = trapez.point(k);