Split MappingCartesian::compute_fill into several functions.

Presently MappingCartesian::compute_fill is called from the various
fill_fe*_values functions to compute quadrature points and normals.
compute_fill is called with dummy-normals in fill_fe_values and with
invalid face indices whenever face_index and subface_index has no
meaning. This makes the code hard to understand and the asserts that
check the arguments complicated. Split the compute_fill function into
several smaller functions to avoid passing around dummy or invalid
data.

diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h
index 6a51c77cc9cd369c3b84d9c7bd3d3703dfe674b7..7b18a9e518cf06ed91e6d2193dcf76c7afb9c730 100644 (file)
--- a/include/deal.II/fe/mapping_cartesian.h
+++ b/include/deal.II/fe/mapping_cartesian.h
@@ -19,6 +19,8 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/qprojector.h>
+
 #include <deal.II/fe/mapping.h>
 
 #include <cmath>
@@ -255,24 +257,77 @@ private:
    * @}
    */
 
+  /**
+   * Update the cell_extents field of the incoming InternalData object with the
+   * size of the incoming cell.
+   */
+  void
+  update_cell_extents(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const CellSimilarity::Similarity                            cell_similarity,
+    const InternalData &                                        data) const;
+
+  /**
+   * Compute the quadrature points if the UpdateFlags of the incoming
+   * InternalData object say that they should be updated.
+   *
+   * Called from fill_fe_values.
+   */
+  void
+  maybe_update_cell_quadrature_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const InternalData &                                        data,
+    std::vector<Point<dim>> &quadrature_points) const;
+
+  /**
+   * Compute the quadrature points if the UpdateFlags of the incoming
+   * InternalData object say that they should be updated.
+   *
+   * Called from fill_fe_face_values.
+   */
+  void
+  maybe_update_face_quadrature_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_no,
+    const InternalData &                                        data,
+    std::vector<Point<dim>> &quadrature_points) const;
 
+  /**
+   * Compute the quadrature points if the UpdateFlags of the incoming
+   * InternalData object say that they should be updated.
+   *
+   * Called from fill_fe_subface_values.
+   */
+  void
+  maybe_update_subface_quadrature_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int                                          face_no,
+    const unsigned int                                          sub_no,
+    const InternalData &                                        data,
+    std::vector<Point<dim>> &quadrature_points) const;
 
   /**
-   * Do the computation for the <tt>fill_*</tt> functions.
+   * Transform quadrature points in InternalData to real space by scaling unit
+   * coordinates with cell_extends in each direction.
+   *
+   * Called from the various maybe_update_*_quadrature_points functions.
    */
   void
-  compute_fill(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-               const unsigned int               face_no,
-               const unsigned int               sub_no,
-               const CellSimilarity::Similarity cell_similarity,
-               const InternalData &             data,
-               std::vector<Point<dim>> &        quadrature_points,
-               std::vector<Tensor<1, dim>> &    normal_vectors) const;
+  transform_quadrature_points(
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const InternalData &                                        data,
+    const typename QProjector<dim>::DataSetDescriptor &         offset,
+    std::vector<Point<dim>> &quadrature_points) const;
 
   /**
-   * Value to indicate that a given face or subface number is invalid.
+   * Compute the normal vectors if the UpdateFlags of the incoming InternalData
+   * object say that they should be updated.
    */
-  static const unsigned int invalid_face_number = numbers::invalid_unsigned_int;
+  void
+  maybe_update_normal_vectors(
+    const unsigned int           face_no,
+    const InternalData &         data,
+    std::vector<Tensor<1, dim>> &normal_vectors) const;
 };
 
 /*@}*/

diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc
index 2b084f528d9b281eb7d1284cc13d5f489a5f0854..96aa3d35741239674b481a9c3fd544c33b89a6e4 100644 (file)
--- a/source/fe/mapping_cartesian.cc
+++ b/source/fe/mapping_cartesian.cc
@@ -38,10 +38,6 @@
 DEAL_II_NAMESPACE_OPEN
 
 
-template <int dim, int spacedim>
-const unsigned int MappingCartesian<dim, spacedim>::invalid_face_number;
-
-
 
 template <int dim, int spacedim>
 MappingCartesian<dim, spacedim>::InternalData::InternalData(
@@ -164,64 +160,18 @@ MappingCartesian<dim, spacedim>::get_subface_data(
 
 template <int dim, int spacedim>
 void
-MappingCartesian<dim, spacedim>::compute_fill(
+MappingCartesian<dim, spacedim>::update_cell_extents(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const unsigned int                                          face_no,
-  const unsigned int                                          sub_no,
   const CellSimilarity::Similarity                            cell_similarity,
-  const InternalData &                                        data,
-  std::vector<Point<dim>> &                                   quadrature_points,
-  std::vector<Tensor<1, dim>> &normal_vectors) const
+  const InternalData &                                        data) const
 {
-  const UpdateFlags update_flags = data.update_each;
-
-  // some more sanity checks
-  if (face_no != invalid_face_number)
-    {
-      // Add 1 on both sides of
-      // assertion to avoid compiler
-      // warning about testing
-      // unsigned int < 0 in 1d.
-      Assert(face_no + 1 < GeometryInfo<dim>::faces_per_cell + 1,
-             ExcIndexRange(face_no, 0, GeometryInfo<dim>::faces_per_cell));
-
-      // We would like to check for
-      // sub_no < cell->face(face_no)->n_children(),
-      // but unfortunately the current
-      // function is also called for
-      // faces without children (see
-      // tests/fe/mapping.cc). Therefore,
-      // we must use following workaround
-      // of two separate assertions
-      Assert((sub_no == invalid_face_number) ||
-               cell->face(face_no)->has_children() ||
-               (sub_no + 1 < GeometryInfo<dim>::max_children_per_face + 1),
-             ExcIndexRange(sub_no,
-                           0,
-                           GeometryInfo<dim>::max_children_per_face));
-      Assert((sub_no == invalid_face_number) ||
-               !cell->face(face_no)->has_children() ||
-               (sub_no < cell->face(face_no)->n_children()),
-             ExcIndexRange(sub_no, 0, cell->face(face_no)->n_children()));
-    }
-  else
-    // invalid face number, so
-    // subface should be invalid as
-    // well
-    Assert(sub_no == invalid_face_number, ExcInternalError());
-
-  // let @p{start} be the origin of a
-  // local coordinate system. it is
-  // chosen as the (lower) left
-  // vertex
-  const Point<dim> start = cell->vertex(0);
-
   // Compute start point and sizes
   // along axes.  Strange vertex
   // numbering makes this complicated
   // again.
   if (cell_similarity != CellSimilarity::translation)
     {
+      const Point<dim> &start = cell->vertex(0);
       switch (dim)
         {
           case 1:
@@ -240,47 +190,121 @@ MappingCartesian<dim, spacedim>::compute_fill(
             Assert(false, ExcNotImplemented());
         }
     }
+}
+
 
 
-  // transform quadrature point. this
-  // is obtained simply by scaling
-  // unit coordinates with lengths in
-  // each direction
-  if (update_flags & update_quadrature_points)
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::maybe_update_cell_quadrature_points(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const InternalData &                                        data,
+  std::vector<Point<dim>> &quadrature_points) const
+{
+  if (data.update_each & update_quadrature_points)
     {
       const typename QProjector<dim>::DataSetDescriptor offset =
-        (face_no == invalid_face_number ?
-           QProjector<dim>::DataSetDescriptor::cell() :
-           (sub_no == invalid_face_number ?
-              // called from FEFaceValues
-              QProjector<dim>::DataSetDescriptor::face(
-                face_no,
-                cell->face_orientation(face_no),
-                cell->face_flip(face_no),
-                cell->face_rotation(face_no),
-                quadrature_points.size()) :
-              // called from FESubfaceValues
-              QProjector<dim>::DataSetDescriptor::subface(
-                face_no,
-                sub_no,
-                cell->face_orientation(face_no),
-                cell->face_flip(face_no),
-                cell->face_rotation(face_no),
-                quadrature_points.size(),
-                cell->subface_case(face_no))));
-
-      for (unsigned int i = 0; i < quadrature_points.size(); ++i)
-        {
-          quadrature_points[i] = start;
-          for (unsigned int d = 0; d < dim; ++d)
-            quadrature_points[i](d) +=
-              data.cell_extents[d] * data.quadrature_points[i + offset](d);
-        }
+        QProjector<dim>::DataSetDescriptor::cell();
+
+      transform_quadrature_points(cell, data, offset, quadrature_points);
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::maybe_update_face_quadrature_points(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const InternalData &                                        data,
+  std::vector<Point<dim>> &quadrature_points) const
+{
+  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
+
+  if (data.update_each & update_quadrature_points)
+    {
+      const typename QProjector<dim>::DataSetDescriptor offset =
+        QProjector<dim>::DataSetDescriptor::face(face_no,
+                                                 cell->face_orientation(
+                                                   face_no),
+                                                 cell->face_flip(face_no),
+                                                 cell->face_rotation(face_no),
+                                                 quadrature_points.size());
+
+
+      transform_quadrature_points(cell, data, offset, quadrature_points);
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::maybe_update_subface_quadrature_points(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const unsigned int                                          face_no,
+  const unsigned int                                          sub_no,
+  const InternalData &                                        data,
+  std::vector<Point<dim>> &quadrature_points) const
+{
+  AssertIndexRange(face_no, GeometryInfo<dim>::faces_per_cell);
+  AssertIndexRange(sub_no, GeometryInfo<dim>::max_children_per_face);
+  if (cell->face(face_no)->has_children())
+    {
+      AssertIndexRange(sub_no, cell->face(face_no)->n_children());
+    }
+
+  if (data.update_each & update_quadrature_points)
+    {
+      const typename QProjector<dim>::DataSetDescriptor offset =
+        QProjector<dim>::DataSetDescriptor::subface(
+          face_no,
+          sub_no,
+          cell->face_orientation(face_no),
+          cell->face_flip(face_no),
+          cell->face_rotation(face_no),
+          quadrature_points.size(),
+          cell->subface_case(face_no));
+
+      transform_quadrature_points(cell, data, offset, quadrature_points);
+    }
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::transform_quadrature_points(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const InternalData &                                        data,
+  const typename QProjector<dim>::DataSetDescriptor &         offset,
+  std::vector<Point<dim>> &quadrature_points) const
+{
+  // let @p{start} be the origin of a local coordinate system. it is chosen
+  // as the (lower) left vertex
+  const Point<dim> &start = cell->vertex(0);
+
+  for (unsigned int i = 0; i < quadrature_points.size(); ++i)
+    {
+      quadrature_points[i] = start;
+      for (unsigned int d = 0; d < dim; ++d)
+        quadrature_points[i](d) +=
+          data.cell_extents[d] * data.quadrature_points[i + offset](d);
     }
+}
 
 
+
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::maybe_update_normal_vectors(
+  const unsigned int           face_no,
+  const InternalData &         data,
+  std::vector<Tensor<1, dim>> &normal_vectors) const
+{
   // compute normal vectors. All normals on a face have the same value.
-  if (update_flags & update_normal_vectors)
+  if (data.update_each & update_normal_vectors)
     {
       Assert(face_no < GeometryInfo<dim>::faces_per_cell, ExcInternalError());
       std::fill(normal_vectors.begin(),
@@ -307,15 +331,12 @@ MappingCartesian<dim, spacedim>::fill_fe_values(
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
 
-  std::vector<Tensor<1, dim>> dummy;
 
-  compute_fill(cell,
-               invalid_face_number,
-               invalid_face_number,
-               cell_similarity,
-               data,
-               output_data.quadrature_points,
-               dummy);
+  update_cell_extents(cell, cell_similarity, data);
+
+  maybe_update_cell_quadrature_points(cell,
+                                      data,
+                                      output_data.quadrature_points);
 
   // compute Jacobian determinant. all values are equal and are the
   // product of the local lengths in each coordinate direction
@@ -420,13 +441,14 @@ MappingCartesian<dim, spacedim>::fill_fe_face_values(
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
 
-  compute_fill(cell,
-               face_no,
-               invalid_face_number,
-               CellSimilarity::none,
-               data,
-               output_data.quadrature_points,
-               output_data.normal_vectors);
+  update_cell_extents(cell, CellSimilarity::none, data);
+
+  maybe_update_face_quadrature_points(cell,
+                                      face_no,
+                                      data,
+                                      output_data.quadrature_points);
+
+  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
 
   // first compute Jacobian determinant, which is simply the product
   // of the local lengths since the jacobian is diagonal
@@ -523,13 +545,12 @@ MappingCartesian<dim, spacedim>::fill_fe_subface_values(
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
 
-  compute_fill(cell,
-               face_no,
-               subface_no,
-               CellSimilarity::none,
-               data,
-               output_data.quadrature_points,
-               output_data.normal_vectors);
+  update_cell_extents(cell, CellSimilarity::none, data);
+
+  maybe_update_subface_quadrature_points(
+    cell, face_no, subface_no, data, output_data.quadrature_points);
+
+  maybe_update_normal_vectors(face_no, data, output_data.normal_vectors);
 
   // first compute Jacobian determinant, which is simply the product
   // of the local lengths since the jacobian is diagonal
@@ -898,6 +919,7 @@ MappingCartesian<dim, spacedim>::transform(
 }
 
 
+
 template <int dim, int spacedim>
 void
 MappingCartesian<dim, spacedim>::transform(
@@ -935,6 +957,8 @@ MappingCartesian<dim, spacedim>::transform(
     }
 }
 
+
+
 template <int dim, int spacedim>
 void
 MappingCartesian<dim, spacedim>::transform(
@@ -1022,6 +1046,7 @@ MappingCartesian<dim, spacedim>::transform(
 }
 
 
+
 template <int dim, int spacedim>
 Point<spacedim>
 MappingCartesian<dim, spacedim>::transform_unit_to_real_cell(
@@ -1090,6 +1115,7 @@ MappingCartesian<dim, spacedim>::transform_real_to_unit_cell(
 }
 
 
+
 template <int dim, int spacedim>
 std::unique_ptr<Mapping<dim, spacedim>>
 MappingCartesian<dim, spacedim>::clone() const