From: Martin Kronbichler <martin.kronbichler@rub.de>
Date: Sat, 6 Apr 2024 10:12:34 +0000 (+0200)
Subject: More efficient access via new Mapping::InternalDataBase::reinit
X-Git-Tag: v9.6.0-rc1~378^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9a3edfead47b355eee7f45d8e5aaca6965240dc9;p=dealii.git

More efficient access via new Mapping::InternalDataBase::reinit
---

diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h
index 74e0b11267..baba788e1c 100644
--- a/include/deal.II/fe/mapping.h
+++ b/include/deal.II/fe/mapping.h
@@ -59,6 +59,8 @@ namespace NonMatching
     template <int dim, int spacedim>
     class ComputeMappingDataHelper;
   }
+  template <int dim, int spacedim, typename Number>
+  class MappingInfo;
 } // namespace NonMatching
 
 
@@ -655,6 +657,18 @@ public:
      */
     InternalDataBase(const InternalDataBase &) = delete;
 
+    /**
+     * This function initializes the data fields related to evaluation of the
+     * mapping on cells, implemented by (derived) classes. This function is
+     * used both when setting up a field of this class for the first time or
+     * when a new Quadrature formula should be considered without creating an
+     * entirely new object. This is used when the number of evaluation points
+     * is different on each cell, e.g. when using FEPointEvaluation for
+     * handling particles or with certain non-matching problem settings.
+     */
+    virtual void
+    reinit(const UpdateFlags update_flags, const Quadrature<dim> &quadrature);
+
     /**
      * Virtual destructor for derived classes
      */
@@ -683,7 +697,6 @@ public:
     memory_consumption() const;
   };
 
-
 protected:
   /**
    * Given a set of update flags, compute which other quantities <i>also</i>
@@ -1323,6 +1336,8 @@ public:
   friend class FESubfaceValues<dim, spacedim>;
   friend class NonMatching::FEImmersedSurfaceValues<dim>;
   friend class NonMatching::internal::ComputeMappingDataHelper<dim, spacedim>;
+  template <int, int, typename>
+  friend class NonMatching::MappingInfo;
 };
 
 
diff --git a/include/deal.II/fe/mapping_cartesian.h b/include/deal.II/fe/mapping_cartesian.h
index 73f1f2006f..14dbaf6788 100644
--- a/include/deal.II/fe/mapping_cartesian.h
+++ b/include/deal.II/fe/mapping_cartesian.h
@@ -218,6 +218,11 @@ public:
      */
     InternalData(const Quadrature<dim> &quadrature);
 
+    // Documentation see Mapping::InternalDataBase.
+    virtual void
+    reinit(const UpdateFlags      update_flags,
+           const Quadrature<dim> &quadrature) override;
+
     /**
      * Return an estimate (in bytes) for the memory consumption of this object.
      */
@@ -243,7 +248,9 @@ public:
     mutable double volume_element;
 
     /**
-     * Vector of all quadrature points. Especially, all points on all faces.
+     * Location of quadrature points of faces or subfaces in 3d with all
+     * possible orientations. Can be accessed with the correct offset provided
+     * via QProjector::DataSetDescriptor. Not needed/used for cells.
      */
     std::vector<Point<dim>> quadrature_points;
   };
@@ -335,7 +342,8 @@ private:
   maybe_update_cell_quadrature_points(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const InternalData                                         &data,
-    std::vector<Point<dim>> &quadrature_points) const;
+    const ArrayView<const Point<dim>> &unit_quadrature_points,
+    std::vector<Point<dim>>           &quadrature_points) const;
 
   /**
    * Compute the quadrature points if the UpdateFlags of the incoming
@@ -364,19 +372,6 @@ private:
     const InternalData                                         &data,
     std::vector<Point<dim>> &quadrature_points) const;
 
-  /**
-   * Transform quadrature points in InternalData to real space by scaling unit
-   * coordinates with cell_extents in each direction.
-   *
-   * Called from the various maybe_update_*_quadrature_points functions.
-   */
-  void
-  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;
-
   /**
    * Compute the normal vectors if the UpdateFlags of the incoming InternalData
    * object say that they should be updated.
diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h
index ff265d1c55..48c9e18bd3 100644
--- a/include/deal.II/fe/mapping_fe.h
+++ b/include/deal.II/fe/mapping_fe.h
@@ -185,23 +185,15 @@ public:
      */
     InternalData(const FiniteElement<dim, spacedim> &fe);
 
-    /**
-     * Initialize the object's member variables related to cell data based on
-     * the given arguments.
-     *
-     * The function also calls compute_shape_function_values() to actually set
-     * the member variables related to the values and derivatives of the
-     * mapping shape functions.
-     */
-    void
-    initialize(const UpdateFlags      update_flags,
-               const Quadrature<dim> &quadrature,
-               const unsigned int     n_original_q_points);
+    // Documentation see Mapping::InternalDataBase.
+    virtual void
+    reinit(const UpdateFlags      update_flags,
+           const Quadrature<dim> &quadrature) override;
 
     /**
      * Initialize the object's member variables related to cell and face data
      * based on the given arguments. In order to initialize cell data, this
-     * function calls initialize().
+     * function calls reinit().
      */
     void
     initialize_face(const UpdateFlags      update_flags,
diff --git a/include/deal.II/fe/mapping_fe_field.h b/include/deal.II/fe/mapping_fe_field.h
index fd842faa55..1f90d82033 100644
--- a/include/deal.II/fe/mapping_fe_field.h
+++ b/include/deal.II/fe/mapping_fe_field.h
@@ -289,6 +289,11 @@ public:
     InternalData(const FiniteElement<dim, spacedim> &fe,
                  const ComponentMask                &mask);
 
+    // Documentation see Mapping::InternalDataBase.
+    virtual void
+    reinit(const UpdateFlags      update_flags,
+           const Quadrature<dim> &quadrature) override;
+
     /**
      * Shape function at quadrature point. Shape functions are in tensor
      * product order, so vertices must be reordered to obtain transformation.
@@ -359,6 +364,11 @@ public:
     virtual std::size_t
     memory_consumption() const override;
 
+    /**
+     * A pointer to the underlying finite element.
+     */
+    SmartPointer<const FiniteElement<dim, spacedim>> fe;
+
     /**
      * Values of shape functions. Access by function @p shape.
      *
@@ -590,7 +600,6 @@ private:
   Point<spacedim>
   do_transform_unit_to_real_cell(const InternalData &mdata) const;
 
-
   /**
    * Transform the point @p p on the real cell to the corresponding point on
    * the unit cell @p cell by a Newton iteration.
@@ -598,10 +607,8 @@ private:
    * Takes a reference to an @p InternalData that is assumed to be previously
    * created by the @p get_data function with @p UpdateFlags including @p
    * update_transformation_values and @p update_transformation_gradients and a
-   * one point Quadrature that includes the given initial guess for the
-   * transformation @p initial_p_unit.  Hence this function assumes that @p
-   * mdata already includes the transformation shape values and gradients
-   * computed at @p initial_p_unit.
+   * one point Quadrature that includes the given initial guess specified
+   * through the given @p point_quadrature object.
    *
    * @p mdata will be changed by this function.
    */
@@ -609,8 +616,8 @@ private:
   do_transform_real_to_unit_cell(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const Point<spacedim>                                      &p,
-    const Point<dim>                                           &initial_p_unit,
-    InternalData                                               &mdata) const;
+    Quadrature<dim> &point_quadrature,
+    InternalData    &mdata) const;
 
   /**
    * Update internal degrees of freedom.
@@ -621,15 +628,6 @@ private:
     const typename MappingFEField<dim, spacedim, VectorType>::InternalData
       &data) const;
 
-  /**
-   * See the documentation of the base class for detailed information.
-   */
-  virtual void
-  compute_shapes_virtual(
-    const std::vector<Point<dim>> &unit_points,
-    typename MappingFEField<dim, spacedim, VectorType>::InternalData &data)
-    const;
-
   /*
    * Which components to use for the mapping.
    */
@@ -658,17 +656,8 @@ private:
   mutable Threads::Mutex fe_values_mutex;
 
   void
-  compute_data(const UpdateFlags      update_flags,
-               const Quadrature<dim> &q,
-               const unsigned int     n_original_q_points,
-               InternalData          &data) const;
-
-  void
-  compute_face_data(const UpdateFlags      update_flags,
-                    const Quadrature<dim> &q,
-                    const unsigned int     n_original_q_points,
-                    InternalData          &data) const;
-
+  compute_face_data(const unsigned int n_original_q_points,
+                    InternalData      &data) const;
 
   // Declare other MappingFEField classes friends.
   template <int, int, class>
diff --git a/include/deal.II/fe/mapping_manifold.h b/include/deal.II/fe/mapping_manifold.h
index db10995ce0..b7de835a49 100644
--- a/include/deal.II/fe/mapping_manifold.h
+++ b/include/deal.II/fe/mapping_manifold.h
@@ -167,18 +167,10 @@ public:
      */
     InternalData() = default;
 
-    /**
-     * Initialize the object's member variables related to cell data based on
-     * the given arguments.
-     *
-     * The function also calls compute_shape_function_values() to actually set
-     * the member variables related to the values and derivatives of the
-     * mapping shape functions.
-     */
-    void
-    initialize(const UpdateFlags      update_flags,
-               const Quadrature<dim> &quadrature,
-               const unsigned int     n_original_q_points);
+    // Documentation see Mapping::InternalDataBase.
+    virtual void
+    reinit(const UpdateFlags      update_flags,
+           const Quadrature<dim> &quadrature) override;
 
     /**
      * Initialize the object's member variables related to cell and face data
@@ -391,37 +383,6 @@ private:
 
 #ifndef DOXYGEN
 
-template <int dim, int spacedim>
-inline void
-MappingManifold<dim, spacedim>::InternalData::store_vertices(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
-{
-  vertices.resize(GeometryInfo<dim>::vertices_per_cell);
-  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-    vertices[i] = cell->vertex(i);
-  this->cell = cell;
-}
-
-
-template <int dim, int spacedim>
-inline void
-MappingManifold<dim, spacedim>::InternalData::
-  compute_manifold_quadrature_weights(const Quadrature<dim> &quad)
-{
-  cell_manifold_quadrature_weights.resize(
-    quad.size(), std::vector<double>(GeometryInfo<dim>::vertices_per_cell));
-  for (unsigned int q = 0; q < quad.size(); ++q)
-    {
-      for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
-        {
-          cell_manifold_quadrature_weights[q][i] =
-            GeometryInfo<dim>::d_linear_shape_function(quad.point(q), i);
-        }
-    }
-}
-
-
-
 template <int dim, int spacedim>
 inline bool
 MappingManifold<dim, spacedim>::preserves_vertex_locations() const
diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h
index e2da64c35c..a440434ae6 100644
--- a/include/deal.II/fe/mapping_q.h
+++ b/include/deal.II/fe/mapping_q.h
@@ -309,14 +309,10 @@ public:
      */
     InternalData(const unsigned int polynomial_degree);
 
-    /**
-     * Initialize the object's member variables related to cell data based on
-     * the given arguments.
-     */
-    void
-    initialize(const UpdateFlags      update_flags,
-               const Quadrature<dim> &quadrature,
-               const unsigned int     n_original_q_points);
+    // Documentation see Mapping::InternalDataBase.
+    virtual void
+    reinit(const UpdateFlags      update_flags,
+           const Quadrature<dim> &quadrature) override;
 
     /**
      * Initialize the object's member variables related to cell and face data
diff --git a/include/deal.II/non_matching/mapping_info.h b/include/deal.II/non_matching/mapping_info.h
index e931f59a1a..93672de532 100644
--- a/include/deal.II/non_matching/mapping_info.h
+++ b/include/deal.II/non_matching/mapping_info.h
@@ -70,24 +70,7 @@ namespace NonMatching
         MappingData &mapping_data)
       {
         mapping_data.initialize(quadrature.size(), update_flags_mapping);
-
-        // reuse internal_mapping_data for MappingQ to avoid memory allocations
-        if (const MappingQ<dim, spacedim> *mapping_q =
-              dynamic_cast<const MappingQ<dim, spacedim> *>(mapping.get()))
-          {
-            (void)mapping_q;
-            auto &data =
-              static_cast<typename MappingQ<dim, spacedim>::InternalData &>(
-                *internal_mapping_data);
-            data.initialize(update_flags_mapping,
-                            quadrature,
-                            quadrature.size());
-          }
-        else
-          {
-            internal_mapping_data =
-              mapping->get_data(update_flags_mapping, quadrature);
-          }
+        internal_mapping_data->reinit(update_flags_mapping, quadrature);
 
         cell_similarity = mapping->fill_fe_values(cell,
                                                   cell_similarity,
@@ -110,21 +93,7 @@ namespace NonMatching
       {
         mapping_data.initialize(quadrature.size(), update_flags_mapping);
 
-        // reuse internal_mapping_data for MappingQ to avoid memory allocations
-        if (dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
-          {
-            auto &data =
-              dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
-                *internal_mapping_data);
-            data.initialize(update_flags_mapping,
-                            quadrature,
-                            quadrature.size());
-          }
-        else
-          {
-            internal_mapping_data =
-              mapping->get_data(update_flags_mapping, quadrature);
-          }
+        internal_mapping_data->reinit(update_flags_mapping, quadrature);
 
         mapping->fill_fe_immersed_surface_values(cell,
                                                  quadrature,
@@ -241,6 +210,12 @@ namespace NonMatching
    * template argument, e.g. VectorizedArray<double>) provides both mapping data
    * and unit points in vectorized format. The Number template parameter of
    * MappingInfo and FEPointEvaluation has to be identical.
+   *
+   * @note This class cannot be copied or copy-constructed. The intention of
+   * this class is to be available as a single use, e.g. inside
+   * FEPointEvaluation or as a storage container for the whole mesh. Use a
+   * shared or unique pointer of this class in case several objects should be
+   * held.
    */
   template <int dim, int spacedim = dim, typename Number = double>
   class MappingInfo : public Subscriptor
@@ -300,6 +275,17 @@ namespace NonMatching
                 const UpdateFlags             update_flags,
                 const AdditionalData additional_data = AdditionalData());
 
+    /**
+     * Do not allow making copies.
+     */
+    MappingInfo(const MappingInfo &) = delete;
+
+    /**
+     * Do not allow copy assignment of this class.
+     */
+    MappingInfo &
+    operator=(const MappingInfo &) = delete;
+
     /**
      * Compute the mapping information for the incoming cell and unit
      * points. This overload is needed to resolve ambiguity.
@@ -665,7 +651,9 @@ namespace NonMatching
     Quadrature<dim> quadrature;
 
     /**
-     * An array containing the data fields of this class.
+     * Helper class that temporarily holds the data requested for one cell by
+     * Mapping::fill_fe_values() before it is filled into appropriate data
+     * structures for consumption by FEPointEvaluation.
      */
     MappingData mapping_data;
 
@@ -837,15 +825,9 @@ namespace NonMatching
       internal::ComputeMappingDataHelper<dim, spacedim>::required_update_flags(
         this->mapping, update_flags_mapping);
 
-    // construct internal_mapping_data for MappingQ to be able to reuse it in
-    // reinit() calls to avoid memory allocations
-    if (const MappingQ<dim, spacedim> *mapping_q =
-          dynamic_cast<const MappingQ<dim, spacedim> *>(&mapping))
-      {
-        internal_mapping_data =
-          std::make_unique<typename MappingQ<dim, spacedim>::InternalData>(
-            mapping_q->get_degree());
-      }
+    // construct internal_mapping_data for mappings for reuse in reinit()
+    // calls to avoid frequent memory allocations
+    internal_mapping_data = mapping.get_data(update_flags, Quadrature<dim>());
   }
 
 
diff --git a/source/fe/fe_hermite.cc b/source/fe/fe_hermite.cc
index 47687f18ed..5e6d97c3ec 100644
--- a/source/fe/fe_hermite.cc
+++ b/source/fe/fe_hermite.cc
@@ -247,8 +247,7 @@ namespace internal
       const typename Mapping<1, spacedim>::InternalDataBase &mapping_data,
       Table<2, Number>                                      &value_list)
     {
-      unsigned int n_q_points;
-      double       cell_extent = 1.0;
+      double cell_extent = 1.0;
 
       // Check mapping_data is associated with a compatible mapping class
       if (dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
@@ -257,7 +256,6 @@ namespace internal
           const typename MappingCartesian<1>::InternalData *data =
             dynamic_cast<const typename MappingCartesian<1>::InternalData *>(
               &mapping_data);
-          n_q_points  = data->quadrature_points.size();
           cell_extent = data->cell_extents[0];
         }
       else
@@ -270,8 +268,6 @@ namespace internal
 
       AssertDimension(value_list.size(0), n_dofs_per_cell);
       AssertDimension(n_dofs_per_cell, 2 * regularity + 2);
-      AssertIndexRange(n_q_points_out, n_q_points + 1);
-      (void)n_q_points;
 
       std::vector<unsigned int> l2h =
         hermite_lexicographic_to_hierarchic_numbering<1>(regularity);
@@ -306,7 +302,6 @@ namespace internal
       const typename Mapping<2, spacedim>::InternalDataBase &mapping_data,
       Table<2, Number>                                      &value_list)
     {
-      unsigned int n_q_points;
       Tensor<1, 2> cell_extents;
 
       // Check mapping_data is associated with a compatible mapping class
@@ -316,7 +311,6 @@ namespace internal
           const typename MappingCartesian<2>::InternalData *data =
             dynamic_cast<const typename MappingCartesian<2>::InternalData *>(
               &mapping_data);
-          n_q_points   = data->quadrature_points.size();
           cell_extents = data->cell_extents;
         }
       else
@@ -330,8 +324,6 @@ namespace internal
 
       AssertDimension(value_list.size(0), n_dofs_per_cell);
       AssertDimension(n_dofs_per_dim * n_dofs_per_dim, n_dofs_per_cell);
-      AssertIndexRange(n_q_points_out, n_q_points + 1);
-      (void)n_q_points;
 
       std::vector<unsigned int> l2h =
         hermite_lexicographic_to_hierarchic_numbering<2>(regularity);
@@ -380,7 +372,6 @@ namespace internal
       const typename Mapping<3, spacedim>::InternalDataBase &mapping_data,
       Table<2, Number>                                      &value_list)
     {
-      unsigned int n_q_points;
       Tensor<1, 3> cell_extents;
 
       // Check mapping_data is associated with a compatible mapping class
@@ -390,7 +381,6 @@ namespace internal
           const typename MappingCartesian<3>::InternalData *data =
             dynamic_cast<const typename MappingCartesian<3>::InternalData *>(
               &mapping_data);
-          n_q_points   = data->quadrature_points.size();
           cell_extents = data->cell_extents;
         }
       else
@@ -405,8 +395,6 @@ namespace internal
 
       AssertDimension(value_list.size(0), n_dofs_per_cell);
       AssertDimension(Utilities::pow(n_dofs_per_dim, 3), n_dofs_per_cell);
-      AssertIndexRange(n_q_points_out, n_q_points + 1);
-      (void)n_q_points;
 
       std::vector<unsigned int> l2h =
         hermite_lexicographic_to_hierarchic_numbering<3>(regularity);
diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc
index 4487a54854..314a2d18a8 100644
--- a/source/fe/mapping.cc
+++ b/source/fe/mapping.cc
@@ -269,6 +269,16 @@ Mapping<dim, spacedim>::InternalDataBase::InternalDataBase()
 
 
 
+template <int dim, int spacedim>
+void
+Mapping<dim, spacedim>::InternalDataBase::reinit(const UpdateFlags,
+                                                 const Quadrature<dim> &)
+{
+  DEAL_II_ASSERT_UNREACHABLE();
+}
+
+
+
 template <int dim, int spacedim>
 std::size_t
 Mapping<dim, spacedim>::InternalDataBase::memory_consumption() const
diff --git a/source/fe/mapping_cartesian.cc b/source/fe/mapping_cartesian.cc
index 07b185eddf..2bb5977b7e 100644
--- a/source/fe/mapping_cartesian.cc
+++ b/source/fe/mapping_cartesian.cc
@@ -98,6 +98,20 @@ MappingCartesian<dim, spacedim>::InternalData::InternalData(
 
 
 
+template <int dim, int spacedim>
+void
+MappingCartesian<dim, spacedim>::InternalData::reinit(
+  const UpdateFlags update_flags,
+  const Quadrature<dim> &)
+{
+  // store the flags in the internal data object so we can access them
+  // in fill_fe_*_values(). use the transitive hull of the required
+  // flags
+  this->update_each = update_flags;
+}
+
+
+
 template <int dim, int spacedim>
 std::size_t
 MappingCartesian<dim, spacedim>::InternalData::memory_consumption() const
@@ -105,8 +119,7 @@ MappingCartesian<dim, spacedim>::InternalData::memory_consumption() const
   return (Mapping<dim, spacedim>::InternalDataBase::memory_consumption() +
           MemoryConsumption::memory_consumption(cell_extents) +
           MemoryConsumption::memory_consumption(inverse_cell_extents) +
-          MemoryConsumption::memory_consumption(volume_element) +
-          MemoryConsumption::memory_consumption(quadrature_points));
+          MemoryConsumption::memory_consumption(volume_element));
 }
 
 
@@ -162,13 +175,8 @@ MappingCartesian<dim, spacedim>::get_data(const UpdateFlags      update_flags,
                                           const Quadrature<dim> &q) const
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
-    std::make_unique<InternalData>(q);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-
-  // store the flags in the internal data object so we can access them
-  // in fill_fe_*_values(). use the transitive hull of the required
-  // flags
-  data.update_each = requires_update_flags(update_flags);
+    std::make_unique<InternalData>();
+  data_ptr->reinit(requires_update_flags(update_flags), q);
 
   return data_ptr;
 }
@@ -252,18 +260,46 @@ MappingCartesian<dim, spacedim>::update_cell_extents(
 
 
 
+namespace
+{
+  template <int dim>
+  void
+  transform_quadrature_points(
+    const Tensor<1, dim>                               first_vertex,
+    const Tensor<1, dim>                               cell_extents,
+    const ArrayView<const Point<dim>>                 &unit_quadrature_points,
+    const typename QProjector<dim>::DataSetDescriptor &offset,
+    std::vector<Point<dim>>                           &quadrature_points)
+  {
+    for (unsigned int i = 0; i < quadrature_points.size(); ++i)
+      {
+        quadrature_points[i] = first_vertex;
+        for (unsigned int d = 0; d < dim; ++d)
+          quadrature_points[i][d] +=
+            cell_extents[d] * unit_quadrature_points[i + offset][d];
+      }
+  }
+} // namespace
+
+
+
 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
+  const ArrayView<const Point<dim>> &unit_quadrature_points,
+  std::vector<Point<dim>>           &quadrature_points) const
 {
   if (data.update_each & update_quadrature_points)
     {
       const auto offset = QProjector<dim>::DataSetDescriptor::cell();
 
-      transform_quadrature_points(cell, data, offset, quadrature_points);
+      transform_quadrature_points(cell->vertex(0),
+                                  data.cell_extents,
+                                  unit_quadrature_points,
+                                  offset,
+                                  quadrature_points);
     }
 }
 
@@ -288,7 +324,11 @@ MappingCartesian<dim, spacedim>::maybe_update_face_quadrature_points(
         quadrature_points.size());
 
 
-      transform_quadrature_points(cell, data, offset, quadrature_points);
+      transform_quadrature_points(cell->vertex(0),
+                                  data.cell_extents,
+                                  make_array_view(data.quadrature_points),
+                                  offset,
+                                  quadrature_points);
     }
 }
 
@@ -320,30 +360,11 @@ MappingCartesian<dim, spacedim>::maybe_update_subface_quadrature_points(
         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];
+      transform_quadrature_points(cell->vertex(0),
+                                  data.cell_extents,
+                                  make_array_view(data.quadrature_points),
+                                  offset,
+                                  quadrature_points);
     }
 }
 
@@ -504,6 +525,7 @@ MappingCartesian<dim, spacedim>::fill_fe_values(
 
   maybe_update_cell_quadrature_points(cell,
                                       data,
+                                      quadrature.get_points(),
                                       output_data.quadrature_points);
 
   // compute Jacobian determinant. all values are equal and are the
@@ -553,13 +575,12 @@ MappingCartesian<dim, spacedim>::fill_mapping_data_for_generic_points(
 
   InternalData data;
   data.update_each = update_flags;
-  data.quadrature_points =
-    std::vector<Point<dim>>(unit_points.begin(), unit_points.end());
 
   update_cell_extents(cell, CellSimilarity::none, data);
 
   maybe_update_cell_quadrature_points(cell,
                                       data,
+                                      unit_points,
                                       output_data.quadrature_points);
 
   maybe_update_jacobians(data, CellSimilarity::none, output_data);
@@ -703,6 +724,7 @@ MappingCartesian<dim, spacedim>::fill_fe_immersed_surface_values(
 
   maybe_update_cell_quadrature_points(cell,
                                       data,
+                                      quadrature.get_points(),
                                       output_data.quadrature_points);
 
   if (data.update_each & update_normal_vectors)
diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc
index 9707bb7a81..7856db406f 100644
--- a/source/fe/mapping_fe.cc
+++ b/source/fe/mapping_fe.cc
@@ -74,12 +74,11 @@ MappingFE<dim, spacedim>::InternalData::memory_consumption() const
 }
 
 
+
 template <int dim, int spacedim>
 void
-MappingFE<dim, spacedim>::InternalData::initialize(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
-  const unsigned int     n_original_q_points)
+MappingFE<dim, spacedim>::InternalData::reinit(const UpdateFlags update_flags,
+                                               const Quadrature<dim> &q)
 {
   // store the flags in the internal data object so we can access them
   // in fill_fe_*_values()
@@ -88,13 +87,13 @@ MappingFE<dim, spacedim>::InternalData::initialize(
   const unsigned int n_q_points = q.size();
 
   if (this->update_each & update_covariant_transformation)
-    covariant.resize(n_original_q_points);
+    covariant.resize(n_q_points);
 
   if (this->update_each & update_contravariant_transformation)
-    contravariant.resize(n_original_q_points);
+    contravariant.resize(n_q_points);
 
   if (this->update_each & update_volume_elements)
-    volume_elements.resize(n_original_q_points);
+    volume_elements.resize(n_q_points);
 
   // see if we need the (transformation) shape function values
   // and/or gradients and resize the necessary arrays
@@ -139,7 +138,7 @@ MappingFE<dim, spacedim>::InternalData::initialize_face(
   const Quadrature<dim> &q,
   const unsigned int     n_original_q_points)
 {
-  initialize(update_flags, q, n_original_q_points);
+  reinit(update_flags, q);
 
   if (this->update_each &
       (update_boundary_forms | update_normal_vectors | update_jacobians |
@@ -1059,9 +1058,7 @@ MappingFE<dim, spacedim>::get_data(const UpdateFlags      update_flags,
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>(*this->fe);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-  data.initialize(this->requires_update_flags(update_flags), q, q.size());
-
+  data_ptr->reinit(this->requires_update_flags(update_flags), q);
   return data_ptr;
 }
 
diff --git a/source/fe/mapping_fe_field.cc b/source/fe/mapping_fe_field.cc
index 5fa44c8548..50fd675c30 100644
--- a/source/fe/mapping_fe_field.cc
+++ b/source/fe/mapping_fe_field.cc
@@ -59,7 +59,8 @@ template <int dim, int spacedim, typename VectorType>
 MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
   const FiniteElement<dim, spacedim> &fe,
   const ComponentMask                &mask)
-  : unit_tangentials()
+  : fe(&fe)
+  , unit_tangentials()
   , n_shape_functions(fe.n_dofs_per_cell())
   , mask(mask)
   , local_dof_indices(fe.n_dofs_per_cell())
@@ -68,6 +69,86 @@ MappingFEField<dim, spacedim, VectorType>::InternalData::InternalData(
 
 
 
+template <int dim, int spacedim, typename VectorType>
+void
+MappingFEField<dim, spacedim, VectorType>::InternalData::reinit(
+  const UpdateFlags      update_flags,
+  const Quadrature<dim> &quadrature)
+{
+  // store the flags in the internal data object so we can access them
+  // in fill_fe_*_values(). use the transitive hull of the required
+  // flags
+  this->update_each = update_flags;
+
+  const unsigned int             n_q_points = quadrature.size();
+  const std::vector<Point<dim>> &points     = quadrature.get_points();
+
+  // see if we need the (transformation) shape function values
+  // and/or gradients and resize the necessary arrays
+  if (update_flags & update_quadrature_points)
+    {
+      shape_values.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          shape(point, i) = fe->shape_value(i, points[point]);
+    }
+
+  if (update_flags &
+      (update_covariant_transformation | update_contravariant_transformation |
+       update_JxW_values | update_boundary_forms | update_normal_vectors |
+       update_jacobians | update_jacobian_grads | update_inverse_jacobians))
+    {
+      shape_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          derivative(point, i) = fe->shape_grad(i, points[point]);
+    }
+
+  if (update_flags & update_covariant_transformation)
+    covariant.resize(n_q_points);
+
+  if (update_flags & update_contravariant_transformation)
+    contravariant.resize(n_q_points);
+
+  if (update_flags & update_volume_elements)
+    volume_elements.resize(n_q_points);
+
+  if (update_flags &
+      (update_jacobian_grads | update_jacobian_pushed_forward_grads))
+    {
+      shape_second_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          second_derivative(point, i) = fe->shape_grad_grad(i, points[point]);
+    }
+
+  if (update_flags & (update_jacobian_2nd_derivatives |
+                      update_jacobian_pushed_forward_2nd_derivatives))
+    {
+      shape_third_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          third_derivative(point, i) =
+            fe->shape_3rd_derivative(i, points[point]);
+    }
+
+  if (update_flags & (update_jacobian_3rd_derivatives |
+                      update_jacobian_pushed_forward_3rd_derivatives))
+    {
+      shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
+      for (unsigned int point = 0; point < n_q_points; ++point)
+        for (unsigned int i = 0; i < n_shape_functions; ++i)
+          fourth_derivative(point, i) =
+            fe->shape_4th_derivative(i, points[point]);
+    }
+
+  // This (for face values and simplices) can be different for different
+  // calls, so always copy
+  quadrature_weights = quadrature.get_weights();
+}
+
+
+
 template <int dim, int spacedim, typename VectorType>
 std::size_t
 MappingFEField<dim, spacedim, VectorType>::InternalData::memory_consumption()
@@ -420,43 +501,6 @@ MappingFEField<dim, spacedim, VectorType>::get_vertices(
 
 
 
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::compute_shapes_virtual(
-  const std::vector<Point<dim>>                                    &unit_points,
-  typename MappingFEField<dim, spacedim, VectorType>::InternalData &data) const
-{
-  const auto         fe       = &euler_dof_handler->get_fe();
-  const unsigned int n_points = unit_points.size();
-
-  for (unsigned int point = 0; point < n_points; ++point)
-    {
-      if (data.shape_values.size() != 0)
-        for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-          data.shape(point, i) = fe->shape_value(i, unit_points[point]);
-
-      if (data.shape_derivatives.size() != 0)
-        for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-          data.derivative(point, i) = fe->shape_grad(i, unit_points[point]);
-
-      if (data.shape_second_derivatives.size() != 0)
-        for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-          data.second_derivative(point, i) =
-            fe->shape_grad_grad(i, unit_points[point]);
-
-      if (data.shape_third_derivatives.size() != 0)
-        for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-          data.third_derivative(point, i) =
-            fe->shape_3rd_derivative(i, unit_points[point]);
-
-      if (data.shape_fourth_derivatives.size() != 0)
-        for (unsigned int i = 0; i < data.n_shape_functions; ++i)
-          data.fourth_derivative(point, i) =
-            fe->shape_4th_derivative(i, unit_points[point]);
-    }
-}
-
-
 template <int dim, int spacedim, typename VectorType>
 UpdateFlags
 MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
@@ -511,33 +555,14 @@ MappingFEField<dim, spacedim, VectorType>::requires_update_flags(
 }
 
 
-
 template <int dim, int spacedim, typename VectorType>
 void
-MappingFEField<dim, spacedim, VectorType>::compute_data(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
-  const unsigned int     n_original_q_points,
-  InternalData          &data) const
+MappingFEField<dim, spacedim, VectorType>::compute_face_data(
+  const unsigned int n_original_q_points,
+  InternalData      &data) const
 {
-  // store the flags in the internal data object so we can access them
-  // in fill_fe_*_values(). use the transitive hull of the required
-  // flags
-  data.update_each = requires_update_flags(update_flags);
-
-  const unsigned int n_q_points = q.size();
-
-  // see if we need the (transformation) shape function values
-  // and/or gradients and resize the necessary arrays
-  if (data.update_each & update_quadrature_points)
-    data.shape_values.resize(data.n_shape_functions * n_q_points);
-
-  if (data.update_each &
-      (update_covariant_transformation | update_contravariant_transformation |
-       update_JxW_values | update_boundary_forms | update_normal_vectors |
-       update_jacobians | update_jacobian_grads | update_inverse_jacobians))
-    data.shape_derivatives.resize(data.n_shape_functions * n_q_points);
-
+  // Set to the size of a single quadrature object for faces, as the size set
+  // in in reinit() is for all points
   if (data.update_each & update_covariant_transformation)
     data.covariant.resize(n_original_q_points);
 
@@ -547,36 +572,6 @@ MappingFEField<dim, spacedim, VectorType>::compute_data(
   if (data.update_each & update_volume_elements)
     data.volume_elements.resize(n_original_q_points);
 
-  if (data.update_each &
-      (update_jacobian_grads | update_jacobian_pushed_forward_grads))
-    data.shape_second_derivatives.resize(data.n_shape_functions * n_q_points);
-
-  if (data.update_each & (update_jacobian_2nd_derivatives |
-                          update_jacobian_pushed_forward_2nd_derivatives))
-    data.shape_third_derivatives.resize(data.n_shape_functions * n_q_points);
-
-  if (data.update_each & (update_jacobian_3rd_derivatives |
-                          update_jacobian_pushed_forward_3rd_derivatives))
-    data.shape_fourth_derivatives.resize(data.n_shape_functions * n_q_points);
-
-  compute_shapes_virtual(q.get_points(), data);
-
-  // This (for face values and simplices) can be different for different calls,
-  // so always copy
-  data.quadrature_weights = q.get_weights();
-}
-
-
-template <int dim, int spacedim, typename VectorType>
-void
-MappingFEField<dim, spacedim, VectorType>::compute_face_data(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
-  const unsigned int     n_original_q_points,
-  InternalData          &data) const
-{
-  compute_data(update_flags, q, n_original_q_points, data);
-
   if (dim > 1)
     {
       if (data.update_each & update_boundary_forms)
@@ -611,6 +606,7 @@ MappingFEField<dim, spacedim, VectorType>::compute_face_data(
 }
 
 
+
 template <int dim, int spacedim, typename VectorType>
 typename std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
 MappingFEField<dim, spacedim, VectorType>::get_data(
@@ -619,8 +615,7 @@ MappingFEField<dim, spacedim, VectorType>::get_data(
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-  this->compute_data(update_flags, quadrature, quadrature.size(), data);
+  data_ptr->reinit(requires_update_flags(update_flags), quadrature);
 
   return data_ptr;
 }
@@ -637,10 +632,12 @@ MappingFEField<dim, spacedim, VectorType>::get_face_data(
 
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  auto                 &data = dynamic_cast<InternalData &>(*data_ptr);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
   const Quadrature<dim> q(
     QProjector<dim>::project_to_all_faces(reference_cell, quadrature[0]));
-  this->compute_face_data(update_flags, q, quadrature[0].size(), data);
+  data.reinit(requires_update_flags(update_flags), q);
+  this->compute_face_data(quadrature[0].size(), data);
 
   return data_ptr;
 }
@@ -654,10 +651,12 @@ MappingFEField<dim, spacedim, VectorType>::get_subface_data(
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>(euler_dof_handler->get_fe(), fe_mask);
-  auto                 &data = dynamic_cast<InternalData &>(*data_ptr);
+  auto &data = dynamic_cast<InternalData &>(*data_ptr);
+
   const Quadrature<dim> q(
     QProjector<dim>::project_to_all_subfaces(reference_cell, quadrature));
-  this->compute_face_data(update_flags, q, quadrature.size(), data);
+  data.reinit(requires_update_flags(update_flags), q);
+  this->compute_face_data(quadrature.size(), data);
 
   return data_ptr;
 }
@@ -2198,9 +2197,9 @@ MappingFEField<dim, spacedim, VectorType>::transform_unit_to_real_cell(
   Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
          ExcInternalError());
 
-  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
+  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
 
-  return do_transform_unit_to_real_cell(dynamic_cast<InternalData &>(*mdata));
+  return do_transform_unit_to_real_cell(static_cast<InternalData &>(*mdata));
 }
 
 
@@ -2250,7 +2249,7 @@ MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
 
   initial_p_unit = cell->reference_cell().closest_point(initial_p_unit);
 
-  const Quadrature<dim> point_quadrature(initial_p_unit);
+  Quadrature<dim> point_quadrature(initial_p_unit);
 
   UpdateFlags update_flags = update_quadrature_points | update_jacobians;
   if (spacedim > dim)
@@ -2260,12 +2259,12 @@ MappingFEField<dim, spacedim, VectorType>::transform_real_to_unit_cell(
   Assert(dynamic_cast<InternalData *>(mdata.get()) != nullptr,
          ExcInternalError());
 
-  update_internal_dofs(cell, dynamic_cast<InternalData &>(*mdata));
+  update_internal_dofs(cell, static_cast<InternalData &>(*mdata));
 
   return do_transform_real_to_unit_cell(cell,
                                         p,
-                                        initial_p_unit,
-                                        dynamic_cast<InternalData &>(*mdata));
+                                        point_quadrature,
+                                        static_cast<InternalData &>(*mdata));
 }
 
 
@@ -2274,7 +2273,7 @@ Point<dim>
 MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const Point<spacedim>                                      &p,
-  const Point<dim>                                           &initial_p_unit,
+  Quadrature<dim>                                            &point_quadrature,
   InternalData                                               &mdata) const
 {
   const unsigned int n_shapes = mdata.shape_values.size();
@@ -2291,10 +2290,12 @@ MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
   // The shape values and derivatives
   // of the mapping at this point are
   // previously computed.
-  // f(x)
-  Point<dim> p_unit = initial_p_unit;
+
+  AssertDimension(point_quadrature.size(), 1);
+  Point<dim> p_unit = point_quadrature.point(0);
   Point<dim> f;
-  compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit), mdata);
+  mdata.reinit(mdata.update_each, point_quadrature);
+
   Point<spacedim>     p_real(do_transform_unit_to_real_cell(mdata));
   Tensor<1, spacedim> p_minus_F              = p - p_real;
   const double        eps                    = 1.e-12 * cell->diameter();
@@ -2342,8 +2343,9 @@ MappingFEField<dim, spacedim, VectorType>::do_transform_real_to_unit_cell(
             p_unit_trial[i] -= step_length * delta[i];
           // shape values and derivatives
           // at new p_unit point
-          compute_shapes_virtual(std::vector<Point<dim>>(1, p_unit_trial),
-                                 mdata);
+          point_quadrature.initialize(
+            ArrayView<const Point<dim>>(&p_unit_trial, 1));
+          mdata.reinit(mdata.update_each, point_quadrature);
           // f(x)
           Point<spacedim> p_real_trial = do_transform_unit_to_real_cell(mdata);
           const Tensor<1, spacedim> f_trial = p - p_real_trial;
diff --git a/source/fe/mapping_manifold.cc b/source/fe/mapping_manifold.cc
index 5ae5ffb3e6..551b83b87d 100644
--- a/source/fe/mapping_manifold.cc
+++ b/source/fe/mapping_manifold.cc
@@ -63,19 +63,21 @@ MappingManifold<dim, spacedim>::InternalData::memory_consumption() const
 }
 
 
+
 template <int dim, int spacedim>
 void
-MappingManifold<dim, spacedim>::InternalData::initialize(
+MappingManifold<dim, spacedim>::InternalData::reinit(
   const UpdateFlags      update_flags,
-  const Quadrature<dim> &q,
-  const unsigned int     n_original_q_points)
+  const Quadrature<dim> &q)
 {
   // store the flags in the internal data object so we can access them
   // in fill_fe_*_values()
   this->update_each = update_flags;
 
+  const unsigned int n_q_points = q.size();
+
   // Store the quadrature
-  this->quad = q;
+  this->quad.initialize(q.get_points(), q.get_weights());
 
   // Resize the weights
   this->vertex_weights.resize(GeometryInfo<dim>::vertices_per_cell);
@@ -87,16 +89,17 @@ MappingManifold<dim, spacedim>::InternalData::initialize(
     compute_manifold_quadrature_weights(q);
 
   if (this->update_each & update_covariant_transformation)
-    covariant.resize(n_original_q_points);
+    covariant.resize(n_q_points);
 
   if (this->update_each & update_contravariant_transformation)
-    contravariant.resize(n_original_q_points);
+    contravariant.resize(n_q_points);
 
   if (this->update_each & update_volume_elements)
-    volume_elements.resize(n_original_q_points);
+    volume_elements.resize(n_q_points);
 }
 
 
+
 template <int dim, int spacedim>
 void
 MappingManifold<dim, spacedim>::InternalData::initialize_face(
@@ -104,7 +107,18 @@ MappingManifold<dim, spacedim>::InternalData::initialize_face(
   const Quadrature<dim> &q,
   const unsigned int     n_original_q_points)
 {
-  initialize(update_flags, q, n_original_q_points);
+  reinit(update_flags, q);
+
+  // Set to the size of a single quadrature object for faces, as the size set
+  // in in reinit() is for all points
+  if (this->update_each & update_covariant_transformation)
+    covariant.resize(n_original_q_points);
+
+  if (this->update_each & update_contravariant_transformation)
+    contravariant.resize(n_original_q_points);
+
+  if (this->update_each & update_volume_elements)
+    volume_elements.resize(n_original_q_points);
 
   if (dim > 1)
     {
@@ -138,6 +152,38 @@ MappingManifold<dim, spacedim>::InternalData::initialize_face(
 
 
 
+template <int dim, int spacedim>
+void
+MappingManifold<dim, spacedim>::InternalData::store_vertices(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
+{
+  vertices.resize(GeometryInfo<dim>::vertices_per_cell);
+  for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+    vertices[i] = cell->vertex(i);
+  this->cell = cell;
+}
+
+
+
+template <int dim, int spacedim>
+void
+MappingManifold<dim, spacedim>::InternalData::
+  compute_manifold_quadrature_weights(const Quadrature<dim> &quad)
+{
+  cell_manifold_quadrature_weights.resize(
+    quad.size(), std::vector<double>(GeometryInfo<dim>::vertices_per_cell));
+  for (unsigned int q = 0; q < quad.size(); ++q)
+    {
+      for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
+        {
+          cell_manifold_quadrature_weights[q][i] =
+            GeometryInfo<dim>::d_linear_shape_function(quad.point(q), i);
+        }
+    }
+}
+
+
+
 template <int dim, int spacedim>
 MappingManifold<dim, spacedim>::MappingManifold(
   const MappingManifold<dim, spacedim> &)
@@ -269,8 +315,7 @@ MappingManifold<dim, spacedim>::get_data(const UpdateFlags      update_flags,
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>();
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-  data.initialize(this->requires_update_flags(update_flags), q, q.size());
+  data_ptr->reinit(this->requires_update_flags(update_flags), q);
 
   return data_ptr;
 }
diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc
index 5b77e81603..c51f40bf7e 100644
--- a/source/fe/mapping_q.cc
+++ b/source/fe/mapping_q.cc
@@ -78,10 +78,8 @@ MappingQ<dim, spacedim>::InternalData::memory_consumption() const
 
 template <int dim, int spacedim>
 void
-MappingQ<dim, spacedim>::InternalData::initialize(
-  const UpdateFlags      update_flags,
-  const Quadrature<dim> &quadrature,
-  const unsigned int     n_original_q_points)
+MappingQ<dim, spacedim>::InternalData::reinit(const UpdateFlags update_flags,
+                                              const Quadrature<dim> &quadrature)
 {
   // store the flags in the internal data object so we can access them
   // in fill_fe_*_values()
@@ -90,7 +88,7 @@ MappingQ<dim, spacedim>::InternalData::initialize(
   const unsigned int n_q_points = quadrature.size();
 
   if (this->update_each & update_volume_elements)
-    volume_elements.resize(n_original_q_points);
+    volume_elements.resize(n_q_points);
 
   tensor_product_quadrature = quadrature.is_tensor_product();
 
@@ -163,7 +161,7 @@ MappingQ<dim, spacedim>::InternalData::initialize_face(
   const Quadrature<dim> &quadrature,
   const unsigned int     n_original_q_points)
 {
-  initialize(update_flags, quadrature, n_original_q_points);
+  reinit(update_flags, quadrature);
 
   quadrature_points = quadrature.get_points();
 
@@ -785,9 +783,7 @@ MappingQ<dim, spacedim>::get_data(const UpdateFlags      update_flags,
 {
   std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase> data_ptr =
     std::make_unique<InternalData>(polynomial_degree);
-  auto &data = dynamic_cast<InternalData &>(*data_ptr);
-  data.initialize(this->requires_update_flags(update_flags), q, q.size());
-
+  data_ptr->reinit(update_flags, q);
   return data_ptr;
 }