Towards making DataOutFaces work with simplices.

source/numerics/data_out.cc

@@ -153,7 +153,7 @@ DataOut<dim, spacedim>::build_one_patch(
       // explicitly given
       patch.points_are_available = true;
 
-      // then resize the patch.data member in order to have enough memory for
+      // then size the patch.data member in order to have enough memory for
       // the quadrature points as well, and copy the quadrature points there
       const std::vector<Point<spacedim>> &q_points =
         fe_patch_values.get_quadrature_points();

source/numerics/data_out_faces.cc

@@ -99,31 +99,41 @@ DataOutFaces<dim, spacedim>::build_one_patch(
   internal::DataOutFacesImplementation::ParallelData<dim, spacedim> &data,
   DataOutBase::Patch<patch_dim, patch_spacedim> &                    patch)
 {
-  Assert(cell_and_face->first->is_locally_owned(), ExcNotImplemented());
+  const cell_iterator cell        = cell_and_face->first;
+  const unsigned int  face_number = cell_and_face->second;
 
-  // we use the mapping to transform the vertices. However, the mapping works
+  Assert(cell->is_locally_owned(), ExcNotImplemented());
+
+  // First set the kind of object we are dealing with here in the 'patch'
+  // object.
+  patch.reference_cell = cell->face(face_number)->reference_cell();
+
+  // We use the mapping to transform the vertices. However, the mapping works
   // on cells, not faces, so transform the face vertex to a cell vertex, that
   // to a unit cell vertex and then, finally, that to the mapped vertex. In
   // most cases this complicated procedure will be the identity.
-  for (unsigned int vertex = 0;
-       vertex < GeometryInfo<dim - 1>::vertices_per_cell;
-       ++vertex)
-    patch.vertices[vertex] =
-      data.mapping_collection[0].transform_unit_to_real_cell(
-        cell_and_face->first,
-        GeometryInfo<dim>::unit_cell_vertex(
-          GeometryInfo<dim>::face_to_cell_vertices(
-            cell_and_face->second,
+  for (const unsigned int vertex : cell->face(face_number)->vertex_indices())
+    {
+      const Point<dim> vertex_reference_coordinates =
+        cell->reference_cell().template vertex<dim>(
+          cell->reference_cell().face_to_cell_vertices(
+            face_number,
             vertex,
-            cell_and_face->first->face_orientation(cell_and_face->second),
-            cell_and_face->first->face_flip(cell_and_face->second),
-            cell_and_face->first->face_rotation(cell_and_face->second))));
+            cell->face_orientation(face_number) +
+              4 * cell->face_flip(face_number) +
+              2 * cell->face_rotation(face_number)));
+
+      const Point<dim> vertex_real_coordinates =
+        data.mapping_collection[0].transform_unit_to_real_cell(
+          cell, vertex_reference_coordinates);
+
+      patch.vertices[vertex] = vertex_real_coordinates;
+    }
+
 
   if (data.n_datasets > 0)
     {
-      data.reinit_all_fe_values(this->dof_data,
-                                cell_and_face->first,
-                                cell_and_face->second);
+      data.reinit_all_fe_values(this->dof_data, cell, face_number);
       const FEValuesBase<dim> &fe_patch_values = data.get_present_fe_values(0);
 
       const unsigned int n_q_points = fe_patch_values.n_quadrature_points;
@@ -132,9 +142,9 @@ DataOutFaces<dim, spacedim>::build_one_patch(
       Assert(patch.space_dim == dim, ExcInternalError());
       const std::vector<Point<dim>> &q_points =
         fe_patch_values.get_quadrature_points();
-      // resize the patch.data member in order to have enough memory for the
+      // size the patch.data member in order to have enough memory for the
       // quadrature points as well
-      patch.data.reinit(data.n_datasets + dim, patch.data.size(1));
+      patch.data.reinit(data.n_datasets + dim, q_points.size());
       // set the flag indicating that for this cell the points are explicitly
       // given
       patch.points_are_available = true;
@@ -194,9 +204,9 @@ DataOutFaces<dim, spacedim>::build_one_patch(
                       this_fe_patch_values.get_normal_vectors();
 
                   const typename DoFHandler<dim, spacedim>::active_cell_iterator
-                    dh_cell(&cell_and_face->first->get_triangulation(),
-                            cell_and_face->first->level(),
-                            cell_and_face->first->index(),
+                    dh_cell(&cell->get_triangulation(),
+                            cell->level(),
+                            cell->index(),
                             this->dof_data[dataset]->dof_handler);
                   data.patch_values_scalar.template set_cell<dim>(dh_cell);
 
@@ -237,9 +247,9 @@ DataOutFaces<dim, spacedim>::build_one_patch(
                       this_fe_patch_values.get_normal_vectors();
 
                   const typename DoFHandler<dim, spacedim>::active_cell_iterator
-                    dh_cell(&cell_and_face->first->get_triangulation(),
-                            cell_and_face->first->level(),
-                            cell_and_face->first->index(),
+                    dh_cell(&cell->get_triangulation(),
+                            cell->level(),
+                            cell->index(),
                             this->dof_data[dataset]->dof_handler);
                   data.patch_values_system.template set_cell<dim>(dh_cell);
 
@@ -294,15 +304,14 @@ DataOutFaces<dim, spacedim>::build_one_patch(
           // belongs in order to access the cell data. this is not readily
           // available, so choose the following rather inefficient way:
           Assert(
-            cell_and_face->first->is_active(),
+            cell->is_active(),
             ExcMessage(
               "The current function is trying to generate cell-data output "
               "for a face that does not belong to an active cell. This is "
               "not supported."));
           const unsigned int cell_number = std::distance(
             this->triangulation->begin_active(),
-            typename Triangulation<dim, spacedim>::active_cell_iterator(
-              cell_and_face->first));
+            typename Triangulation<dim, spacedim>::active_cell_iterator(cell));
 
           const double value = this->cell_data[dataset]->get_cell_data_value(
             cell_number,
@@ -395,9 +404,6 @@ DataOutFaces<dim, spacedim>::build_patches(
     update_flags);
   DataOutBase::Patch<patch_dim, patch_spacedim> sample_patch;
   sample_patch.n_subdivisions = n_subdivisions;
-  sample_patch.data.reinit(n_datasets,
-                           Utilities::fixed_power<patch_dim>(n_subdivisions +
-                                                             1));
 
   // now build the patches in parallel
   WorkStream::run(