]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Remove tabulated shape functions from MappingQ 15125/head
authorMartin Kronbichler <martin.kronbichler@uni-a.de>
Fri, 21 Apr 2023 16:54:19 +0000 (18:54 +0200)
committerMartin Kronbichler <martin.kronbichler@uni-a.de>
Fri, 21 Apr 2023 21:19:33 +0000 (23:19 +0200)
include/deal.II/fe/fe_values.h
include/deal.II/fe/mapping_q.h
include/deal.II/fe/mapping_q_internal.h
source/fe/mapping_q.cc

index 8b5e7affffe97124b8edc032ce8c1489468b25f0..b3d385f635cd0b27c6198d6c225d7f96882af87c 100644 (file)
@@ -3954,7 +3954,6 @@ protected:
   internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
     mapping_output;
 
-
   /**
    * A pointer to the finite element object associated with this FEValues
    * object.
index 28f34de1bf96fd46865421311a0fa6bf48534b05..0643b17e24dbaff51042c80d3db75c749bcfe29f 100644 (file)
@@ -334,88 +334,6 @@ public:
                     const Quadrature<dim> &quadrature,
                     const unsigned int     n_original_q_points);
 
-    /**
-     * Compute the values and/or derivatives of the shape functions used for
-     * the mapping.
-     *
-     * Which values, derivatives, or higher order derivatives are computed is
-     * determined by which of the member arrays have nonzero sizes. They are
-     * typically set to their appropriate sizes by the initialize() and
-     * initialize_face() functions, which indeed call this function
-     * internally. However, it is possible (and at times useful) to do the
-     * resizing by hand and then call this function directly. An example is in
-     * a Newton iteration where we update the location of a quadrature point
-     * (e.g., in MappingQ::transform_real_to_uni_cell()) and need to re-
-     * compute the mapping and its derivatives at this location, but have
-     * already sized all internal arrays correctly.
-     */
-    void
-    compute_shape_function_values(const std::vector<Point<dim>> &unit_points);
-
-    /**
-     * Shape function at quadrature point. Shape functions are in tensor
-     * product order, so vertices must be reordered to obtain transformation.
-     */
-    const double &
-    shape(const unsigned int qpoint, const unsigned int shape_nr) const;
-
-    /**
-     * Shape function at quadrature point. See above.
-     */
-    double &
-    shape(const unsigned int qpoint, const unsigned int shape_nr);
-
-    /**
-     * Gradient of shape function in quadrature point. See above.
-     */
-    const Tensor<1, dim> &
-    derivative(const unsigned int qpoint, const unsigned int shape_nr) const;
-
-    /**
-     * Gradient of shape function in quadrature point. See above.
-     */
-    Tensor<1, dim> &
-    derivative(const unsigned int qpoint, const unsigned int shape_nr);
-
-    /**
-     * Second derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<2, dim> &
-    second_derivative(const unsigned int qpoint,
-                      const unsigned int shape_nr) const;
-
-    /**
-     * Second derivative of shape function in quadrature point. See above.
-     */
-    Tensor<2, dim> &
-    second_derivative(const unsigned int qpoint, const unsigned int shape_nr);
-
-    /**
-     * third derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<3, dim> &
-    third_derivative(const unsigned int qpoint,
-                     const unsigned int shape_nr) const;
-
-    /**
-     * third derivative of shape function in quadrature point. See above.
-     */
-    Tensor<3, dim> &
-    third_derivative(const unsigned int qpoint, const unsigned int shape_nr);
-
-    /**
-     * fourth derivative of shape function in quadrature point. See above.
-     */
-    const Tensor<4, dim> &
-    fourth_derivative(const unsigned int qpoint,
-                      const unsigned int shape_nr) const;
-
-    /**
-     * fourth derivative of shape function in quadrature point. See above.
-     */
-    Tensor<4, dim> &
-    fourth_derivative(const unsigned int qpoint, const unsigned int shape_nr);
-
     /**
      * Return an estimate (in bytes) for the memory consumption of this object.
      */
@@ -427,45 +345,7 @@ public:
      * possible orientations. Can be accessed with the correct offset provided
      * via QProjector::DataSetDescriptor. Not needed/used for cells.
      */
-    AlignedVector<Point<dim>> quadrature_points;
-
-    /**
-     * Values of shape functions. Access by function @p shape.
-     *
-     * Computed once.
-     */
-    AlignedVector<double> shape_values;
-
-    /**
-     * Values of shape function derivatives. Access by function @p derivative.
-     *
-     * Computed once.
-     */
-    AlignedVector<Tensor<1, dim>> shape_derivatives;
-
-    /**
-     * Values of shape function second derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    AlignedVector<Tensor<2, dim>> shape_second_derivatives;
-
-    /**
-     * Values of shape function third derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    AlignedVector<Tensor<3, dim>> shape_third_derivatives;
-
-    /**
-     * Values of shape function fourth derivatives. Access by function @p
-     * second_derivative.
-     *
-     * Computed once.
-     */
-    AlignedVector<Tensor<4, dim>> shape_fourth_derivatives;
+    std::vector<Point<dim>> quadrature_points;
 
     /**
      * Unit tangential vectors. Used for the computation of boundary forms and
@@ -548,26 +428,6 @@ public:
      */
     bool tensor_product_quadrature;
 
-    /**
-     * Tensors of covariant transformation at each of the quadrature points.
-     * The matrix stored is the Jacobian * G^{-1}, where G = Jacobian^{t} *
-     * Jacobian, is the first fundamental form of the map; if dim=spacedim
-     * then it reduces to the transpose of the inverse of the Jacobian matrix,
-     * which itself is stored in the @p contravariant field of this structure.
-     *
-     * Computed on each cell.
-     */
-    mutable AlignedVector<DerivativeForm<1, dim, spacedim>> covariant;
-
-    /**
-     * Tensors of contravariant transformation at each of the quadrature
-     * points. The contravariant matrix is the Jacobian of the transformation,
-     * i.e. $J_{ij}=dx_i/d\hat x_j$.
-     *
-     * Computed on each cell.
-     */
-    mutable AlignedVector<DerivativeForm<1, dim, spacedim>> contravariant;
-
     /**
      * Auxiliary vectors for internal use.
      */
@@ -590,6 +450,14 @@ public:
      * #update_volume_elements.
      */
     mutable AlignedVector<double> volume_elements;
+
+    /**
+     * Pointer to the mapping output data that holds most of the arrays,
+     * including the Jacobians representing the covariant and contravariant
+     * transformations.
+     */
+    mutable internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+      *output_data;
   };
 
 protected:
@@ -840,123 +708,6 @@ using MappingQGeneric = MappingQ<dim, spacedim>;
 
 #ifndef DOXYGEN
 
-template <int dim, int spacedim>
-inline const double &
-MappingQ<dim, spacedim>::InternalData::shape(const unsigned int qpoint,
-                                             const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
-  return shape_values[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim>
-inline double &
-MappingQ<dim, spacedim>::InternalData::shape(const unsigned int qpoint,
-                                             const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr, shape_values.size());
-  return shape_values[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline const Tensor<1, dim> &
-MappingQ<dim, spacedim>::InternalData::derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_derivatives.size());
-  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
-template <int dim, int spacedim>
-inline Tensor<1, dim> &
-MappingQ<dim, spacedim>::InternalData::derivative(const unsigned int qpoint,
-                                                  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_derivatives.size());
-  return shape_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline const Tensor<2, dim> &
-MappingQ<dim, spacedim>::InternalData::second_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_second_derivatives.size());
-  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline Tensor<2, dim> &
-MappingQ<dim, spacedim>::InternalData::second_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_second_derivatives.size());
-  return shape_second_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-template <int dim, int spacedim>
-inline const Tensor<3, dim> &
-MappingQ<dim, spacedim>::InternalData::third_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_third_derivatives.size());
-  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline Tensor<3, dim> &
-MappingQ<dim, spacedim>::InternalData::third_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_third_derivatives.size());
-  return shape_third_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline const Tensor<4, dim> &
-MappingQ<dim, spacedim>::InternalData::fourth_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr) const
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_fourth_derivatives.size());
-  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-template <int dim, int spacedim>
-inline Tensor<4, dim> &
-MappingQ<dim, spacedim>::InternalData::fourth_derivative(
-  const unsigned int qpoint,
-  const unsigned int shape_nr)
-{
-  AssertIndexRange(qpoint * n_shape_functions + shape_nr,
-                   shape_fourth_derivatives.size());
-  return shape_fourth_derivatives[qpoint * n_shape_functions + shape_nr];
-}
-
-
-
 template <int dim, int spacedim>
 inline bool
 MappingQ<dim, spacedim>::preserves_vertex_locations() const
index e144e330966f93a7221b6ae751cebab709e9ec44..e39202634c6a55542ac8439d5805a27eeb30dc6a 100644 (file)
@@ -733,111 +733,53 @@ namespace internal
     template <int dim>
     inline Point<dim>
     do_transform_real_to_unit_cell_internal_codim1(
-      const typename dealii::Triangulation<dim, dim + 1>::cell_iterator &cell,
-      const Point<dim + 1> &                                             p,
-      const Point<dim> &                                     initial_p_unit,
-      typename dealii::MappingQ<dim, dim + 1>::InternalData &mdata)
+      const Point<dim + 1> &                              p,
+      const Point<dim> &                                  initial_p_unit,
+      const std::vector<Point<dim + 1>> &                 points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> &                   renumber)
     {
       const unsigned int spacedim = dim + 1;
 
-      const unsigned int n_shapes = mdata.shape_values.size();
-      (void)n_shapes;
-      Assert(n_shapes != 0, ExcInternalError());
-      Assert(mdata.shape_derivatives.size() == n_shapes, ExcInternalError());
-      Assert(mdata.shape_second_derivatives.size() == n_shapes,
-             ExcInternalError());
-
-      std::vector<Point<spacedim>> &points = mdata.mapping_support_points;
-      Assert(points.size() == n_shapes, ExcInternalError());
-
-      Point<spacedim> p_minus_F;
+      AssertDimension(points.size(),
+                      Utilities::pow(polynomials_1d.size(), dim));
 
-      Tensor<1, spacedim> DF[dim];
-      Tensor<1, spacedim> D2F[dim][dim];
+      Point<dim> p_unit = initial_p_unit;
 
-      Point<dim>     p_unit = initial_p_unit;
-      Point<dim>     f;
-      Tensor<2, dim> df;
+      const double       eps        = 1.e-12;
+      const unsigned int loop_limit = 10;
 
-      // Evaluate first and second derivatives
-      mdata.compute_shape_function_values(std::vector<Point<dim>>(1, p_unit));
+      unsigned int loop                   = 0;
+      double       f_weighted_norm_square = 1.;
 
-      for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
+      while (f_weighted_norm_square > eps * eps && loop++ < loop_limit)
         {
-          const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
-          const Tensor<2, dim> & hessian_k  = mdata.second_derivative(0, k);
-          const Point<spacedim> &point_k    = points[k];
-
+          const auto p_real =
+            internal::evaluate_tensor_product_value_and_gradient(
+              polynomials_1d,
+              points,
+              p_unit,
+              polynomials_1d.size() == 2,
+              renumber);
+          Tensor<1, spacedim>                    p_minus_F = p - p_real.first;
+          const DerivativeForm<1, spacedim, dim> DF        = p_real.second;
+
+          const auto hessian = internal::evaluate_tensor_product_hessian(
+            polynomials_1d, points, p_unit, renumber);
+          Point<dim>     f;
+          Tensor<2, dim> df;
           for (unsigned int j = 0; j < dim; ++j)
             {
-              DF[j] += grad_phi_k[j] * point_k;
+              f[j] = DF[j] * p_minus_F;
               for (unsigned int l = 0; l < dim; ++l)
-                D2F[j][l] += hessian_k[j][l] * point_k;
+                df[j][l] = -DF[j] * DF[l] + hessian[j][l] * p_minus_F;
             }
-        }
-
-      p_minus_F = p;
-      p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
-
-
-      for (unsigned int j = 0; j < dim; ++j)
-        f[j] = DF[j] * p_minus_F;
-
-      for (unsigned int j = 0; j < dim; ++j)
-        {
-          f[j] = DF[j] * p_minus_F;
-          for (unsigned int l = 0; l < dim; ++l)
-            df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
-        }
-
 
-      const double       eps        = 1.e-12 * cell->diameter();
-      const unsigned int loop_limit = 10;
-
-      unsigned int loop = 0;
-
-      while (f.norm() > eps && loop++ < loop_limit)
-        {
           // Solve  [df(x)]d=f(x)
           const Tensor<1, dim> d =
             invert(df) * static_cast<const Tensor<1, dim> &>(f);
+          f_weighted_norm_square = d.norm_square();
           p_unit -= d;
-
-          for (unsigned int j = 0; j < dim; ++j)
-            {
-              DF[j].clear();
-              for (unsigned int l = 0; l < dim; ++l)
-                D2F[j][l].clear();
-            }
-
-          mdata.compute_shape_function_values(
-            std::vector<Point<dim>>(1, p_unit));
-
-          for (unsigned int k = 0; k < mdata.n_shape_functions; ++k)
-            {
-              const Tensor<1, dim> & grad_phi_k = mdata.derivative(0, k);
-              const Tensor<2, dim> & hessian_k  = mdata.second_derivative(0, k);
-              const Point<spacedim> &point_k    = points[k];
-
-              for (unsigned int j = 0; j < dim; ++j)
-                {
-                  DF[j] += grad_phi_k[j] * point_k;
-                  for (unsigned int l = 0; l < dim; ++l)
-                    D2F[j][l] += hessian_k[j][l] * point_k;
-                }
-            }
-
-          // TODO: implement a line search here in much the same way as for
-          // the corresponding function above that does so for dim==spacedim
-          p_minus_F = p;
-          p_minus_F -= compute_mapped_location_of_point<dim, spacedim>(mdata);
-
-          for (unsigned int j = 0; j < dim; ++j)
-            {
-              f[j] = DF[j] * p_minus_F;
-              for (unsigned int l = 0; l < dim; ++l)
-                df[j][l] = -DF[j] * DF[l] + D2F[j][l] * p_minus_F;
-            }
         }
 
 
@@ -1131,6 +1073,15 @@ namespace internal
       constexpr unsigned int n_comp         = 1 + (spacedim - 1) / n_lanes;
       constexpr unsigned int n_hessians     = (dim * (dim + 1)) / 2;
 
+      // Since MappingQ::InternalData does not have separate arrays for the
+      // covariant and contravariant transformations, but uses the arrays in
+      // the `MappingRelatedData`, it can happen that vectors do not have the
+      // right size
+      if (update_flags & update_contravariant_transformation)
+        jacobians.resize(n_q_points);
+      if (update_flags & update_covariant_transformation)
+        inverse_jacobians.resize(n_q_points);
+
       EvaluationFlags::EvaluationFlags evaluation_flag =
         (update_flags & update_quadrature_points ? EvaluationFlags::values :
                                                    EvaluationFlags::nothing) |
@@ -1151,9 +1102,6 @@ namespace internal
       Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
                data.n_shape_functions > 0,
              ExcInternalError());
-      Assert(!(evaluation_flag & EvaluationFlags::gradients) ||
-               n_q_points == data.contravariant.size(),
-             ExcDimensionMismatch(n_q_points, data.contravariant.size()));
       Assert(!(evaluation_flag & EvaluationFlags::hessians) ||
                n_q_points == jacobian_grads.size(),
              ExcDimensionMismatch(n_q_points, jacobian_grads.size()));
@@ -1214,9 +1162,6 @@ namespace internal
 
       if (evaluation_flag & EvaluationFlags::gradients)
         {
-          std::fill(data.contravariant.begin(),
-                    data.contravariant.end(),
-                    DerivativeForm<1, dim, spacedim>());
           // We need to reinterpret the data after evaluate has been applied.
           for (unsigned int out_comp = 0; out_comp < n_comp; ++out_comp)
             for (unsigned int point = 0; point < n_q_points; ++point)
@@ -1229,43 +1174,28 @@ namespace internal
                     const unsigned int total_number = point * dim + j;
                     const unsigned int new_comp     = total_number / n_q_points;
                     const unsigned int new_point    = total_number % n_q_points;
-                    data.contravariant[new_point][out_comp * n_lanes + in_comp]
-                                      [new_comp] =
+                    jacobians[new_point][out_comp * n_lanes +
+                                         in_comp][new_comp] =
                       eval.begin_gradients()[(out_comp * n_q_points + point) *
                                                dim +
                                              j][in_comp];
                   }
         }
-      if (update_flags & update_covariant_transformation)
-        if (cell_similarity != CellSimilarity::translation)
-          for (unsigned int point = 0; point < n_q_points; ++point)
-            data.covariant[point] =
-              (data.contravariant[point]).covariant_form();
 
       if (update_flags & update_volume_elements)
         if (cell_similarity != CellSimilarity::translation)
           for (unsigned int point = 0; point < n_q_points; ++point)
-            data.volume_elements[point] =
-              data.contravariant[point].determinant();
+            data.volume_elements[point] = jacobians[point].determinant();
 
       // copy values from InternalData to vector given by reference
-      if (update_flags & update_jacobians)
+      if (update_flags & update_covariant_transformation)
         {
-          const unsigned int n_q_points = data.contravariant.size();
           AssertDimension(jacobians.size(), n_q_points);
-          if (cell_similarity != CellSimilarity::translation)
-            for (unsigned int point = 0; point < n_q_points; ++point)
-              jacobians[point] = data.contravariant[point];
-        }
-
-      // copy values from InternalData to vector given by reference
-      if (update_flags & update_inverse_jacobians)
-        {
-          const unsigned int n_q_points = data.contravariant.size();
           AssertDimension(inverse_jacobians.size(), n_q_points);
           if (cell_similarity != CellSimilarity::translation)
             for (unsigned int point = 0; point < n_q_points; ++point)
-              inverse_jacobians[point] = data.covariant[point].transpose();
+              inverse_jacobians[point] =
+                jacobians[point].covariant_form().transpose();
         }
 
       if (evaluation_flag & EvaluationFlags::hessians)
@@ -1309,6 +1239,40 @@ namespace internal
 
 
 
+    template <int dim, int spacedim>
+    inline DEAL_II_ALWAYS_INLINE void
+    store_vectorized_tensor(
+      const unsigned int n_points,
+      const unsigned int cur_index,
+      const DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+        &                                            derivative,
+      std::vector<DerivativeForm<1, dim, spacedim>> &result_array)
+    {
+      AssertDimension(result_array.size(), n_points);
+      constexpr unsigned int n_lanes = VectorizedArray<double>::size();
+      if (cur_index + n_lanes <= n_points)
+        {
+          std::array<unsigned int, n_lanes> indices;
+          for (unsigned int j = 0; j < n_lanes; ++j)
+            indices[j] = j * dim * spacedim;
+          const unsigned int even       = (dim * spacedim) / 4 * 4;
+          double *           result_ptr = &result_array[cur_index][0][0];
+          const VectorizedArray<double> *derivative_ptr = &derivative[0][0];
+          vectorized_transpose_and_store(
+            false, even, derivative_ptr, indices.data(), result_ptr);
+          for (unsigned int d = even; d < dim * spacedim; ++d)
+            for (unsigned int j = 0; j < n_lanes; ++j)
+              result_ptr[j * dim * spacedim + d] = derivative_ptr[d][j];
+        }
+      else
+        for (unsigned int j = 0; j < n_lanes && cur_index + j < n_points; ++j)
+          for (unsigned int d = 0; d < spacedim; ++d)
+            for (unsigned int e = 0; e < dim; ++e)
+              result_array[cur_index + j][d][e] = derivative[d][e][j];
+    }
+
+
+
     template <int dim, int spacedim>
     inline void
     maybe_update_q_points_Jacobians_generic(
@@ -1316,7 +1280,6 @@ namespace internal
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
       const ArrayView<const Point<dim>> &                           unit_points,
       const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
-      const unsigned int                                  polynomial_degree,
       const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       std::vector<Point<spacedim>> &   quadrature_points,
       std::vector<DerivativeForm<1, dim, spacedim>> &jacobians,
@@ -1329,6 +1292,15 @@ namespace internal
       const unsigned int n_points = unit_points.size();
       const unsigned int n_lanes  = VectorizedArray<double>::size();
 
+      // Since MappingQ::InternalData does not have separate arrays for the
+      // covariant and contravariant transformations, but uses the arrays in
+      // the `MappingRelatedData`, it can happen that vectors do not have the
+      // right size
+      if (update_flags & update_contravariant_transformation)
+        jacobians.resize(n_points);
+      if (update_flags & update_covariant_transformation)
+        inverse_jacobians.resize(n_points);
+
       // Use the more heavy VectorizedArray code path if there is more than
       // one point left to compute
       for (unsigned int i = 0; i < n_points; i += n_lanes)
@@ -1348,8 +1320,13 @@ namespace internal
                 polynomials_1d,
                 support_points,
                 p_vec,
-                polynomial_degree == 1,
+                polynomials_1d.size() == 2,
                 renumber_lexicographic_to_hierarchic);
+            DerivativeForm<1, dim, spacedim, VectorizedArray<double>>
+              derivative;
+            for (unsigned int d = 0; d < spacedim; ++d)
+              for (unsigned int e = 0; e < dim; ++e)
+                derivative[d][e] = result.second[e][d];
 
             if (update_flags & update_quadrature_points)
               for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
@@ -1360,28 +1337,23 @@ namespace internal
               continue;
 
             if (update_flags & update_contravariant_transformation)
-              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-                for (unsigned int d = 0; d < spacedim; ++d)
-                  for (unsigned int e = 0; e < dim; ++e)
-                    data.contravariant[i + j][d][e] = result.second[e][d][j];
+              store_vectorized_tensor(n_points, i, derivative, jacobians);
 
             if (update_flags & update_volume_elements)
-              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-                data.volume_elements[i + j] =
-                  data.contravariant[i + j].determinant();
-
-            if (update_flags & update_jacobians)
-              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-                jacobians[i + j] = data.contravariant[i + j];
+              {
+                const auto determinant = derivative.determinant();
+                for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
+                  data.volume_elements[i + j] = determinant[j];
+              }
 
             if (update_flags & update_covariant_transformation)
-              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-                data.covariant[i + j] =
-                  data.contravariant[i + j].covariant_form();
-
-            if (update_flags & update_inverse_jacobians)
-              for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-                inverse_jacobians[i + j] = data.covariant[i + j].transpose();
+              {
+                const auto covariant = derivative.covariant_form();
+                store_vectorized_tensor(n_points,
+                                        i,
+                                        covariant.transpose(),
+                                        inverse_jacobians);
+              }
           }
         else
           {
@@ -1390,7 +1362,7 @@ namespace internal
                 polynomials_1d,
                 support_points,
                 unit_points[i],
-                polynomial_degree == 1,
+                polynomials_1d.size() == 2,
                 renumber_lexicographic_to_hierarchic);
 
             if (update_flags & update_quadrature_points)
@@ -1399,23 +1371,16 @@ namespace internal
             if (cell_similarity == CellSimilarity::translation)
               continue;
 
-            if (update_flags & update_contravariant_transformation)
-              {
-                DerivativeForm<1, spacedim, dim> jac_transposed = result.second;
-                data.contravariant[i] = jac_transposed.transpose();
-              }
-
-            if (update_flags & update_volume_elements)
-              data.volume_elements[i] = data.contravariant[i].determinant();
+            if (dim == spacedim && update_flags & update_volume_elements)
+              data.volume_elements[i] =
+                DerivativeForm<1, spacedim, dim>(result.second).determinant();
 
-            if (update_flags & update_jacobians)
-              jacobians[i] = data.contravariant[i];
+            if (update_flags & update_contravariant_transformation)
+              jacobians[i] =
+                DerivativeForm<1, spacedim, dim>(result.second).transpose();
 
             if (update_flags & update_covariant_transformation)
-              data.covariant[i] = data.contravariant[i].covariant_form();
-
-            if (update_flags & update_inverse_jacobians)
-              inverse_jacobians[i] = data.covariant[i].transpose();
+              inverse_jacobians[i] = jacobians[i].covariant_form().transpose();
           }
     }
 
@@ -1430,38 +1395,33 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_grads(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       std::vector<DerivativeForm<2, dim, spacedim>> &jacobian_grads)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_grads)
+      if (data.update_each & update_jacobian_grads)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points = jacobian_grads.size();
 
           if (cell_similarity != CellSimilarity::translation)
             for (unsigned int point = 0; point < n_q_points; ++point)
               {
-                const Tensor<2, dim> *second =
-                  &data.second_derivative(point + data_set, 0);
-                double result[spacedim][dim][dim];
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  for (unsigned int j = 0; j < dim; ++j)
-                    for (unsigned int l = 0; l < dim; ++l)
-                      result[i][j][l] =
-                        (second[0][j][l] * data.mapping_support_points[0][i]);
-                for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        result[i][j][l] +=
-                          (second[k][j][l] * data.mapping_support_points[k][i]);
+                const SymmetricTensor<2, dim, Tensor<1, spacedim>> second =
+                  internal::evaluate_tensor_product_hessian(
+                    polynomials_1d,
+                    support_points,
+                    unit_points[point],
+                    renumber_lexicographic_to_hierarchic);
 
                 for (unsigned int i = 0; i < spacedim; ++i)
                   for (unsigned int j = 0; j < dim; ++j)
                     for (unsigned int l = 0; l < dim; ++l)
-                      jacobian_grads[point][i][j][l] = result[i][j][l];
+                      jacobian_grads[point][i][j][l] = second[j][l][i];
               }
         }
     }
@@ -1477,14 +1437,17 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_pushed_forward_grads(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
       std::vector<Tensor<3, spacedim>> &jacobian_pushed_forward_grads)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_grads)
+      if (data.update_each & update_jacobian_pushed_forward_grads)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points = jacobian_pushed_forward_grads.size();
 
           if (cell_similarity != CellSimilarity::translation)
@@ -1492,33 +1455,25 @@ namespace internal
               double tmp[spacedim][spacedim][spacedim];
               for (unsigned int point = 0; point < n_q_points; ++point)
                 {
-                  const Tensor<2, dim> *second =
-                    &data.second_derivative(point + data_set, 0);
-                  double result[spacedim][dim][dim];
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        result[i][j][l] =
-                          (second[0][j][l] * data.mapping_support_points[0][i]);
-                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                    for (unsigned int i = 0; i < spacedim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          result[i][j][l] +=
-                            (second[k][j][l] *
-                             data.mapping_support_points[k][i]);
+                  const SymmetricTensor<2, dim, Tensor<1, spacedim>> second =
+                    internal::evaluate_tensor_product_hessian(
+                      polynomials_1d,
+                      support_points,
+                      unit_points[point],
+                      renumber_lexicographic_to_hierarchic);
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[point].transpose();
 
                   // first push forward the j-components
                   for (unsigned int i = 0; i < spacedim; ++i)
                     for (unsigned int j = 0; j < spacedim; ++j)
                       for (unsigned int l = 0; l < dim; ++l)
                         {
-                          tmp[i][j][l] =
-                            result[i][0][l] * data.covariant[point][j][0];
+                          tmp[i][j][l] = second[0][l][i] * covariant[j][0];
                           for (unsigned int jr = 1; jr < dim; ++jr)
                             {
                               tmp[i][j][l] +=
-                                result[i][jr][l] * data.covariant[point][j][jr];
+                                second[jr][l][i] * covariant[j][jr];
                             }
                         }
 
@@ -1528,11 +1483,11 @@ namespace internal
                       for (unsigned int l = 0; l < spacedim; ++l)
                         {
                           jacobian_pushed_forward_grads[point][i][j][l] =
-                            tmp[i][j][0] * data.covariant[point][l][0];
+                            tmp[i][j][0] * covariant[l][0];
                           for (unsigned int lr = 1; lr < dim; ++lr)
                             {
                               jacobian_pushed_forward_grads[point][i][j][l] +=
-                                tmp[i][j][lr] * data.covariant[point][l][lr];
+                                tmp[i][j][lr] * covariant[l][lr];
                             }
                         }
                 }
@@ -1542,6 +1497,52 @@ namespace internal
 
 
 
+    template <int dim, int spacedim, int length_tensor>
+    inline DerivativeForm<3, dim, spacedim>
+    expand_3rd_derivatives(
+      const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
+    {
+      Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
+      DerivativeForm<3, dim, spacedim> result;
+      for (unsigned int i = 0; i < spacedim; ++i)
+        {
+          if (dim == 1)
+            result[i][0][0][0] = compressed[0][i];
+          else if (dim >= 2)
+            {
+              for (unsigned int d = 0; d < 2; ++d)
+                for (unsigned int e = 0; e < 2; ++e)
+                  for (unsigned int f = 0; f < 2; ++f)
+                    result[i][d][e][f] = compressed[d + e + f][i];
+              if (dim == 3)
+                {
+                  AssertDimension(length_tensor, 10);
+
+                  // the derivatives are ordered in rows by increasing z
+                  // derivative, and in each row we have x^{(n-j)} y^{(j)} as
+                  // j walks through the columns
+                  for (unsigned int d = 0; d < 2; ++d)
+                    for (unsigned int e = 0; e < 2; ++e)
+                      {
+                        result[i][d][e][2] = compressed[4 + d + e][i];
+                        result[i][d][2][e] = compressed[4 + d + e][i];
+                        result[i][2][d][e] = compressed[4 + d + e][i];
+                      }
+                  for (unsigned int d = 0; d < 2; ++d)
+                    {
+                      result[i][d][2][2] = compressed[7 + d][i];
+                      result[i][2][d][2] = compressed[7 + d][i];
+                      result[i][2][2][d] = compressed[7 + d][i];
+                    }
+                  result[i][2][2][2] = compressed[9][i];
+                }
+            }
+        }
+      return result;
+    }
+
+
+
     /**
      * Update the third derivatives of the transformation from unit to real
      * cell, the Jacobian hessians.
@@ -1551,45 +1552,29 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_2nd_derivatives(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       std::vector<DerivativeForm<3, dim, spacedim>> &jacobian_2nd_derivatives)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_2nd_derivatives)
+      if (data.update_each & update_jacobian_2nd_derivatives)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points = jacobian_2nd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
             {
               for (unsigned int point = 0; point < n_q_points; ++point)
                 {
-                  const Tensor<3, dim> *third =
-                    &data.third_derivative(point + data_set, 0);
-                  double result[spacedim][dim][dim][dim];
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          result[i][j][l][m] =
-                            (third[0][j][l][m] *
-                             data.mapping_support_points[0][i]);
-                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                    for (unsigned int i = 0; i < spacedim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            result[i][j][l][m] +=
-                              (third[k][j][l][m] *
-                               data.mapping_support_points[k][i]);
-
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          jacobian_2nd_derivatives[point][i][j][l][m] =
-                            result[i][j][l][m];
+                  jacobian_2nd_derivatives[point] = expand_3rd_derivatives<dim>(
+                    internal::evaluate_tensor_product_higher_derivatives<3>(
+                      polynomials_1d,
+                      support_points,
+                      unit_points[point],
+                      renumber_lexicographic_to_hierarchic));
                 }
             }
         }
@@ -1607,40 +1592,35 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_pushed_forward_2nd_derivatives(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
       std::vector<Tensor<4, spacedim>> &jacobian_pushed_forward_2nd_derivatives)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_2nd_derivatives)
+      if (data.update_each & update_jacobian_pushed_forward_2nd_derivatives)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points =
             jacobian_pushed_forward_2nd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
             {
-              double tmp[spacedim][spacedim][spacedim][spacedim];
+              dealii::ndarray<double, spacedim, spacedim, dim, dim>      tmp;
+              dealii::ndarray<double, spacedim, spacedim, spacedim, dim> tmp2;
               for (unsigned int point = 0; point < n_q_points; ++point)
                 {
-                  const Tensor<3, dim> *third =
-                    &data.third_derivative(point + data_set, 0);
-                  double result[spacedim][dim][dim][dim];
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          result[i][j][l][m] =
-                            (third[0][j][l][m] *
-                             data.mapping_support_points[0][i]);
-                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                    for (unsigned int i = 0; i < spacedim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            result[i][j][l][m] +=
-                              (third[k][j][l][m] *
-                               data.mapping_support_points[k][i]);
+                  const DerivativeForm<3, dim, spacedim> third =
+                    expand_3rd_derivatives<dim>(
+                      internal::evaluate_tensor_product_higher_derivatives<3>(
+                        polynomials_1d,
+                        support_points,
+                        unit_points[point],
+                        renumber_lexicographic_to_hierarchic));
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[point].transpose();
 
                   // push forward the j-coordinate
                   for (unsigned int i = 0; i < spacedim; ++i)
@@ -1648,15 +1628,11 @@ namespace internal
                       for (unsigned int l = 0; l < dim; ++l)
                         for (unsigned int m = 0; m < dim; ++m)
                           {
-                            jacobian_pushed_forward_2nd_derivatives
-                              [point][i][j][l][m] = result[i][0][l][m] *
-                                                    data.covariant[point][j][0];
+                            tmp[i][j][l][m] =
+                              third[i][0][l][m] * covariant[j][0];
                             for (unsigned int jr = 1; jr < dim; ++jr)
-                              jacobian_pushed_forward_2nd_derivatives[point][i]
-                                                                     [j][l]
-                                                                     [m] +=
-                                result[i][jr][l][m] *
-                                data.covariant[point][j][jr];
+                              tmp[i][j][l][m] +=
+                                third[i][jr][l][m] * covariant[j][jr];
                           }
 
                   // push forward the l-coordinate
@@ -1665,16 +1641,11 @@ namespace internal
                       for (unsigned int l = 0; l < spacedim; ++l)
                         for (unsigned int m = 0; m < dim; ++m)
                           {
-                            tmp[i][j][l][m] =
-                              jacobian_pushed_forward_2nd_derivatives[point][i]
-                                                                     [j][0][m] *
-                              data.covariant[point][l][0];
+                            tmp2[i][j][l][m] =
+                              tmp[i][j][0][m] * covariant[l][0];
                             for (unsigned int lr = 1; lr < dim; ++lr)
-                              tmp[i][j][l][m] +=
-                                jacobian_pushed_forward_2nd_derivatives[point]
-                                                                       [i][j]
-                                                                       [lr][m] *
-                                data.covariant[point][l][lr];
+                              tmp2[i][j][l][m] +=
+                                tmp[i][j][lr][m] * covariant[l][lr];
                           }
 
                   // push forward the m-coordinate
@@ -1685,12 +1656,12 @@ namespace internal
                           {
                             jacobian_pushed_forward_2nd_derivatives
                               [point][i][j][l][m] =
-                                tmp[i][j][l][0] * data.covariant[point][m][0];
+                                tmp2[i][j][l][0] * covariant[m][0];
                             for (unsigned int mr = 1; mr < dim; ++mr)
                               jacobian_pushed_forward_2nd_derivatives[point][i]
                                                                      [j][l]
                                                                      [m] +=
-                                tmp[i][j][l][mr] * data.covariant[point][m][mr];
+                                tmp2[i][j][l][mr] * covariant[m][mr];
                           }
                 }
             }
@@ -1699,6 +1670,66 @@ namespace internal
 
 
 
+    template <int dim, int spacedim, int length_tensor>
+    inline DerivativeForm<4, dim, spacedim>
+    expand_4th_derivatives(
+      const Tensor<1, length_tensor, Tensor<1, spacedim>> &compressed)
+    {
+      Assert(dim >= 1 && dim <= 3, ExcNotImplemented());
+      DerivativeForm<4, dim, spacedim> result;
+      for (unsigned int i = 0; i < spacedim; ++i)
+        {
+          if (dim == 1)
+            result[i][0][0][0][0] = compressed[0][i];
+          else if (dim >= 2)
+            {
+              for (unsigned int d = 0; d < 2; ++d)
+                for (unsigned int e = 0; e < 2; ++e)
+                  for (unsigned int f = 0; f < 2; ++f)
+                    for (unsigned int g = 0; g < 2; ++g)
+                      result[i][d][e][f][g] = compressed[d + e + f + g][i];
+              if (dim == 3)
+                {
+                  AssertDimension(length_tensor, 15);
+
+                  // the derivatives are ordered in rows by increasing z
+                  // derivative, and in each row we have x^{(n-j)} y^{(j)} as
+                  // j walks through the columns
+                  for (unsigned int d = 0; d < 2; ++d)
+                    for (unsigned int e = 0; e < 2; ++e)
+                      for (unsigned int f = 0; f < 2; ++f)
+                        {
+                          result[i][d][e][f][2] = compressed[5 + d + e + f][i];
+                          result[i][d][e][2][f] = compressed[5 + d + e + f][i];
+                          result[i][d][2][e][f] = compressed[5 + d + e + f][i];
+                          result[i][2][d][e][f] = compressed[5 + d + e + f][i];
+                        }
+                  for (unsigned int d = 0; d < 2; ++d)
+                    for (unsigned int e = 0; e < 2; ++e)
+                      {
+                        result[i][d][e][2][2] = compressed[9 + d + e][i];
+                        result[i][d][2][e][2] = compressed[9 + d + e][i];
+                        result[i][d][2][2][e] = compressed[9 + d + e][i];
+                        result[i][2][d][e][2] = compressed[9 + d + e][i];
+                        result[i][2][d][2][e] = compressed[9 + d + e][i];
+                        result[i][2][2][d][e] = compressed[9 + d + e][i];
+                      }
+                  for (unsigned int d = 0; d < 2; ++d)
+                    {
+                      result[i][d][2][2][2] = compressed[12 + d][i];
+                      result[i][2][d][2][2] = compressed[12 + d][i];
+                      result[i][2][2][d][2] = compressed[12 + d][i];
+                      result[i][2][2][2][d] = compressed[12 + d][i];
+                    }
+                  result[i][2][2][2][2] = compressed[14][i];
+                }
+            }
+        }
+      return result;
+    }
+
+
+
     /**
      * Update the fourth derivatives of the transformation from unit to real
      * cell, the Jacobian hessian gradients.
@@ -1708,48 +1739,29 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_3rd_derivatives(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       std::vector<DerivativeForm<4, dim, spacedim>> &jacobian_3rd_derivatives)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_3rd_derivatives)
+      if (data.update_each & update_jacobian_3rd_derivatives)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points = jacobian_3rd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
             {
               for (unsigned int point = 0; point < n_q_points; ++point)
                 {
-                  const Tensor<4, dim> *fourth =
-                    &data.fourth_derivative(point + data_set, 0);
-                  double result[spacedim][dim][dim][dim][dim];
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          for (unsigned int n = 0; n < dim; ++n)
-                            result[i][j][l][m][n] =
-                              (fourth[0][j][l][m][n] *
-                               data.mapping_support_points[0][i]);
-                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                    for (unsigned int i = 0; i < spacedim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[i][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.mapping_support_points[k][i]);
-
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          for (unsigned int n = 0; n < dim; ++n)
-                            jacobian_3rd_derivatives[point][i][j][l][m][n] =
-                              result[i][j][l][m][n];
+                  jacobian_3rd_derivatives[point] = expand_4th_derivatives<dim>(
+                    internal::evaluate_tensor_product_higher_derivatives<4>(
+                      polynomials_1d,
+                      support_points,
+                      unit_points[point],
+                      renumber_lexicographic_to_hierarchic));
                 }
             }
         }
@@ -1767,42 +1779,39 @@ namespace internal
     template <int dim, int spacedim>
     inline void
     maybe_update_jacobian_pushed_forward_3rd_derivatives(
-      const CellSimilarity::Similarity                  cell_similarity,
-      const typename QProjector<dim>::DataSetDescriptor data_set,
+      const CellSimilarity::Similarity cell_similarity,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
+      const ArrayView<const Point<dim>> &                           unit_points,
+      const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
+      const std::vector<unsigned int> & renumber_lexicographic_to_hierarchic,
       std::vector<Tensor<5, spacedim>> &jacobian_pushed_forward_3rd_derivatives)
     {
-      const UpdateFlags update_flags = data.update_each;
-      if (update_flags & update_jacobian_pushed_forward_3rd_derivatives)
+      if (data.update_each & update_jacobian_pushed_forward_3rd_derivatives)
         {
+          const std::vector<Point<spacedim>> &support_points =
+            data.mapping_support_points;
           const unsigned int n_q_points =
             jacobian_pushed_forward_3rd_derivatives.size();
 
           if (cell_similarity != CellSimilarity::translation)
             {
-              double tmp[spacedim][spacedim][spacedim][spacedim][spacedim];
+              dealii::
+                ndarray<double, spacedim, spacedim, spacedim, spacedim, dim>
+                  tmp;
+              dealii::ndarray<double, spacedim, spacedim, spacedim, dim, dim>
+                tmp2;
               for (unsigned int point = 0; point < n_q_points; ++point)
                 {
-                  const Tensor<4, dim> *fourth =
-                    &data.fourth_derivative(point + data_set, 0);
-                  double result[spacedim][dim][dim][dim][dim];
-                  for (unsigned int i = 0; i < spacedim; ++i)
-                    for (unsigned int j = 0; j < dim; ++j)
-                      for (unsigned int l = 0; l < dim; ++l)
-                        for (unsigned int m = 0; m < dim; ++m)
-                          for (unsigned int n = 0; n < dim; ++n)
-                            result[i][j][l][m][n] =
-                              (fourth[0][j][l][m][n] *
-                               data.mapping_support_points[0][i]);
-                  for (unsigned int k = 1; k < data.n_shape_functions; ++k)
-                    for (unsigned int i = 0; i < spacedim; ++i)
-                      for (unsigned int j = 0; j < dim; ++j)
-                        for (unsigned int l = 0; l < dim; ++l)
-                          for (unsigned int m = 0; m < dim; ++m)
-                            for (unsigned int n = 0; n < dim; ++n)
-                              result[i][j][l][m][n] +=
-                                (fourth[k][j][l][m][n] *
-                                 data.mapping_support_points[k][i]);
+                  const DerivativeForm<4, dim, spacedim> fourth =
+                    expand_4th_derivatives<dim>(
+                      internal::evaluate_tensor_product_higher_derivatives<4>(
+                        polynomials_1d,
+                        support_points,
+                        unit_points[point],
+                        renumber_lexicographic_to_hierarchic));
+
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[point].transpose();
 
                   // push-forward the j-coordinate
                   for (unsigned int i = 0; i < spacedim; ++i)
@@ -1811,12 +1820,11 @@ namespace internal
                         for (unsigned int m = 0; m < dim; ++m)
                           for (unsigned int n = 0; n < dim; ++n)
                             {
-                              tmp[i][j][l][m][n] = result[i][0][l][m][n] *
-                                                   data.covariant[point][j][0];
+                              tmp[i][j][l][m][n] =
+                                fourth[i][0][l][m][n] * covariant[j][0];
                               for (unsigned int jr = 1; jr < dim; ++jr)
                                 tmp[i][j][l][m][n] +=
-                                  result[i][jr][l][m][n] *
-                                  data.covariant[point][j][jr];
+                                  fourth[i][jr][l][m][n] * covariant[j][jr];
                             }
 
                   // push-forward the l-coordinate
@@ -1826,16 +1834,11 @@ namespace internal
                         for (unsigned int m = 0; m < dim; ++m)
                           for (unsigned int n = 0; n < dim; ++n)
                             {
-                              jacobian_pushed_forward_3rd_derivatives
-                                [point][i][j][l][m][n] =
-                                  tmp[i][j][0][m][n] *
-                                  data.covariant[point][l][0];
+                              tmp2[i][j][l][m][n] =
+                                tmp[i][j][0][m][n] * covariant[l][0];
                               for (unsigned int lr = 1; lr < dim; ++lr)
-                                jacobian_pushed_forward_3rd_derivatives[point]
-                                                                       [i][j][l]
-                                                                       [m][n] +=
-                                  tmp[i][j][lr][m][n] *
-                                  data.covariant[point][l][lr];
+                                tmp2[i][j][l][m][n] +=
+                                  tmp[i][j][lr][m][n] * covariant[l][lr];
                             }
 
                   // push-forward the m-coordinate
@@ -1846,15 +1849,10 @@ namespace internal
                           for (unsigned int n = 0; n < dim; ++n)
                             {
                               tmp[i][j][l][m][n] =
-                                jacobian_pushed_forward_3rd_derivatives[point]
-                                                                       [i][j][l]
-                                                                       [0][n] *
-                                data.covariant[point][m][0];
+                                tmp2[i][j][l][0][n] * covariant[m][0];
                               for (unsigned int mr = 1; mr < dim; ++mr)
                                 tmp[i][j][l][m][n] +=
-                                  jacobian_pushed_forward_3rd_derivatives
-                                    [point][i][j][l][mr][n] *
-                                  data.covariant[point][m][mr];
+                                  tmp2[i][j][l][mr][n] * covariant[m][mr];
                             }
 
                   // push-forward the n-coordinate
@@ -1866,14 +1864,12 @@ namespace internal
                             {
                               jacobian_pushed_forward_3rd_derivatives
                                 [point][i][j][l][m][n] =
-                                  tmp[i][j][l][m][0] *
-                                  data.covariant[point][n][0];
+                                  tmp[i][j][l][m][0] * covariant[n][0];
                               for (unsigned int nr = 1; nr < dim; ++nr)
                                 jacobian_pushed_forward_3rd_derivatives[point]
                                                                        [i][j][l]
                                                                        [m][n] +=
-                                  tmp[i][j][l][m][nr] *
-                                  data.covariant[point][n][nr];
+                                  tmp[i][j][l][m][nr] * covariant[n][nr];
                             }
                 }
             }
@@ -1907,8 +1903,7 @@ namespace internal
       const UpdateFlags update_flags = data.update_each;
 
       if (update_flags &
-          (update_boundary_forms | update_normal_vectors | update_jacobians |
-           update_JxW_values | update_inverse_jacobians))
+          (update_boundary_forms | update_normal_vectors | update_JxW_values))
         {
           if (update_flags & update_boundary_forms)
             AssertDimension(output_data.boundary_forms.size(), n_q_points);
@@ -1927,24 +1922,18 @@ namespace internal
           // 0.
           for (unsigned int d = 0; d != dim - 1; ++d)
             {
-              Assert(face_no + GeometryInfo<dim>::faces_per_cell * d <
-                       data.unit_tangentials.size(),
+              const unsigned int vector_index =
+                face_no + GeometryInfo<dim>::faces_per_cell * d;
+              Assert(vector_index < data.unit_tangentials.size(),
                      ExcInternalError());
-              Assert(
-                data.aux[d].size() <=
-                  data
-                    .unit_tangentials[face_no +
-                                      GeometryInfo<dim>::faces_per_cell * d]
-                    .size(),
-                ExcInternalError());
-
-              mapping.transform(
-                make_array_view(
-                  data.unit_tangentials[face_no +
-                                        GeometryInfo<dim>::faces_per_cell * d]),
-                mapping_contravariant,
-                data,
-                make_array_view(data.aux[d].begin(), data.aux[d].end()));
+              Assert(data.aux[d].size() <=
+                       data.unit_tangentials[vector_index].size(),
+                     ExcInternalError());
+              mapping.transform(make_array_view(
+                                  data.unit_tangentials[vector_index]),
+                                mapping_contravariant,
+                                data,
+                                make_array_view(data.aux[d]));
             }
 
           if (update_flags & update_boundary_forms)
@@ -1985,15 +1974,18 @@ namespace internal
                   //
                   // to compute the cell normal, use the same method used in
                   // fill_fe_values for cells above
-                  AssertDimension(data.contravariant.size(), n_q_points);
+                  AssertDimension(data.output_data->jacobians.size(),
+                                  n_q_points);
 
                   for (unsigned int point = 0; point < n_q_points; ++point)
                     {
+                      const DerivativeForm<1, dim, spacedim> contravariant =
+                        data.output_data->jacobians[point];
                       if (dim == 1)
                         {
                           // J is a tangent vector
                           output_data.boundary_forms[point] =
-                            data.contravariant[point].transpose()[0];
+                            contravariant.transpose()[0];
                           output_data.boundary_forms[point] /=
                             (face_no == 0 ? -1. : +1.) *
                             output_data.boundary_forms[point].norm();
@@ -2002,7 +1994,7 @@ namespace internal
                       if (dim == 2)
                         {
                           const DerivativeForm<1, spacedim, dim> DX_t =
-                            data.contravariant[point].transpose();
+                            contravariant.transpose();
 
                           Tensor<1, spacedim> cell_normal =
                             cross_product_3d(DX_t[0], DX_t[1]);
@@ -2058,11 +2050,12 @@ namespace internal
       const Quadrature<dim - 1> &                                   quadrature,
       const typename dealii::MappingQ<dim, spacedim>::InternalData &data,
       const std::vector<Polynomials::Polynomial<double>> &polynomials_1d,
-      const unsigned int                                  polynomial_degree,
       const std::vector<unsigned int> &renumber_lexicographic_to_hierarchic,
       internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
         &output_data)
     {
+      const ArrayView<const Point<dim>> quadrature_points(
+        &data.quadrature_points[data_set], quadrature.size());
       if (dim > 1 && data.tensor_product_quadrature)
         {
           maybe_update_q_points_Jacobians_and_grads_tensor<dim, spacedim>(
@@ -2079,42 +2072,54 @@ namespace internal
             maybe_update_q_points_Jacobians_generic(
               CellSimilarity::none,
               data,
-              make_array_view(&data.quadrature_points[data_set],
-                              &data.quadrature_points[data_set] +
-                                quadrature.size()),
+              quadrature_points,
               polynomials_1d,
-              polynomial_degree,
               renumber_lexicographic_to_hierarchic,
               output_data.quadrature_points,
               output_data.jacobians,
               output_data.inverse_jacobians);
           maybe_update_jacobian_grads<dim, spacedim>(
-            CellSimilarity::none, data_set, data, output_data.jacobian_grads);
+            CellSimilarity::none,
+            data,
+            quadrature_points,
+            polynomials_1d,
+            renumber_lexicographic_to_hierarchic,
+            output_data.jacobian_grads);
         }
       maybe_update_jacobian_pushed_forward_grads<dim, spacedim>(
         CellSimilarity::none,
-        data_set,
         data,
+        quadrature_points,
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_pushed_forward_grads);
       maybe_update_jacobian_2nd_derivatives<dim, spacedim>(
         CellSimilarity::none,
-        data_set,
         data,
+        quadrature_points,
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_2nd_derivatives);
       maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
         CellSimilarity::none,
-        data_set,
         data,
+        quadrature_points,
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_pushed_forward_2nd_derivatives);
       maybe_update_jacobian_3rd_derivatives<dim, spacedim>(
         CellSimilarity::none,
-        data_set,
         data,
+        quadrature_points,
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_3rd_derivatives);
       maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
         CellSimilarity::none,
-        data_set,
         data,
+        quadrature_points,
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_pushed_forward_3rd_derivatives);
 
       maybe_compute_face_data(mapping,
@@ -2159,8 +2164,8 @@ namespace internal
                        "update_contravariant_transformation"));
 
               for (unsigned int i = 0; i < output.size(); ++i)
-                output[i] =
-                  apply_transformation(data.contravariant[i], input[i]);
+                output[i] = apply_transformation(data.output_data->jacobians[i],
+                                                 input[i]);
 
               return;
             }
@@ -2180,7 +2185,8 @@ namespace internal
               for (unsigned int i = 0; i < output.size(); ++i)
                 {
                   output[i] =
-                    apply_transformation(data.contravariant[i], input[i]);
+                    apply_transformation(data.output_data->jacobians[i],
+                                         input[i]);
                   output[i] /= data.volume_elements[i];
                 }
               return;
@@ -2195,7 +2201,8 @@ namespace internal
                        "update_covariant_transformation"));
 
               for (unsigned int i = 0; i < output.size(); ++i)
-                output[i] = apply_transformation(data.covariant[i], input[i]);
+                output[i] = apply_transformation(
+                  data.output_data->inverse_jacobians[i].transpose(), input[i]);
 
               return;
             }
@@ -2243,10 +2250,11 @@ namespace internal
               for (unsigned int i = 0; i < output.size(); ++i)
                 {
                   const DerivativeForm<1, spacedim, dim> A =
-                    apply_transformation(data.contravariant[i],
+                    apply_transformation(data.output_data->jacobians[i],
                                          transpose(input[i]));
-                  output[i] =
-                    apply_transformation(data.covariant[i], A.transpose());
+                  output[i] = apply_transformation(
+                    data.output_data->inverse_jacobians[i].transpose(),
+                    A.transpose());
                 }
 
               return;
@@ -2261,11 +2269,11 @@ namespace internal
 
               for (unsigned int i = 0; i < output.size(); ++i)
                 {
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[i].transpose();
                   const DerivativeForm<1, spacedim, dim> A =
-                    apply_transformation(data.covariant[i],
-                                         transpose(input[i]));
-                  output[i] =
-                    apply_transformation(data.covariant[i], A.transpose());
+                    apply_transformation(covariant, transpose(input[i]));
+                  output[i] = apply_transformation(covariant, A.transpose());
                 }
 
               return;
@@ -2286,10 +2294,13 @@ namespace internal
 
               for (unsigned int i = 0; i < output.size(); ++i)
                 {
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[i].transpose();
                   const DerivativeForm<1, spacedim, dim> A =
-                    apply_transformation(data.covariant[i], input[i]);
+                    apply_transformation(covariant, input[i]);
                   const Tensor<2, spacedim> T =
-                    apply_transformation(data.contravariant[i], A.transpose());
+                    apply_transformation(data.output_data->jacobians[i],
+                                         A.transpose());
 
                   output[i] = transpose(T);
                   output[i] /= data.volume_elements[i];
@@ -2338,37 +2349,42 @@ namespace internal
                        "update_contravariant_transformation"));
 
               for (unsigned int q = 0; q < output.size(); ++q)
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  {
-                    double tmp1[dim][dim];
-                    for (unsigned int J = 0; J < dim; ++J)
-                      for (unsigned int K = 0; K < dim; ++K)
-                        {
-                          tmp1[J][K] =
-                            data.contravariant[q][i][0] * input[q][0][J][K];
-                          for (unsigned int I = 1; I < dim; ++I)
-                            tmp1[J][K] +=
-                              data.contravariant[q][i][I] * input[q][I][J][K];
-                        }
-                    for (unsigned int j = 0; j < spacedim; ++j)
-                      {
-                        double tmp2[dim];
+                {
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[q].transpose();
+                  const DerivativeForm<1, dim, spacedim> contravariant =
+                    data.output_data->jacobians[q];
+
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    {
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
                         for (unsigned int K = 0; K < dim; ++K)
                           {
-                            tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
-                            for (unsigned int J = 1; J < dim; ++J)
-                              tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+                            tmp1[J][K] =
+                              contravariant[i][0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] +=
+                                contravariant[i][I] * input[q][I][J][K];
                           }
-                        for (unsigned int k = 0; k < spacedim; ++k)
-                          {
-                            output[q][i][j][k] =
-                              data.covariant[q][k][0] * tmp2[0];
-                            for (unsigned int K = 1; K < dim; ++K)
-                              output[q][i][j][k] +=
-                                data.covariant[q][k][K] * tmp2[K];
-                          }
-                      }
-                  }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = covariant[j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += covariant[j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] = covariant[k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] += covariant[k][K] * tmp2[K];
+                            }
+                        }
+                    }
+                }
               return;
             }
 
@@ -2379,37 +2395,38 @@ namespace internal
                        "update_covariant_transformation"));
 
               for (unsigned int q = 0; q < output.size(); ++q)
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  {
-                    double tmp1[dim][dim];
-                    for (unsigned int J = 0; J < dim; ++J)
-                      for (unsigned int K = 0; K < dim; ++K)
-                        {
-                          tmp1[J][K] =
-                            data.covariant[q][i][0] * input[q][0][J][K];
-                          for (unsigned int I = 1; I < dim; ++I)
-                            tmp1[J][K] +=
-                              data.covariant[q][i][I] * input[q][I][J][K];
-                        }
-                    for (unsigned int j = 0; j < spacedim; ++j)
-                      {
-                        double tmp2[dim];
+                {
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[q].transpose();
+
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    {
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
                         for (unsigned int K = 0; K < dim; ++K)
                           {
-                            tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
-                            for (unsigned int J = 1; J < dim; ++J)
-                              tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
-                          }
-                        for (unsigned int k = 0; k < spacedim; ++k)
-                          {
-                            output[q][i][j][k] =
-                              data.covariant[q][k][0] * tmp2[0];
-                            for (unsigned int K = 1; K < dim; ++K)
-                              output[q][i][j][k] +=
-                                data.covariant[q][k][K] * tmp2[K];
+                            tmp1[J][K] = covariant[i][0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] += covariant[i][I] * input[q][I][J][K];
                           }
-                      }
-                  }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = covariant[j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += covariant[j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] = covariant[k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] += covariant[k][K] * tmp2[K];
+                            }
+                        }
+                    }
+                }
 
               return;
             }
@@ -2427,39 +2444,43 @@ namespace internal
                        "update_volume_elements"));
 
               for (unsigned int q = 0; q < output.size(); ++q)
-                for (unsigned int i = 0; i < spacedim; ++i)
-                  {
-                    double factor[dim];
-                    for (unsigned int I = 0; I < dim; ++I)
-                      factor[I] =
-                        data.contravariant[q][i][I] / data.volume_elements[q];
-                    double tmp1[dim][dim];
-                    for (unsigned int J = 0; J < dim; ++J)
-                      for (unsigned int K = 0; K < dim; ++K)
-                        {
-                          tmp1[J][K] = factor[0] * input[q][0][J][K];
-                          for (unsigned int I = 1; I < dim; ++I)
-                            tmp1[J][K] += factor[I] * input[q][I][J][K];
-                        }
-                    for (unsigned int j = 0; j < spacedim; ++j)
-                      {
-                        double tmp2[dim];
+                {
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.output_data->inverse_jacobians[q].transpose();
+                  const DerivativeForm<1, dim, spacedim> contravariant =
+                    data.output_data->jacobians[q];
+                  for (unsigned int i = 0; i < spacedim; ++i)
+                    {
+                      double factor[dim];
+                      for (unsigned int I = 0; I < dim; ++I)
+                        factor[I] =
+                          contravariant[i][I] * (1. / data.volume_elements[q]);
+                      double tmp1[dim][dim];
+                      for (unsigned int J = 0; J < dim; ++J)
                         for (unsigned int K = 0; K < dim; ++K)
                           {
-                            tmp2[K] = data.covariant[q][j][0] * tmp1[0][K];
-                            for (unsigned int J = 1; J < dim; ++J)
-                              tmp2[K] += data.covariant[q][j][J] * tmp1[J][K];
+                            tmp1[J][K] = factor[0] * input[q][0][J][K];
+                            for (unsigned int I = 1; I < dim; ++I)
+                              tmp1[J][K] += factor[I] * input[q][I][J][K];
                           }
-                        for (unsigned int k = 0; k < spacedim; ++k)
-                          {
-                            output[q][i][j][k] =
-                              data.covariant[q][k][0] * tmp2[0];
-                            for (unsigned int K = 1; K < dim; ++K)
-                              output[q][i][j][k] +=
-                                data.covariant[q][k][K] * tmp2[K];
-                          }
-                      }
-                  }
+                      for (unsigned int j = 0; j < spacedim; ++j)
+                        {
+                          double tmp2[dim];
+                          for (unsigned int K = 0; K < dim; ++K)
+                            {
+                              tmp2[K] = covariant[j][0] * tmp1[0][K];
+                              for (unsigned int J = 1; J < dim; ++J)
+                                tmp2[K] += covariant[j][J] * tmp1[J][K];
+                            }
+                          for (unsigned int k = 0; k < spacedim; ++k)
+                            {
+                              output[q][i][j][k] = covariant[k][0] * tmp2[0];
+                              for (unsigned int K = 1; K < dim; ++K)
+                                output[q][i][j][k] += covariant[k][K] * tmp2[K];
+                            }
+                        }
+                    }
+                }
 
               return;
             }
@@ -2502,7 +2523,8 @@ namespace internal
                        "update_covariant_transformation"));
 
               for (unsigned int i = 0; i < output.size(); ++i)
-                output[i] = apply_transformation(data.covariant[i], input[i]);
+                output[i] = apply_transformation(
+                  data.output_data->inverse_jacobians[i].transpose(), input[i]);
 
               return;
             }
index 2961f7302611fecd64508fc78c0db04c7c722c17..b1e8eb1676db3b30361c6ac2af648f62e137c52c 100644 (file)
@@ -54,6 +54,7 @@ MappingQ<dim, spacedim>::InternalData::InternalData(
   , n_shape_functions(Utilities::fixed_power<dim>(polynomial_degree + 1))
   , line_support_points(QGaussLobatto<1>(polynomial_degree + 1))
   , tensor_product_quadrature(false)
+  , output_data(nullptr)
 {}
 
 
@@ -64,10 +65,7 @@ MappingQ<dim, spacedim>::InternalData::memory_consumption() const
 {
   return (
     Mapping<dim, spacedim>::InternalDataBase::memory_consumption() +
-    MemoryConsumption::memory_consumption(shape_values) +
-    MemoryConsumption::memory_consumption(shape_derivatives) +
-    MemoryConsumption::memory_consumption(covariant) +
-    MemoryConsumption::memory_consumption(contravariant) +
+    MemoryConsumption::memory_consumption(quadrature_points) +
     MemoryConsumption::memory_consumption(unit_tangentials) +
     MemoryConsumption::memory_consumption(aux) +
     MemoryConsumption::memory_consumption(mapping_support_points) +
@@ -92,12 +90,6 @@ MappingQ<dim, spacedim>::InternalData::initialize(
 
   const unsigned int n_q_points = quadrature.size();
 
-  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);
 
@@ -161,49 +153,6 @@ MappingQ<dim, spacedim>::InternalData::initialize(
             }
         }
     }
-
-  const bool needs_higher_order_terms =
-    this->update_each &
-    (update_jacobian_pushed_forward_grads | update_jacobian_2nd_derivatives |
-     update_jacobian_pushed_forward_2nd_derivatives |
-     update_jacobian_3rd_derivatives |
-     update_jacobian_pushed_forward_3rd_derivatives);
-
-  const bool needs_higher_order_terms_generic =
-    !tensor_product_quadrature &&
-    (needs_higher_order_terms || this->update_each & update_jacobian_grads);
-
-  // Only fill the big arrays on demand in case we cannot use the tensor
-  // product quadrature code path
-  if (dim == 1 || needs_higher_order_terms_generic || needs_higher_order_terms)
-    {
-      // compute shapes and derivatives for codim1 (for
-      // do_transform_real_to_unit_cell_internal_codim1)
-      if (dim == (spacedim - 1))
-        {
-          // see if we need the (transformation) shape function values
-          // and/or gradients and resize the necessary arrays
-          if (this->update_each & update_quadrature_points)
-            shape_values.resize(n_shape_functions * n_q_points);
-          if (this->update_each & update_jacobians)
-            shape_derivatives.resize(n_shape_functions * n_q_points);
-        }
-
-      if (this->update_each &
-          (update_jacobian_grads | update_jacobian_pushed_forward_grads))
-        shape_second_derivatives.resize(n_shape_functions * n_q_points);
-
-      if (this->update_each & (update_jacobian_2nd_derivatives |
-                               update_jacobian_pushed_forward_2nd_derivatives))
-        shape_third_derivatives.resize(n_shape_functions * n_q_points);
-
-      if (this->update_each & (update_jacobian_3rd_derivatives |
-                               update_jacobian_pushed_forward_3rd_derivatives))
-        shape_fourth_derivatives.resize(n_shape_functions * n_q_points);
-
-      // now also fill the various fields with their correct values
-      compute_shape_function_values(quadrature.get_points());
-    }
 }
 
 
@@ -217,10 +166,7 @@ MappingQ<dim, spacedim>::InternalData::initialize_face(
 {
   initialize(update_flags, quadrature, n_original_q_points);
 
-  const unsigned int n_q_points = quadrature.size();
-  quadrature_points.resize(n_q_points);
-  for (unsigned int q = 0; q < n_q_points; ++q)
-    quadrature_points[q] = quadrature.get_points()[q];
+  quadrature_points = quadrature.get_points();
 
   if (dim > 1 && tensor_product_quadrature)
     {
@@ -238,8 +184,7 @@ MappingQ<dim, spacedim>::InternalData::initialize_face(
   if (dim > 1)
     {
       if (this->update_each &
-          (update_boundary_forms | update_normal_vectors | update_jacobians |
-           update_JxW_values | update_inverse_jacobians))
+          (update_boundary_forms | update_normal_vectors | update_JxW_values))
         {
           aux.resize(dim - 1,
                      AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
@@ -269,98 +214,6 @@ MappingQ<dim, spacedim>::InternalData::initialize_face(
 
 
 
-template <int dim, int spacedim>
-void
-MappingQ<dim, spacedim>::InternalData::compute_shape_function_values(
-  const std::vector<Point<dim>> &unit_points)
-{
-  const unsigned int n_points = unit_points.size();
-
-  // Construct the tensor product polynomials used as shape functions for
-  // the Qp mapping of cells at the boundary.
-  const TensorProductPolynomials<dim> tensor_pols(
-    Polynomials::generate_complete_Lagrange_basis(
-      line_support_points.get_points()));
-  Assert(n_shape_functions == tensor_pols.n(), ExcInternalError());
-
-  // then also construct the mapping from lexicographic to the Qp shape
-  // function numbering
-  const std::vector<unsigned int> renumber =
-    FETools::hierarchic_to_lexicographic_numbering<dim>(polynomial_degree);
-
-  std::vector<double>         values;
-  std::vector<Tensor<1, dim>> grads;
-  if (shape_values.size() != 0)
-    {
-      Assert(shape_values.size() == n_shape_functions * n_points,
-             ExcInternalError());
-      values.resize(n_shape_functions);
-    }
-  if (shape_derivatives.size() != 0)
-    {
-      Assert(shape_derivatives.size() == n_shape_functions * n_points,
-             ExcInternalError());
-      grads.resize(n_shape_functions);
-    }
-
-  std::vector<Tensor<2, dim>> grad2;
-  if (shape_second_derivatives.size() != 0)
-    {
-      Assert(shape_second_derivatives.size() == n_shape_functions * n_points,
-             ExcInternalError());
-      grad2.resize(n_shape_functions);
-    }
-
-  std::vector<Tensor<3, dim>> grad3;
-  if (shape_third_derivatives.size() != 0)
-    {
-      Assert(shape_third_derivatives.size() == n_shape_functions * n_points,
-             ExcInternalError());
-      grad3.resize(n_shape_functions);
-    }
-
-  std::vector<Tensor<4, dim>> grad4;
-  if (shape_fourth_derivatives.size() != 0)
-    {
-      Assert(shape_fourth_derivatives.size() == n_shape_functions * n_points,
-             ExcInternalError());
-      grad4.resize(n_shape_functions);
-    }
-
-
-  if (shape_values.size() != 0 || shape_derivatives.size() != 0 ||
-      shape_second_derivatives.size() != 0 ||
-      shape_third_derivatives.size() != 0 ||
-      shape_fourth_derivatives.size() != 0)
-    for (unsigned int point = 0; point < n_points; ++point)
-      {
-        tensor_pols.evaluate(
-          unit_points[point], values, grads, grad2, grad3, grad4);
-
-        if (shape_values.size() != 0)
-          for (unsigned int i = 0; i < n_shape_functions; ++i)
-            shape(point, i) = values[renumber[i]];
-
-        if (shape_derivatives.size() != 0)
-          for (unsigned int i = 0; i < n_shape_functions; ++i)
-            derivative(point, i) = grads[renumber[i]];
-
-        if (shape_second_derivatives.size() != 0)
-          for (unsigned int i = 0; i < n_shape_functions; ++i)
-            second_derivative(point, i) = grad2[renumber[i]];
-
-        if (shape_third_derivatives.size() != 0)
-          for (unsigned int i = 0; i < n_shape_functions; ++i)
-            third_derivative(point, i) = grad3[renumber[i]];
-
-        if (shape_fourth_derivatives.size() != 0)
-          for (unsigned int i = 0; i < n_shape_functions; ++i)
-            fourth_derivative(point, i) = grad4[renumber[i]];
-      }
-}
-
-
-
 template <int dim, int spacedim>
 MappingQ<dim, spacedim>::MappingQ(const unsigned int p)
   : polynomial_degree(p)
@@ -392,30 +245,8 @@ MappingQ<dim, spacedim>::MappingQ(const unsigned int p)
 
 template <int dim, int spacedim>
 MappingQ<dim, spacedim>::MappingQ(const unsigned int p, const bool)
-  : polynomial_degree(p)
-  , line_support_points(
-      QGaussLobatto<1>(this->polynomial_degree + 1).get_points())
-  , polynomials_1d(
-      Polynomials::generate_complete_Lagrange_basis(line_support_points))
-  , renumber_lexicographic_to_hierarchic(
-      FETools::lexicographic_to_hierarchic_numbering<dim>(p))
-  , unit_cell_support_points(
-      internal::MappingQImplementation::unit_support_points<dim>(
-        line_support_points,
-        renumber_lexicographic_to_hierarchic))
-  , support_point_weights_perimeter_to_interior(
-      internal::MappingQImplementation::
-        compute_support_point_weights_perimeter_to_interior(
-          this->polynomial_degree,
-          dim))
-  , support_point_weights_cell(
-      internal::MappingQImplementation::compute_support_point_weights_cell<dim>(
-        this->polynomial_degree))
-{
-  Assert(p >= 1,
-         ExcMessage("It only makes sense to create polynomial mappings "
-                    "with a polynomial degree greater or equal to one."));
-}
+  : MappingQ<dim, spacedim>(p)
+{}
 
 
 
@@ -579,10 +410,12 @@ MappingQ<1, 2>::transform_real_to_unit_cell_internal(
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
   return internal::MappingQImplementation::
-    do_transform_real_to_unit_cell_internal_codim1<1>(cell,
-                                                      p,
-                                                      initial_p_unit,
-                                                      *mdata);
+    do_transform_real_to_unit_cell_internal_codim1<1>(
+      p,
+      initial_p_unit,
+      mdata->mapping_support_points,
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic);
 }
 
 
@@ -610,12 +443,16 @@ MappingQ<2, 3>::transform_real_to_unit_cell_internal(
   // dispatch to the various specializations for spacedim=dim,
   // spacedim=dim+1, etc
   return internal::MappingQImplementation::
-    do_transform_real_to_unit_cell_internal_codim1<2>(cell,
-                                                      p,
-                                                      initial_p_unit,
-                                                      *mdata);
+    do_transform_real_to_unit_cell_internal_codim1<2>(
+      p,
+      initial_p_unit,
+      mdata->mapping_support_points,
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic);
 }
 
+
+
 template <>
 Point<1>
 MappingQ<1, 3>::transform_real_to_unit_cell_internal(
@@ -857,9 +694,9 @@ MappingQ<dim, spacedim>::requires_update_flags(const UpdateFlags in) const
       if (out & (update_JxW_values | update_normal_vectors))
         out |= update_boundary_forms;
 
-      if (out & (update_covariant_transformation | update_JxW_values |
-                 update_jacobians | update_jacobian_grads |
-                 update_boundary_forms | update_normal_vectors))
+      if (out &
+          (update_covariant_transformation | update_jacobian_grads |
+           update_jacobians | update_boundary_forms | update_normal_vectors))
         out |= update_contravariant_transformation;
 
       if (out &
@@ -957,6 +794,7 @@ MappingQ<dim, spacedim>::fill_fe_values(
   Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
+  data.output_data         = &output_data;
 
   const unsigned int n_q_points = quadrature.size();
 
@@ -1000,7 +838,6 @@ MappingQ<dim, spacedim>::fill_fe_values(
         data,
         make_array_view(quadrature.get_points()),
         polynomials_1d,
-        polynomial_degree,
         renumber_lexicographic_to_hierarchic,
         output_data.quadrature_points,
         output_data.jacobians,
@@ -1009,44 +846,56 @@ MappingQ<dim, spacedim>::fill_fe_values(
       internal::MappingQImplementation::maybe_update_jacobian_grads<dim,
                                                                     spacedim>(
         computed_cell_similarity,
-        QProjector<dim>::DataSetDescriptor::cell(),
         data,
+        make_array_view(quadrature.get_points()),
+        polynomials_1d,
+        renumber_lexicographic_to_hierarchic,
         output_data.jacobian_grads);
     }
 
   internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads<
     dim,
     spacedim>(computed_cell_similarity,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_pushed_forward_grads);
 
   internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives<
     dim,
     spacedim>(computed_cell_similarity,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_2nd_derivatives);
 
   internal::MappingQImplementation::
     maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
       computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell(),
       data,
+      make_array_view(quadrature.get_points()),
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic,
       output_data.jacobian_pushed_forward_2nd_derivatives);
 
   internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives<
     dim,
     spacedim>(computed_cell_similarity,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_3rd_derivatives);
 
   internal::MappingQImplementation::
     maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
       computed_cell_similarity,
-      QProjector<dim>::DataSetDescriptor::cell(),
       data,
+      make_array_view(quadrature.get_points()),
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic,
       output_data.jacobian_pushed_forward_3rd_derivatives);
 
   const UpdateFlags          update_flags = data.update_each;
@@ -1070,7 +919,7 @@ MappingQ<dim, spacedim>::fill_fe_values(
           {
             if (dim == spacedim)
               {
-                const double det = data.contravariant[point].determinant();
+                const double det = data.volume_elements[point];
 
                 // check for distorted cells.
 
@@ -1093,7 +942,7 @@ MappingQ<dim, spacedim>::fill_fe_values(
                 Tensor<1, spacedim> DX_t[dim];
                 for (unsigned int i = 0; i < spacedim; ++i)
                   for (unsigned int j = 0; j < dim; ++j)
-                    DX_t[j][i] = data.contravariant[point][i][j];
+                    DX_t[j][i] = output_data.jacobians[point][i][j];
 
                 Tensor<2, dim> G; // First fundamental form
                 for (unsigned int i = 0; i < dim; ++i)
@@ -1163,11 +1012,12 @@ MappingQ<dim, spacedim>::fill_fe_face_values(
   Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
+  data.output_data         = &output_data;
 
   // if necessary, recompute the support points of the transformation of this
   // cell (note that we need to first check the triangulation pointer, since
-  // otherwise the second test might trigger an exception if the triangulations
-  // are not the same)
+  // otherwise the second test might trigger an exception if the
+  // triangulations are not the same)
   if ((data.mapping_support_points.size() == 0) ||
       (&cell->get_triangulation() !=
        &data.cell_of_current_support_points->get_triangulation()) ||
@@ -1192,7 +1042,6 @@ MappingQ<dim, spacedim>::fill_fe_face_values(
     quadrature[0],
     data,
     polynomials_1d,
-    polynomial_degree,
     renumber_lexicographic_to_hierarchic,
     output_data);
 }
@@ -1214,11 +1063,12 @@ MappingQ<dim, spacedim>::fill_fe_subface_values(
   Assert((dynamic_cast<const InternalData *>(&internal_data) != nullptr),
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
+  data.output_data         = &output_data;
 
   // if necessary, recompute the support points of the transformation of this
   // cell (note that we need to first check the triangulation pointer, since
-  // otherwise the second test might trigger an exception if the triangulations
-  // are not the same)
+  // otherwise the second test might trigger an exception if the
+  // triangulations are not the same)
   if ((data.mapping_support_points.size() == 0) ||
       (&cell->get_triangulation() !=
        &data.cell_of_current_support_points->get_triangulation()) ||
@@ -1245,7 +1095,6 @@ MappingQ<dim, spacedim>::fill_fe_subface_values(
     quadrature,
     data,
     polynomials_1d,
-    polynomial_degree,
     renumber_lexicographic_to_hierarchic,
     output_data);
 }
@@ -1267,6 +1116,7 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
   Assert(dynamic_cast<const InternalData *>(&internal_data) != nullptr,
          ExcInternalError());
   const InternalData &data = static_cast<const InternalData &>(internal_data);
+  data.output_data         = &output_data;
 
   const unsigned int n_q_points = quadrature.size();
 
@@ -1278,7 +1128,6 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
     data,
     make_array_view(quadrature.get_points()),
     polynomials_1d,
-    polynomial_degree,
     renumber_lexicographic_to_hierarchic,
     output_data.quadrature_points,
     output_data.jacobians,
@@ -1286,43 +1135,55 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
 
   internal::MappingQImplementation::maybe_update_jacobian_grads<dim, spacedim>(
     CellSimilarity::none,
-    QProjector<dim>::DataSetDescriptor::cell(),
     data,
+    make_array_view(quadrature.get_points()),
+    polynomials_1d,
+    renumber_lexicographic_to_hierarchic,
     output_data.jacobian_grads);
 
   internal::MappingQImplementation::maybe_update_jacobian_pushed_forward_grads<
     dim,
     spacedim>(CellSimilarity::none,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_pushed_forward_grads);
 
   internal::MappingQImplementation::maybe_update_jacobian_2nd_derivatives<
     dim,
     spacedim>(CellSimilarity::none,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_2nd_derivatives);
 
   internal::MappingQImplementation::
     maybe_update_jacobian_pushed_forward_2nd_derivatives<dim, spacedim>(
       CellSimilarity::none,
-      QProjector<dim>::DataSetDescriptor::cell(),
       data,
+      make_array_view(quadrature.get_points()),
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic,
       output_data.jacobian_pushed_forward_2nd_derivatives);
 
   internal::MappingQImplementation::maybe_update_jacobian_3rd_derivatives<
     dim,
     spacedim>(CellSimilarity::none,
-              QProjector<dim>::DataSetDescriptor::cell(),
               data,
+              make_array_view(quadrature.get_points()),
+              polynomials_1d,
+              renumber_lexicographic_to_hierarchic,
               output_data.jacobian_3rd_derivatives);
 
   internal::MappingQImplementation::
     maybe_update_jacobian_pushed_forward_3rd_derivatives<dim, spacedim>(
       CellSimilarity::none,
-      QProjector<dim>::DataSetDescriptor::cell(),
       data,
+      make_array_view(quadrature.get_points()),
+      polynomials_1d,
+      renumber_lexicographic_to_hierarchic,
       output_data.jacobian_pushed_forward_3rd_derivatives);
 
   const UpdateFlags          update_flags = data.update_each;
@@ -1340,7 +1201,7 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
 
       for (unsigned int point = 0; point < n_q_points; ++point)
         {
-          const double det = data.contravariant[point].determinant();
+          const double det = data.volume_elements[point];
 
           // check for distorted cells.
 
@@ -1355,8 +1216,8 @@ MappingQ<dim, spacedim>::fill_fe_immersed_surface_values(
           // The normals are n = J^{-T} * \hat{n} before normalizing.
           Tensor<1, spacedim> normal;
           for (unsigned int d = 0; d < spacedim; d++)
-            normal[d] =
-              data.covariant[point][d] * quadrature.normal_vector(point);
+            normal[d] = output_data.inverse_jacobians[point].transpose()[d] *
+                        quadrature.normal_vector(point);
 
           output_data.JxW_values[point] = weights[point] * det * normal.norm();
 
@@ -1395,6 +1256,7 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_generic_points(
                    Quadrature<dim>(std::vector<Point<dim>>(unit_points.begin(),
                                                            unit_points.end())));
   const InternalData &data = static_cast<const InternalData &>(*internal_data);
+  data.output_data         = &output_data;
   data.mapping_support_points = this->compute_mapping_support_points(cell);
 
   internal::MappingQImplementation::maybe_update_q_points_Jacobians_generic(
@@ -1402,7 +1264,6 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_generic_points(
     data,
     unit_points,
     polynomials_1d,
-    polynomial_degree,
     renumber_lexicographic_to_hierarchic,
     output_data.quadrature_points,
     output_data.jacobians,
@@ -1430,6 +1291,7 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_face_quadrature(
   const InternalData &data = static_cast<const InternalData &>(internal_data);
 
   data.mapping_support_points = this->compute_mapping_support_points(cell);
+  data.output_data            = &output_data;
 
   internal::MappingQImplementation::do_fill_fe_face_values(
     *this,
@@ -1440,7 +1302,6 @@ MappingQ<dim, spacedim>::fill_mapping_data_for_face_quadrature(
     face_quadrature,
     data,
     polynomials_1d,
-    polynomial_degree,
     renumber_lexicographic_to_hierarchic,
     output_data);
 }
@@ -1522,13 +1383,14 @@ MappingQ<dim, spacedim>::transform(
   AssertDimension(input.size(), output.size());
   Assert(dynamic_cast<const InternalData *>(&mapping_data) != nullptr,
          ExcInternalError());
-  const InternalData &data = static_cast<const InternalData &>(mapping_data);
+  const internal::FEValuesImplementation::MappingRelatedData<dim, spacedim>
+    &data = *static_cast<const InternalData &>(mapping_data).output_data;
 
   switch (mapping_kind)
     {
       case mapping_covariant_gradient:
         {
-          Assert(data.update_each & update_contravariant_transformation,
+          Assert(!data.inverse_jacobians.empty(),
                  typename FEValuesBase<dim>::ExcAccessToUninitializedField(
                    "update_covariant_transformation"));
 
@@ -1536,18 +1398,20 @@ MappingQ<dim, spacedim>::transform(
             for (unsigned int i = 0; i < spacedim; ++i)
               for (unsigned int j = 0; j < spacedim; ++j)
                 {
-                  double tmp[dim];
+                  double                                 tmp[dim];
+                  const DerivativeForm<1, dim, spacedim> covariant =
+                    data.inverse_jacobians[q].transpose();
                   for (unsigned int K = 0; K < dim; ++K)
                     {
-                      tmp[K] = data.covariant[q][j][0] * input[q][i][0][K];
+                      tmp[K] = covariant[j][0] * input[q][i][0][K];
                       for (unsigned int J = 1; J < dim; ++J)
-                        tmp[K] += data.covariant[q][j][J] * input[q][i][J][K];
+                        tmp[K] += covariant[j][J] * input[q][i][J][K];
                     }
                   for (unsigned int k = 0; k < spacedim; ++k)
                     {
-                      output[q][i][j][k] = data.covariant[q][k][0] * tmp[0];
+                      output[q][i][j][k] = covariant[k][0] * tmp[0];
                       for (unsigned int K = 1; K < dim; ++K)
-                        output[q][i][j][k] += data.covariant[q][k][K] * tmp[K];
+                        output[q][i][j][k] += covariant[k][K] * tmp[K];
                     }
                 }
           return;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.