]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Split out mapping-related part from FEEvaluationBase into base class 10901/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Thu, 10 Sep 2020 14:57:31 +0000 (16:57 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sat, 12 Sep 2020 20:09:50 +0000 (22:09 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h

index 19055df97d0c7b48ea82314941957e1154b54a11..57513e8977547a2b1c6bca9bbbdc6eb225d61bee 100644 (file)
 DEAL_II_NAMESPACE_OPEN
 
 
+// forward declaration
+template <int, typename, bool, typename>
+class FEEvaluationBaseData;
+
+
 
 namespace internal
 {
@@ -3165,12 +3170,14 @@ namespace internal
   template <int dim, int fe_degree, typename Number>
   struct CellwiseInverseMassMatrixImpl
   {
-    template <typename FEEvaluationType>
     static void
-    apply(const unsigned int      n_components,
-          const FEEvaluationType &fe_eval,
-          const Number *          in_array,
-          Number *                out_array)
+    apply(const unsigned int                  n_components,
+          const FEEvaluationBaseData<dim,
+                                     typename Number::value_type,
+                                     false,
+                                     Number> &fe_eval,
+          const Number *                      in_array,
+          Number *                            out_array)
     {
       constexpr unsigned int dofs_per_component =
         Utilities::pow(fe_degree + 1, dim);
index 4e4d7db7270008640552e4be32ff32573b13466d..38bb982967a2d4fd0dab81f2d5496647191733a7 100644 (file)
@@ -60,26 +60,22 @@ class FEEvaluation;
 
 
 /**
- * This is the base class for the FEEvaluation classes. This class is a base
- * class and needs usually not be called in user code. It does not have any
- * public constructor. The usage is through the class FEEvaluation instead. It
- * implements a reinit method that is used to set pointers so that operations
- * on quadrature points can be performed quickly, access functions to vectors
- * for the @p read_dof_values, @p set_dof_values, and @p
- * distribute_local_to_global functions, as well as methods to access values
- * and gradients of finite element functions.
+ * This base class of the FEEvaluation and FEFaceEvaluation classes handles
+ * mapping-related information independent of the degrees of freedom and
+ * finite element in use. This class provides access functionality for user
+ * code but is otherwise invisible without any public constructor. The usage
+ * is through the class FEEvaluation instead.
  *
- * This class has three template arguments:
+ * This class has four template arguments:
  *
  * @tparam dim Dimension in which this class is to be used
  *
- * @tparam n_components Number of vector components when solving a system of
- * PDEs. If the same operation is applied to several components of a PDE (e.g.
- * a vector Laplace equation), they can be applied simultaneously with one
- * call (and often more efficiently)
- *
  * @tparam Number Number format, usually @p double or @p float
  *
+ * @tparam is_face Whether the class is used for a cell integrator (with
+ * quadrature dimension the same as the space dimension) or for a face
+ * integrator (with quadrature dimension one less)
+ *
  * @tparam VectorizedArrayType Type of array to be woked on in a vectorized
  *                             fashion, defaults to VectorizedArray<Number>
  *
@@ -90,32 +86,22 @@ class FEEvaluation;
  * @ingroup matrixfree
  */
 template <int dim,
-          int n_components_,
           typename Number,
           bool is_face                 = false,
           typename VectorizedArrayType = VectorizedArray<Number>>
-class FEEvaluationBase
+class FEEvaluationBaseData
 {
   static_assert(
     std::is_same<Number, typename VectorizedArrayType::value_type>::value,
     "Type of Number and of VectorizedArrayType do not match.");
 
 public:
-  using number_type = Number;
-  using value_type  = Tensor<1, n_components_, VectorizedArrayType>;
-  using gradient_type =
-    Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>;
-  static constexpr unsigned int dimension    = dim;
-  static constexpr unsigned int n_components = n_components_;
+  static constexpr unsigned int dimension = dim;
 
-  /**
-   * @name 1: General operations
-   */
-  //@{
   /**
    * Destructor.
    */
-  ~FEEvaluationBase();
+  ~FEEvaluationBaseData();
 
   /**
    * Return the index offset within the geometry fields for the cell the @p
@@ -143,10 +129,365 @@ public:
   const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
   get_shape_info() const;
 
-  //@}
+  /**
+   * Return the determinant of the Jacobian from the unit to the real cell
+   * times the quadrature weight.
+   */
+  VectorizedArrayType
+  JxW(const unsigned int q_point) const;
+
+  /**
+   * Return the inverse and transposed version of Jacobian of the mapping
+   * between the unit to the real cell (representing the covariant
+   * transformation). This is exactly the matrix used internally to transform
+   * the unit cell gradients to gradients on the real cell.
+   */
+  Tensor<2, dim, VectorizedArrayType>
+  inverse_jacobian(const unsigned int q_point) const;
+
+  /**
+   * Return the unit normal vector on a face. Note that both sides of a face
+   * use the same orientation of the normal vector: For the faces enumerated
+   * as `interior` in FaceToCellTopology and selected with the
+   * `is_interior_face=true` flag of the constructor, this corresponds to the
+   * outer normal vector, whereas for faces enumerated as `exterior` in
+   * FaceToCellTopology and selected with the `is_interior_face=false` flag of
+   * the constructor, the normal points into the element as a consequence of
+   * the single normal vector.
+   *
+   * @note Only implemented in case `is_face == true`.
+   */
+  Tensor<1, dim, VectorizedArrayType>
+  get_normal_vector(const unsigned int q_point) const;
+
+  /**
+   * Provides a unified interface to access data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_macro_cells() +
+   * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces
+   * (indirect addressing).
+   */
+  VectorizedArrayType
+  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
+
+  /**
+   * Provides a unified interface to set data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_macro_cells() +
+   * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces
+   * (indirect addressing).
+   */
+  void
+  set_cell_data(AlignedVector<VectorizedArrayType> &array,
+                const VectorizedArrayType &         value) const;
+
+  /**
+   * The same as above, just for std::array of length of VectorizedArrayType for
+   * arbitrary data type.
+   */
+  template <typename T>
+  std::array<T, VectorizedArrayType::size()>
+  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
+                   &array) const;
+
+  /**
+   * The same as above, just for std::array of length of VectorizedArrayType for
+   * arbitrary data type.
+   */
+  template <typename T>
+  void
+  set_cell_data(
+    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+    const std::array<T, VectorizedArrayType::size()> &         value) const;
+
+  /**
+   * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
+   * associated with.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()>
+  get_cell_ids() const;
+
+  /**
+   * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
+   * associated with.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()>
+  get_cell_or_face_ids() const;
+
+
+  /**
+   * Return the numbering of local degrees of freedom within the evaluation
+   * routines of FEEvaluation in terms of the standard numbering on finite
+   * elements.
+   */
+  const std::vector<unsigned int> &
+  get_internal_dof_numbering() const;
+
+  /**
+   * Return an ArrayView to internal memory for temporary use. Note that some
+   * of this memory is overwritten during evaluate() and integrate() calls so
+   * do not assume it to be stable over those calls. The maximum size you can
+   * write into is 3*dofs_per_cell+2*n_q_points.
+   */
+  ArrayView<VectorizedArrayType>
+  get_scratch_data() const;
+
+protected:
+  /**
+   * Constructor. Made protected to prevent users from directly using this
+   * class. Takes all data stored in MatrixFree. If applied to problems with
+   * more than one quadrature formula selected during construction of
+   * `matrix_free`, `quad_no` allows to select the appropriate formula.
+   */
+  FEEvaluationBaseData(
+    const MatrixFree<dim, Number, VectorizedArrayType> &matrix_free,
+    const unsigned int                                  dof_no,
+    const unsigned int first_selected_component,
+    const unsigned int quad_no,
+    const unsigned int fe_degree,
+    const unsigned int n_q_points,
+    const bool         is_interior_face);
+
+  /**
+   * Constructor that comes with reduced functionality and works similar as
+   * FEValues.
+   */
+  FEEvaluationBaseData(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other);
+
+  /**
+   * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
+   * quadrature, and update flags, the underlying geometry evaluation based on
+   * FEValues will be deep-copied in order to allow for using in parallel with
+   * threads.
+   */
+  FEEvaluationBaseData(const FEEvaluationBaseData &other);
+
+  /**
+   * Copy assignment operator. If FEEvaluationBase was constructed from a
+   * mapping, fe, quadrature, and update flags, the underlying geometry
+   * evaluation based on FEValues will be deep-copied in order to allow for
+   * using in parallel with threads.
+   */
+  FEEvaluationBaseData &
+  operator=(const FEEvaluationBaseData &other);
+
+  /**
+   * This is the general array for all data fields.
+   */
+  AlignedVector<VectorizedArrayType> *scratch_data_array;
+
+  /**
+   * This is the user-visible part of scratch_data_array, only showing the
+   * last part of scratch_data_array. The first part is consumed by
+   * values_dofs, values_quad, etc.
+   */
+  VectorizedArrayType *scratch_data;
+
+  /**
+   * The number of the quadrature formula of the present cell.
+   */
+  const unsigned int quad_no;
 
   /**
-   * @name 2: Reading from and writing to vectors
+   * The active fe index for this class for efficient indexing in the hp case.
+   */
+  const unsigned int active_fe_index;
+
+  /**
+   * The active quadrature index for this class for efficient indexing in the
+   * hp case.
+   */
+  const unsigned int active_quad_index;
+
+  /**
+   * The number of quadrature points in the current evaluation context.
+   */
+  const unsigned int n_quadrature_points;
+
+  /**
+   * A pointer to the underlying data.
+   */
+  const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
+
+  /**
+   * A pointer to the underlying DoF indices and constraint description
+   * for the component specified at construction. Also contained in
+   * matrix_info, but it simplifies code if we store a reference to it.
+   */
+  const internal::MatrixFreeFunctions::DoFInfo *dof_info;
+
+  /**
+   * A pointer to the underlying transformation data from unit to real cells
+   * for the given quadrature formula specified at construction. Also
+   * contained in matrix_info, but it simplifies code if we store a reference
+   * to it.
+   */
+  const internal::MatrixFreeFunctions::MappingInfoStorage<
+    (is_face ? dim - 1 : dim),
+    dim,
+    Number,
+    VectorizedArrayType> *mapping_data;
+
+  /**
+   * A pointer to the unit cell shape data, i.e., values, gradients and
+   * Hessians in 1D at the quadrature points that constitute the tensor
+   * product. Also contained in matrix_info, but it simplifies code if we
+   * store a reference to it.
+   */
+  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data;
+
+  /**
+   * A pointer to the Jacobian information of the present cell. Only set to a
+   * useful value if on a non-Cartesian cell.
+   */
+  const Tensor<2, dim, VectorizedArrayType> *jacobian;
+
+  /**
+   * A pointer to the Jacobian determinant of the present cell. If on a
+   * Cartesian cell or on a cell with constant Jacobian, this is just the
+   * Jacobian determinant, otherwise the Jacobian determinant times the
+   * quadrature weight.
+   */
+  const VectorizedArrayType *J_value;
+
+  /**
+   * A pointer to the normal vectors at faces.
+   */
+  const Tensor<1, dim, VectorizedArrayType> *normal_vectors;
+
+  /**
+   * A pointer to the normal vectors times the jacobian at faces.
+   */
+  const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian;
+
+  /**
+   * A pointer to the quadrature weights of the underlying quadrature formula.
+   */
+  const Number *quadrature_weights;
+
+  /**
+   * After a call to reinit(), stores the number of the cell we are currently
+   * working with.
+   */
+  unsigned int cell;
+
+  /**
+   * Flag holding information whether a face is an interior or exterior face
+   * according to the defined direction of the normal.  Not used for cells.
+   */
+  bool is_interior_face;
+
+  /**
+   * Stores the index an FEFaceEvaluation object is currently pointing into
+   * (interior face, exterior face, data associated with cell).
+   */
+  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
+
+  /**
+   * Stores the current number of a face within the given cell in case
+   * `is_face==true`, using values between `0` and `2*dim`.
+   */
+  unsigned int face_no;
+
+  /**
+   * Stores the orientation of the given face with respect to the standard
+   * orientation, 0 if in standard orientation.
+   */
+  unsigned int face_orientation;
+
+  /**
+   * Stores the subface index of the given face. Usually, this variable takes
+   * the value numbers::invalid_unsigned_int to indicate integration over the
+   * full face, but in case the current physical face has a neighbor that is
+   * more refined, it is a subface and must scale the entries in ShapeInfo
+   * appropriately.
+   */
+  unsigned int subface_index;
+
+  /**
+   * Stores the type of the cell we are currently working with after a call to
+   * reinit(). Valid values are @p cartesian, @p affine and @p general, which
+   * have different implications on how the Jacobian transformations are
+   * stored internally in MappingInfo.
+   */
+  internal::MatrixFreeFunctions::GeometryType cell_type;
+
+  /**
+   * Geometry data that can be generated FEValues on the fly with the
+   * respective constructor.
+   */
+  std::shared_ptr<internal::MatrixFreeFunctions::
+                    MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
+    mapped_geometry;
+
+  // Make FEEvaluation objects friends for access to protected member
+  // mapped_geometry.
+  template <int, int, int, int, typename, typename>
+  friend class FEEvaluation;
+};
+
+
+
+/**
+ * This is the base class for the FEEvaluation classes.  This class needs
+ * usually not be called in user code and does not have any public
+ * constructor. The usage is through the class FEEvaluation instead. It
+ * implements a reinit method that is used to set pointers so that operations
+ * on quadrature points can be performed quickly, access functions to vectors
+ * for the FEEvaluationBase::read_dof_values(),
+ * FEEvaluationBase::set_dof_values(), and
+ * FEEvaluationBase::distribute_local_to_global() functions, as well as
+ * methods to access values and gradients of finite element functions. It also
+ * inherits the geometry access functions provided by the class
+ * FEEvaluationBaseData.
+ *
+ * This class has five template arguments:
+ *
+ * @tparam dim Dimension in which this class is to be used
+ *
+ * @tparam n_components Number of vector components when solving a system of
+ * PDEs. If the same operation is applied to several components of a PDE (e.g.
+ * a vector Laplace equation), they can be applied simultaneously with one
+ * call (and often more efficiently)
+ *
+ * @tparam Number Number format, usually @p double or @p float
+ *
+ * @tparam is_face Whether the class is used for a cell integrator (with
+ * quadrature dimension the same as the space dimension) or for a face
+ * integrator (with quadrature dimension one less)
+ *
+ * @tparam VectorizedArrayType Type of array to be woked on in a vectorized
+ *                             fashion, defaults to VectorizedArray<Number>
+ *
+ * @note Currently only VectorizedArray<Number, width> is supported as
+ *       VectorizedArrayType.
+ *
+ *
+ * @ingroup matrixfree
+ */
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face                 = false,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+class FEEvaluationBase
+  : public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+{
+public:
+  using number_type = Number;
+  using value_type  = Tensor<1, n_components_, VectorizedArrayType>;
+  using gradient_type =
+    Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>>;
+  static constexpr unsigned int dimension    = dim;
+  static constexpr unsigned int n_components = n_components_;
+
+  /**
+   * @name 1: Reading from and writing to vectors
    */
   //@{
   /**
@@ -305,7 +646,7 @@ public:
   //@}
 
   /**
-   * @name 3: Data access
+   * @name 2: Access to data at quadrature points or the gather vector data
    */
   //@{
   /**
@@ -480,7 +821,7 @@ public:
 #ifdef DOXYGEN
   // doxygen does not anyhow mention functions coming from partial template
   // specialization of the base class, in this case FEEvaluationAccess<dim,dim>.
-  // For now, hack-in those functions manually only to fix documentation:
+  // For now, hack in those functions manually only to fix documentation:
 
   /**
    * Return the divergence of a vector-valued finite element at quadrature
@@ -570,7 +911,7 @@ public:
 #endif
 
   /**
-   * Takes values on quadrature points, multiplies by the Jacobian determinant
+   * Takes values at quadrature points, multiplies by the Jacobian determinant
    * and quadrature weights (JxW) and sums the values for all quadrature
    * points on the cell. The result is a scalar, representing the integral
    * over the function over the cell. If a vector-element is used, the
@@ -580,102 +921,18 @@ public:
    *
    * @note In case the FEEvaluation object is initialized with a batch of
    * cells where not all lanes in the SIMD vector VectorizedArray are
-   * representing actual data, this method performs computations on dummy data
-   * (that is copied from the last valid lane) that will not make sense. Thus,
-   * the user needs to make sure that it is not used in any computation
-   * explicitly, like when summing the results of several cells.
-   */
-  value_type
-  integrate_value() const;
-
-  /**
-   * Return the determinant of the Jacobian from the unit to the real cell
-   * times the quadrature weight.
-   */
-  VectorizedArrayType
-  JxW(const unsigned int q_point) const;
-
-  /**
-   * Return the inverse and transposed version of Jacobian of the mapping
-   * between the unit to the real cell (representing the covariant
-   * transformation). This is exactly the matrix used internally to transform
-   * the unit cell gradients to gradients on the real cell.
-   */
-  Tensor<2, dim, VectorizedArrayType>
-  inverse_jacobian(const unsigned int q_point) const;
-
-  /**
-   * Return the unit normal vector on a face. Note that both sides of a face
-   * use the same orientation of the normal vector: For the faces enumerated
-   * as `interior` in FaceToCellTopology and selected with the
-   * `is_interior_face=true` flag of the constructor, this corresponds to the
-   * outer normal vector, whereas for faces enumerated as `exterior` in
-   * FaceToCellTopology and selected with the `is_interior_face=false` flag of
-   * the constructor, the normal points into the element as a consequence of
-   * the single normal vector.
-   *
-   * @note Only implemented in case `is_face == true`.
-   */
-  Tensor<1, dim, VectorizedArrayType>
-  get_normal_vector(const unsigned int q_point) const;
-
-  /**
-   * Provides a unified interface to access data in a vector of
-   * VectorizedArray fields of length MatrixFree::n_macro_cells() +
-   * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces
-   * (indirect addressing).
-   */
-  VectorizedArrayType
-  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const;
-
-  /**
-   * Provides a unified interface to set data in a vector of
-   * VectorizedArray fields of length MatrixFree::n_macro_cells() +
-   * MatrixFree::n_macro_ghost_cells() for both cells (plain read) and faces
-   * (indirect addressing).
-   */
-  void
-  set_cell_data(AlignedVector<VectorizedArrayType> &array,
-                const VectorizedArrayType &         value) const;
-
-  /**
-   * The same as above, just for std::array of length of VectorizedArrayType for
-   * arbitrary data type.
-   */
-  template <typename T>
-  std::array<T, VectorizedArrayType::size()>
-  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
-                   &array) const;
-
-  /**
-   * The same as above, just for std::array of length of VectorizedArrayType for
-   * arbitrary data type.
-   */
-  template <typename T>
-  void
-  set_cell_data(
-    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
-    const std::array<T, VectorizedArrayType::size()> &         value) const;
-
-  /**
-   * Return the id of the cells this FEEvaluation or FEFaceEvaluation is
-   * associated with.
-   */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  get_cell_ids() const;
-
-
-  /**
-   * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
-   * associated with.
+   * representing actual data, this method performs computations on dummy data
+   * (that is copied from the last valid lane) that will not make sense. Thus,
+   * the user needs to make sure that it is not used in any computation
+   * explicitly, like when summing the results of several cells.
    */
-  std::array<unsigned int, VectorizedArrayType::size()>
-  get_cell_or_face_ids() const;
+  value_type
+  integrate_value() const;
 
   //@}
 
   /**
-   * @name 4: Access to internal data
+   * @name 3: Access to internal data
    */
   //@{
   /**
@@ -784,23 +1041,6 @@ public:
   VectorizedArrayType *
   begin_hessians();
 
-  /**
-   * Return the numbering of local degrees of freedom within the evaluation
-   * routines of FEEvaluation in terms of the standard numbering on finite
-   * elements.
-   */
-  const std::vector<unsigned int> &
-  get_internal_dof_numbering() const;
-
-  /**
-   * Return an ArrayView to internal memory for temporary use. Note that some
-   * of this memory is overwritten during evaluate() and integrate() calls so
-   * do not assume it to be stable over those calls. The maximum size you can
-   * write into is 3*dofs_per_cell+2*n_q_points.
-   */
-  ArrayView<VectorizedArrayType>
-  get_scratch_data() const;
-
   //@}
 
 protected:
@@ -847,25 +1087,23 @@ protected:
    * possibly within the element if the evaluate/integrate routines are
    * combined inside user code (e.g. for computing cell matrices).
    *
-   * The optional FEEvaluationBase object allows several FEEvaluation objects
-   * to share the geometry evaluation, i.e., the underlying mapping and
-   * quadrature points do only need to be evaluated once. This only works if
-   * the quadrature formulas are the same. Otherwise, a new evaluation object
-   * is created. Make sure to not pass an optional object around when you
-   * intend to use the FEEvaluation object in %parallel with another one
-   * because otherwise the intended sharing may create race conditions.
-   */
-  template <int n_components_other>
-  FEEvaluationBase(const Mapping<dim> &      mapping,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<1> &     quadrature,
-                   const UpdateFlags         update_flags,
-                   const unsigned int        first_selected_component,
-                   const FEEvaluationBase<dim,
-                                          n_components_other,
-                                          Number,
-                                          is_face,
-                                          VectorizedArrayType> *other);
+   * The optional FEEvaluationBaseData object allows several
+   * FEEvaluation objects to share the geometry evaluation, i.e., the
+   * underlying mapping and quadrature points do only need to be evaluated
+   * once. This only works if the quadrature formulas are the same. Otherwise,
+   * a new evaluation object is created. Make sure to not pass an optional
+   * object around when you intend to use the FEEvaluation object in %parallel
+   * with another one because otherwise the intended sharing may create race
+   * conditions.
+   */
+  FEEvaluationBase(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other);
 
   /**
    * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
@@ -923,18 +1161,6 @@ protected:
   read_write_operation_global(const VectorOperation &operation,
                               VectorType *           vectors[]) const;
 
-  /**
-   * This is the general array for all data fields.
-   */
-  AlignedVector<VectorizedArrayType> *scratch_data_array;
-
-  /**
-   * This is the user-visible part of scratch_data_array, only showing the
-   * last part of scratch_data_array. The first part is consumed by
-   * values_dofs, values_quad, etc.
-   */
-  VectorizedArrayType *scratch_data;
-
   /**
    * This field stores the values for local degrees of freedom (e.g. after
    * reading out from a vector but before applying unit cell transformations
@@ -990,142 +1216,12 @@ protected:
    */
   VectorizedArrayType *hessians_quad;
 
-  /**
-   * Stores the number of the quadrature formula of the present cell.
-   */
-  const unsigned int quad_no;
-
   /**
    * Stores the number of components in the finite element as detected in the
    * MatrixFree storage class for comparison with the template argument.
    */
   const unsigned int n_fe_components;
 
-  /**
-   * Stores the active fe index for this class for efficient indexing in the
-   * hp case.
-   */
-  const unsigned int active_fe_index;
-
-  /**
-   * Stores the active quadrature index for this class for efficient indexing
-   * in the hp case.
-   */
-  const unsigned int active_quad_index;
-
-  /**
-   * Stores the number of quadrature points in the current evaluation context.
-   */
-  const unsigned int n_quadrature_points;
-
-  /**
-   * Stores a pointer to the underlying data.
-   */
-  const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
-
-  /**
-   * Stores a pointer to the underlying DoF indices and constraint description
-   * for the component specified at construction. Also contained in
-   * matrix_info, but it simplifies code if we store a reference to it.
-   */
-  const internal::MatrixFreeFunctions::DoFInfo *dof_info;
-
-  /**
-   * Stores a pointer to the underlying transformation data from unit to
-   * real cells for the given quadrature formula specified at construction.
-   * Also contained in matrix_info, but it simplifies code if we store a
-   * reference to it.
-   */
-  const internal::MatrixFreeFunctions::MappingInfoStorage<
-    (is_face ? dim - 1 : dim),
-    dim,
-    Number,
-    VectorizedArrayType> *mapping_data;
-
-  /**
-   * Stores a pointer to the unit cell shape data, i.e., values, gradients and
-   * Hessians in 1D at the quadrature points that constitute the tensor
-   * product. Also contained in matrix_info, but it simplifies code if we
-   * store a reference to it.
-   */
-  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> *data;
-
-  /**
-   * A pointer to the Jacobian information of the present cell. Only set to a
-   * useful value if on a non-Cartesian cell.
-   */
-  const Tensor<2, dim, VectorizedArrayType> *jacobian;
-
-  /**
-   * A pointer to the Jacobian determinant of the present cell. If on a
-   * Cartesian cell or on a cell with constant Jacobian, this is just the
-   * Jacobian determinant, otherwise the Jacobian determinant times the
-   * quadrature weight.
-   */
-  const VectorizedArrayType *J_value;
-
-  /**
-   * A pointer to the normal vectors at faces.
-   */
-  const Tensor<1, dim, VectorizedArrayType> *normal_vectors;
-
-  /**
-   * A pointer to the normal vectors times the jacobian at faces.
-   */
-  const Tensor<1, dim, VectorizedArrayType> *normal_x_jacobian;
-
-  /**
-   * A pointer to the quadrature weights of the underlying quadrature formula.
-   */
-  const Number *quadrature_weights;
-
-  /**
-   * After a call to reinit(), stores the number of the cell we are currently
-   * working with.
-   */
-  unsigned int cell;
-
-  /**
-   * Flag holding information whether a face is an interior or exterior face
-   * according to the defined direction of the normal.  Not used for cells.
-   */
-  bool is_interior_face;
-
-  /**
-   * Stores the index an FEFaceEvaluation object is currently pointing into
-   * (interior face, exterior face, data associated with cell).
-   */
-  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
-
-  /**
-   * Stores the current number of a face within the given cell in case
-   * `is_face==true`, using values between `0` and `2*dim`.
-   */
-  unsigned int face_no;
-
-  /**
-   * Stores the orientation of the given face with respect to the standard
-   * orientation, 0 if in standard orientation.
-   */
-  unsigned int face_orientation;
-
-  /**
-   * Stores the subface index of the given face. Usually, this variable takes
-   * the value numbers::invalid_unsigned_int to indicate integration over the
-   * full face, but in case the current physical face has a neighbor that is
-   * more refined, it is a subface and must scale the entries in ShapeInfo
-   * appropriately.
-   */
-  unsigned int subface_index;
-
-  /**
-   * Stores the type of the cell we are currently working with after a call to
-   * reinit(). Valid values are @p cartesian, @p affine and @p general, which
-   * have different implications on how the Jacobian transformations are
-   * stored internally in MappingInfo.
-   */
-  internal::MatrixFreeFunctions::GeometryType cell_type;
-
   /**
    * Debug information to track whether dof values have been initialized
    * before accessed. Used to control exceptions when uninitialized data is
@@ -1168,14 +1264,6 @@ protected:
    */
   bool gradients_quad_submitted;
 
-  /**
-   * Geometry data that can be generated FEValues on the fly with the
-   * respective constructor.
-   */
-  std::shared_ptr<internal::MatrixFreeFunctions::
-                    MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
-    mapped_geometry;
-
   /**
    * For a FiniteElement with more than one base element, select at which
    * component this data structure should start.
@@ -1191,16 +1279,10 @@ protected:
 private:
   /**
    * Sets the pointers for values, gradients, hessians to the central
-   * scratch_data_array.
+   * scratch_data_array of the base class.
    */
   void
   set_data_pointers();
-
-  // Make other FEEvaluationBase as well as FEEvaluation objects friends.
-  template <int, int, typename, bool, typename>
-  friend class FEEvaluationBase;
-  template <int, int, int, int, typename, typename>
-  friend class FEEvaluation;
 };
 
 
@@ -1258,17 +1340,14 @@ protected:
    * Constructor with reduced functionality for similar usage of FEEvaluation
    * as FEValues, including matrix assembly.
    */
-  template <int n_components_other>
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other);
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other);
 
   /**
    * Copy constructor
@@ -1396,17 +1475,14 @@ protected:
    * Constructor with reduced functionality for similar usage of FEEvaluation
    * as FEValues, including matrix assembly.
    */
-  template <int n_components_other>
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other);
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other);
 
   /**
    * Copy constructor
@@ -1558,17 +1634,14 @@ protected:
    * Constructor with reduced functionality for similar usage of FEEvaluation
    * as FEValues, including matrix assembly.
    */
-  template <int n_components_other>
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other);
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other);
 
   /**
    * Copy constructor
@@ -1707,17 +1780,13 @@ protected:
    * Constructor with reduced functionality for similar usage of FEEvaluation
    * as FEValues, including matrix assembly.
    */
-  template <int n_components_other>
-  FEEvaluationAccess(const Mapping<1> &      mapping,
-                     const FiniteElement<1> &fe,
-                     const Quadrature<1> &   quadrature,
-                     const UpdateFlags       update_flags,
-                     const unsigned int      first_selected_component,
-                     const FEEvaluationBase<1,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other);
+  FEEvaluationAccess(
+    const Mapping<1> &      mapping,
+    const FiniteElement<1> &fe,
+    const Quadrature<1> &   quadrature,
+    const UpdateFlags       update_flags,
+    const unsigned int      first_selected_component,
+    const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other);
 
   /**
    * Copy constructor
@@ -2450,14 +2519,10 @@ public:
    * use the FEEvaluation object in %parallel to the given one because
    * otherwise the intended sharing may create race conditions.
    */
-  template <int n_components_other>
-  FEEvaluation(const FiniteElement<dim> &                   fe,
-               const FEEvaluationBase<dim,
-                                      n_components_other,
-                                      Number,
-                                      false,
-                                      VectorizedArrayType> &other,
-               const unsigned int first_selected_component = 0);
+  FEEvaluation(
+    const FiniteElement<dim> &                                           fe,
+    const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &other,
+    const unsigned int first_selected_component = 0);
 
   /**
    * Copy constructor. If FEEvaluationBase was constructed from a mapping, fe,
@@ -3107,29 +3172,20 @@ namespace internal
 #ifndef DOXYGEN
 
 
+/*----------------------- FEEvaluationBaseData ------------------------*/
 
-/*----------------------- FEEvaluationBase ----------------------------------*/
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType>::
-  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
-                   const unsigned int                                  dof_no,
-                   const unsigned int first_selected_component,
-                   const unsigned int quad_no_in,
-                   const unsigned int fe_degree,
-                   const unsigned int n_q_points,
-                   const bool         is_interior_face)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(
+    const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
+    const unsigned int                                  dof_no,
+    const unsigned int first_selected_component,
+    const unsigned int quad_no_in,
+    const unsigned int fe_degree,
+    const unsigned int n_q_points,
+    const bool         is_interior_face)
   : scratch_data_array(data_in.acquire_scratch_data())
   , quad_no(quad_no_in)
-  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
   , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
                       data_in.get_dof_info(dof_no).fe_index_from_degree(
                         first_selected_component,
@@ -3185,15 +3241,7 @@ inline FEEvaluationBase<dim,
            internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
   , cell_type(internal::MatrixFreeFunctions::general)
-  , dof_values_initialized(false)
-  , values_quad_initialized(false)
-  , gradients_quad_initialized(false)
-  , hessians_quad_initialized(false)
-  , values_quad_submitted(false)
-  , gradients_quad_submitted(false)
-  , first_selected_component(first_selected_component)
 {
-  set_data_pointers();
   Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
   AssertDimension(matrix_info->get_task_info().vectorization_length,
                   VectorizedArrayType::size());
@@ -3201,59 +3249,22 @@ inline FEEvaluationBase<dim,
                   n_quadrature_points);
   AssertDimension(n_quadrature_points,
                   mapping_data->descriptor[active_quad_index].n_q_points);
-  Assert(
-    dof_info->start_components.back() == 1 ||
-      static_cast<int>(n_components_) <=
-        static_cast<int>(
-          dof_info->start_components
-            [dof_info->component_to_base_index[first_selected_component] + 1]) -
-          first_selected_component,
-    ExcMessage(
-      "You tried to construct a vector-valued evaluator with " +
-      std::to_string(n_components) +
-      " components. However, "
-      "the current base element has only " +
-      std::to_string(
-        dof_info->start_components
-          [dof_info->component_to_base_index[first_selected_component] + 1] -
-        first_selected_component) +
-      " components left when starting from local element index " +
-      std::to_string(
-        first_selected_component -
-        dof_info->start_components
-          [dof_info->component_to_base_index[first_selected_component]]) +
-      " (global index " + std::to_string(first_selected_component) + ")"));
-
-  // do not check for correct dimensions of data fields here, should be done
-  // in derived classes
 }
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-template <int n_components_other>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType>::
-  FEEvaluationBase(const Mapping<dim> &      mapping,
-                   const FiniteElement<dim> &fe,
-                   const Quadrature<1> &     quadrature,
-                   const UpdateFlags         update_flags,
-                   const unsigned int        first_selected_component,
-                   const FEEvaluationBase<dim,
-                                          n_components_other,
-                                          Number,
-                                          is_face,
-                                          VectorizedArrayType> *other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other)
   : scratch_data_array(new AlignedVector<VectorizedArrayType>())
   , quad_no(numbers::invalid_unsigned_int)
-  , n_fe_components(n_components_)
   , active_fe_index(numbers::invalid_unsigned_int)
   , active_quad_index(numbers::invalid_unsigned_int)
   , n_quadrature_points(
@@ -3271,24 +3282,12 @@ inline FEEvaluationBase<dim,
   , J_value(nullptr)
   , normal_vectors(nullptr)
   , normal_x_jacobian(nullptr)
-  , quadrature_weights(nullptr)
-  , cell(0)
-  , cell_type(internal::MatrixFreeFunctions::general)
-  , is_interior_face(true)
-  , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
-  , dof_values_initialized(false)
-  , values_quad_initialized(false)
-  , gradients_quad_initialized(false)
-  , hessians_quad_initialized(false)
-  , values_quad_submitted(false)
-  , gradients_quad_submitted(false)
-  ,
-  // keep the number of the selected component within the current base element
-  // for reading dof values
-  first_selected_component(first_selected_component)
+  , quadrature_weights(nullptr)
+  , cell(0)
+  , cell_type(internal::MatrixFreeFunctions::general)
+  , is_interior_face(true)
+  , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
 {
-  set_data_pointers();
-
   Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
          ExcInternalError());
   if (other != nullptr &&
@@ -3304,40 +3303,19 @@ inline FEEvaluationBase<dim,
   mapping_data = &mapped_geometry->get_data_storage();
   jacobian     = mapped_geometry->get_data_storage().jacobians[0].begin();
   J_value      = mapped_geometry->get_data_storage().JxW_values.begin();
-
-  const unsigned int base_element_number =
-    fe.component_to_base_index(first_selected_component).first;
-  Assert(fe.element_multiplicity(base_element_number) == 1 ||
-           fe.element_multiplicity(base_element_number) -
-               first_selected_component >=
-             n_components_,
-         ExcMessage("The underlying element must at least contain as many "
-                    "components as requested by this class"));
-  (void)base_element_number;
 }
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType>::
-  FEEvaluationBase(const FEEvaluationBase<dim,
-                                          n_components_,
-                                          Number,
-                                          is_face,
-                                          VectorizedArrayType> &other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      &other)
   : scratch_data_array(other.matrix_info == nullptr ?
                          new AlignedVector<VectorizedArrayType>() :
                          other.matrix_info->acquire_scratch_data())
   , quad_no(other.quad_no)
-  , n_fe_components(other.n_fe_components)
   , active_fe_index(other.active_fe_index)
   , active_quad_index(other.active_quad_index)
   , n_quadrature_points(other.n_quadrature_points)
@@ -3360,16 +3338,7 @@ inline FEEvaluationBase<dim,
   , cell_type(internal::MatrixFreeFunctions::general)
   , is_interior_face(other.is_interior_face)
   , dof_access_index(other.dof_access_index)
-  , dof_values_initialized(false)
-  , values_quad_initialized(false)
-  , gradients_quad_initialized(false)
-  , hessians_quad_initialized(false)
-  , values_quad_submitted(false)
-  , gradients_quad_submitted(false)
-  , first_selected_component(other.first_selected_component)
 {
-  set_data_pointers();
-
   // Create deep copy of mapped geometry for use in parallel...
   if (other.mapped_geometry.get() != nullptr)
     {
@@ -3389,28 +3358,14 @@ inline FEEvaluationBase<dim,
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType> &
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-operator=(const FEEvaluationBase<dim,
-                                 n_components_,
-                                 Number,
-                                 is_face,
-                                 VectorizedArrayType> &other)
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::operator=(
+  const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &other)
 {
   AssertDimension(quad_no, other.quad_no);
-  AssertDimension(n_fe_components, other.n_fe_components);
   AssertDimension(active_fe_index, other.active_fe_index);
   AssertDimension(active_quad_index, other.active_quad_index);
-  AssertDimension(first_selected_component, other.first_selected_component);
 
   // release old memory
   if (matrix_info == nullptr)
@@ -3437,7 +3392,6 @@ operator=(const FEEvaluationBase<dim,
       data               = other.data;
       scratch_data_array = matrix_info->acquire_scratch_data();
     }
-  set_data_pointers();
 
   quadrature_weights =
     (mapping_data != nullptr ?
@@ -3468,16 +3422,9 @@ operator=(const FEEvaluationBase<dim,
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType>::~FEEvaluationBase()
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  ~FEEvaluationBaseData()
 {
   if (matrix_info != nullptr)
     {
@@ -3499,61 +3446,9 @@ inline FEEvaluationBase<dim,
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_data_pointers()
-{
-  Assert(scratch_data_array != nullptr, ExcInternalError());
-
-  const unsigned int tensor_dofs_per_component =
-    Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
-  const unsigned int dofs_per_component =
-    this->data->dofs_per_component_on_cell;
-  const unsigned int n_quadrature_points =
-    is_face ? this->data->n_q_points_face : this->data->n_q_points;
-
-  const unsigned int shift =
-    std::max(tensor_dofs_per_component + 1, dofs_per_component) *
-      n_components_ * 3 +
-    2 * n_quadrature_points;
-  const unsigned int allocated_size =
-    shift + n_components_ * dofs_per_component +
-    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
-  scratch_data_array->resize_fast(allocated_size);
-
-  // set the pointers to the correct position in the data array
-  for (unsigned int c = 0; c < n_components_; ++c)
-    {
-      this->values_dofs[c] =
-        scratch_data_array->begin() + c * dofs_per_component;
-    }
-  this->values_quad =
-    scratch_data_array->begin() + n_components * dofs_per_component;
-  this->gradients_quad =
-    scratch_data_array->begin() +
-    n_components * (dofs_per_component + n_quadrature_points);
-  this->hessians_quad =
-    scratch_data_array->begin() +
-    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
-  scratch_data =
-    scratch_data_array->begin() + n_components_ * dofs_per_component +
-    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
-}
-
-
-
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline unsigned int
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   get_mapping_data_index_offset() const
 {
   if (matrix_info == nullptr)
@@ -3567,14 +3462,10 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline internal::MatrixFreeFunctions::GeometryType
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  get_cell_type() const
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
+  const
 {
   Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
   return cell_type;
@@ -3582,13 +3473,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   get_shape_info() const
 {
   Assert(data != nullptr, ExcInternalError());
@@ -3597,13 +3484,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
-                             FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+                             FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   get_normal_vector(const unsigned int q_point) const
 {
   AssertIndexRange(q_point, n_quadrature_points);
@@ -3616,13 +3499,9 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
-                             FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::JxW(
+                             FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
   const unsigned int q_point) const
 {
   AssertIndexRange(q_point, n_quadrature_points);
@@ -3638,13 +3517,9 @@ inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
-                             FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+                             FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   inverse_jacobian(const unsigned int q_point) const
 {
   AssertIndexRange(q_point, n_quadrature_points);
@@ -3657,15 +3532,13 @@ inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 inline std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  get_cell_ids() const
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
+  const
 {
+  Assert(this->matrix_info != nullptr, ExcNotInitialized());
+
   const unsigned int                n_lanes = VectorizedArrayType::size();
   std::array<unsigned int, n_lanes> cells;
 
@@ -3737,7 +3610,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 namespace internal
 {
   template <int dim,
-            int n_components_,
             typename Number,
             bool is_face,
             typename VectorizedArrayType,
@@ -3746,11 +3618,7 @@ namespace internal
             typename FU>
   inline void
   process_cell_data(
-    const FEEvaluationBase<dim,
-                           n_components_,
-                           Number,
-                           is_face,
-                           VectorizedArrayType> &       phi,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &phi,
     const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info,
     GlobalVectorType &                                  array,
     VectorizedArrayType2 &                              out,
@@ -3775,13 +3643,9 @@ namespace internal
 
 
 
-template <int dim,
-          int n_components_,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
 std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
   get_cell_or_face_ids() const
 {
   const unsigned int v_len = VectorizedArrayType::size();
@@ -3807,29 +3671,178 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                                                                   this->face_no,
                                                                   i);
 
-          if (fn == numbers::invalid_unsigned_int)
-            continue; // invalid face ID: no neighbor on boundary
+          if (fn == numbers::invalid_unsigned_int)
+            continue; // invalid face ID: no neighbor on boundary
+
+          // get cell ID on both sides of face
+          auto cell_m = this->matrix_info->get_face_info(fn / v_len)
+                          .cells_interior[fn % v_len];
+          auto cell_p = this->matrix_info->get_face_info(fn / v_len)
+                          .cells_exterior[fn % v_len];
+
+          // compare the IDs with the given cell ID
+          if (cell_m == cell_this)
+            cells[i] = cell_p; // neighbor has the other ID
+          else if (cell_p == cell_this)
+            cells[i] = cell_m;
+        }
+    }
+  else
+    {
+      for (unsigned int i = 0; i < v_len; ++i)
+        cells[i] = cell * v_len + i;
+    }
+
+  return cells;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline VectorizedArrayType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::read_cell_data(
+  const AlignedVector<VectorizedArrayType> &array) const
+{
+  VectorizedArrayType out = Number(1.);
+  internal::process_cell_data(
+    *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
+      local = global;
+    });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::set_cell_data(
+  AlignedVector<VectorizedArrayType> &array,
+  const VectorizedArrayType &         in) const
+{
+  internal::process_cell_data(
+    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+      global = local;
+    });
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <typename T>
+inline std::array<T, VectorizedArrayType::size()>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::read_cell_data(
+  const AlignedVector<std::array<T, VectorizedArrayType::size()>> &array) const
+{
+  std::array<T, VectorizedArrayType::size()> out;
+  internal::process_cell_data(
+    *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
+      local = global;
+    });
+  return out;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <typename T>
+inline void
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::set_cell_data(
+  AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+  const std::array<T, VectorizedArrayType::size()> &         in) const
+{
+  internal::process_cell_data(
+    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+      global = local;
+    });
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const std::vector<unsigned int> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_internal_dof_numbering() const
+{
+  return data->lexicographic_numbering;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline ArrayView<VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_scratch_data() const
+{
+  return ArrayView<VectorizedArrayType>(
+    const_cast<VectorizedArrayType *>(scratch_data),
+    scratch_data_array->end() - scratch_data);
+}
+
 
-          // get cell ID on both sides of face
-          auto cell_m = this->matrix_info->get_face_info(fn / v_len)
-                          .cells_interior[fn % v_len];
-          auto cell_p = this->matrix_info->get_face_info(fn / v_len)
-                          .cells_exterior[fn % v_len];
+/*----------------------- FEEvaluationBase ----------------------------------*/
 
-          // compare the IDs with the given cell ID
-          if (cell_m == cell_this)
-            cells[i] = cell_p; // neighbor has the other ID
-          else if (cell_p == cell_this)
-            cells[i] = cell_m;
-        }
-    }
-  else
-    {
-      for (unsigned int i = 0; i < v_len; ++i)
-        cells[i] = cell * v_len + i;
-    }
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType>::
+  FEEvaluationBase(const MatrixFree<dim, Number, VectorizedArrayType> &data_in,
+                   const unsigned int                                  dof_no,
+                   const unsigned int first_selected_component,
+                   const unsigned int quad_no_in,
+                   const unsigned int fe_degree,
+                   const unsigned int n_q_points,
+                   const bool         is_interior_face)
+  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+      data_in,
+      dof_no,
+      first_selected_component,
+      quad_no_in,
+      fe_degree,
+      n_q_points,
+      is_interior_face)
+  , n_fe_components(data_in.get_dof_info(dof_no).start_components.back())
+  , dof_values_initialized(false)
+  , values_quad_initialized(false)
+  , gradients_quad_initialized(false)
+  , hessians_quad_initialized(false)
+  , values_quad_submitted(false)
+  , gradients_quad_submitted(false)
+  , first_selected_component(first_selected_component)
+{
+  set_data_pointers();
+  Assert(
+    this->dof_info->start_components.back() == 1 ||
+      static_cast<int>(n_components_) <=
+        static_cast<int>(
+          this->dof_info->start_components
+            [this->dof_info->component_to_base_index[first_selected_component] +
+             1]) -
+          first_selected_component,
+    ExcMessage(
+      "You tried to construct a vector-valued evaluator with " +
+      std::to_string(n_components) +
+      " components. However, "
+      "the current base element has only " +
+      std::to_string(
+        this->dof_info->start_components
+          [this->dof_info->component_to_base_index[first_selected_component] +
+           1] -
+        first_selected_component) +
+      " components left when starting from local element index " +
+      std::to_string(
+        first_selected_component -
+        this->dof_info->start_components
+          [this->dof_info->component_to_base_index[first_selected_component]]) +
+      " (global index " + std::to_string(first_selected_component) + ")"));
 
-  return cells;
+  // do not check for correct dimensions of data fields here, should be done
+  // in derived classes
 }
 
 
@@ -3839,16 +3852,48 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-inline VectorizedArrayType
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType>::
+  FEEvaluationBase(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other)
+  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+      mapping,
+      fe,
+      quadrature,
+      update_flags,
+      first_selected_component,
+      other)
+  , n_fe_components(n_components_)
+  , dof_values_initialized(false)
+  , values_quad_initialized(false)
+  , gradients_quad_initialized(false)
+  , hessians_quad_initialized(false)
+  , values_quad_submitted(false)
+  , gradients_quad_submitted(false)
+  // keep the number of the selected component within the current base element
+  // for reading dof values
+  , first_selected_component(first_selected_component)
 {
-  VectorizedArrayType out = Number(1.);
-  internal::process_cell_data(
-    *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
-      local = global;
-    });
-  return out;
+  set_data_pointers();
+
+  const unsigned int base_element_number =
+    fe.component_to_base_index(first_selected_component).first;
+  Assert(fe.element_multiplicity(base_element_number) == 1 ||
+           fe.element_multiplicity(base_element_number) -
+               first_selected_component >=
+             n_components_,
+         ExcMessage("The underlying element must at least contain as many "
+                    "components as requested by this class"));
+  (void)base_element_number;
 }
 
 
@@ -3858,15 +3903,27 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-inline void
-FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_cell_data(AlignedVector<VectorizedArrayType> &array,
-                const VectorizedArrayType &         in) const
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType>::
+  FEEvaluationBase(const FEEvaluationBase<dim,
+                                          n_components_,
+                                          Number,
+                                          is_face,
+                                          VectorizedArrayType> &other)
+  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(other)
+  , n_fe_components(other.n_fe_components)
+  , dof_values_initialized(false)
+  , values_quad_initialized(false)
+  , gradients_quad_initialized(false)
+  , hessians_quad_initialized(false)
+  , values_quad_submitted(false)
+  , gradients_quad_submitted(false)
+  , first_selected_component(other.first_selected_component)
 {
-  internal::process_cell_data(
-    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
-      global = local;
-    });
+  set_data_pointers();
 }
 
 
@@ -3876,18 +3933,24 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-template <typename T>
-inline std::array<T, VectorizedArrayType::size()>
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType> &
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
-                   &array) const
+operator=(const FEEvaluationBase<dim,
+                                 n_components_,
+                                 Number,
+                                 is_face,
+                                 VectorizedArrayType> &other)
 {
-  std::array<T, VectorizedArrayType::size()> out;
-  internal::process_cell_data(
-    *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
-      local = global;
-    });
-  return out;
+  this->FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  operator=(other);
+  AssertDimension(n_fe_components, other.n_fe_components);
+  AssertDimension(first_selected_component, other.first_selected_component);
+
+  return *this;
 }
 
 
@@ -3897,17 +3960,44 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-template <typename T>
 inline void
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_cell_data(
-    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
-    const std::array<T, VectorizedArrayType::size()> &         in) const
+  set_data_pointers()
 {
-  internal::process_cell_data(
-    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
-      global = local;
-    });
+  Assert(this->scratch_data_array != nullptr, ExcInternalError());
+
+  const unsigned int tensor_dofs_per_component =
+    Utilities::fixed_power<dim>(this->data->data.front().fe_degree + 1);
+  const unsigned int dofs_per_component =
+    this->data->dofs_per_component_on_cell;
+  const unsigned int n_quadrature_points =
+    is_face ? this->data->n_q_points_face : this->data->n_q_points;
+
+  const unsigned int shift =
+    std::max(tensor_dofs_per_component + 1, dofs_per_component) *
+      n_components_ * 3 +
+    2 * n_quadrature_points;
+  const unsigned int allocated_size =
+    shift + n_components_ * dofs_per_component +
+    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
+  this->scratch_data_array->resize_fast(allocated_size);
+
+  // set the pointers to the correct position in the data array
+  for (unsigned int c = 0; c < n_components_; ++c)
+    {
+      values_dofs[c] =
+        this->scratch_data_array->begin() + c * dofs_per_component;
+    }
+  values_quad =
+    this->scratch_data_array->begin() + n_components * dofs_per_component;
+  gradients_quad = this->scratch_data_array->begin() +
+                   n_components * (dofs_per_component + n_quadrature_points);
+  hessians_quad =
+    this->scratch_data_array->begin() +
+    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
+  this->scratch_data =
+    this->scratch_data_array->begin() + n_components_ * dofs_per_component +
+    (n_components_ * ((dim * (dim + 1)) / 2 + dim + 1) * n_quadrature_points);
 }
 
 
@@ -4006,14 +4096,14 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   // Case 1: No MatrixFree object given, simple case because we do not need to
   // process constraints and need not care about vectorization -> go to
   // separate function
-  if (matrix_info == nullptr)
+  if (this->matrix_info == nullptr)
     {
       read_write_operation_global(operation, src);
       return;
     }
 
-  Assert(dof_info != nullptr, ExcNotInitialized());
-  Assert(matrix_info->indices_initialized() == true, ExcNotInitialized());
+  Assert(this->dof_info != nullptr, ExcNotInitialized());
+  Assert(this->matrix_info->indices_initialized() == true, ExcNotInitialized());
   if (n_fe_components == 1)
     for (unsigned int comp = 0; comp < n_components; ++comp)
       {
@@ -4025,21 +4115,22 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                           "FEEvaluation. In that case, you must pass an "
                           "std::vector<VectorType> or a BlockVector to " +
                           "read_dof_values and distribute_local_to_global."));
-        internal::check_vector_compatibility(*src[comp], *dof_info);
+        internal::check_vector_compatibility(*src[comp], *this->dof_info);
       }
   else
     {
-      internal::check_vector_compatibility(*src[0], *dof_info);
+      internal::check_vector_compatibility(*src[0], *this->dof_info);
     }
 
   // Case 2: contiguous indices which use reduced storage of indices and can
   // use vectorized load/store operations -> go to separate function
-  AssertIndexRange(cell,
-                   dof_info->index_storage_variants[dof_access_index].size());
-  if (dof_info->index_storage_variants
-        [is_face ? dof_access_index :
+  AssertIndexRange(
+    this->cell,
+    this->dof_info->index_storage_variants[this->dof_access_index].size());
+  if (this->dof_info->index_storage_variants
+        [is_face ? this->dof_access_index :
                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-        [cell] >=
+        [this->cell] >=
       internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::contiguous)
     {
       read_write_operation_contiguous(operation, src, mask);
@@ -4059,17 +4150,19 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
   const unsigned int dofs_per_component =
     this->data->dofs_per_component_on_cell;
-  if (dof_info->index_storage_variants
-        [is_face ? dof_access_index :
+  if (this->dof_info->index_storage_variants
+        [is_face ? this->dof_access_index :
                    internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
-        [cell] ==
+        [this->cell] ==
       internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::interleaved)
     {
       const unsigned int *dof_indices =
-        dof_info->dof_indices_interleaved.data() +
-        dof_info->row_starts[cell * n_fe_components * n_lanes].first +
-        dof_info->component_dof_indices_offset[active_fe_index]
-                                              [first_selected_component] *
+        this->dof_info->dof_indices_interleaved.data() +
+        this->dof_info->row_starts[this->cell * n_fe_components * n_lanes]
+          .first +
+        this->dof_info
+            ->component_dof_indices_offset[this->active_fe_index]
+                                          [this->first_selected_component] *
           n_lanes;
       if (n_components == 1 || n_fe_components == 1)
         for (unsigned int i = 0; i < dofs_per_component;
@@ -4098,27 +4191,29 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   unsigned int        cells_copied[n_lanes];
   const unsigned int *cells;
   unsigned int        n_vectorization_actual =
-    dof_info->n_vectorization_lanes_filled[dof_access_index][cell];
+    this->dof_info
+      ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
   bool has_constraints = false;
   if (is_face)
     {
-      if (dof_access_index ==
+      if (this->dof_access_index ==
           internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
         for (unsigned int v = 0; v < n_vectorization_actual; ++v)
-          cells_copied[v] = cell * VectorizedArrayType::size() + v;
-      cells = dof_access_index ==
-                  internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ?
-                &cells_copied[0] :
-                (is_interior_face ?
-                   &this->matrix_info->get_face_info(cell).cells_interior[0] :
-                   &this->matrix_info->get_face_info(cell).cells_exterior[0]);
+          cells_copied[v] = this->cell * VectorizedArrayType::size() + v;
+      cells =
+        this->dof_access_index ==
+            internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ?
+          &cells_copied[0] :
+          (this->is_interior_face ?
+             &this->matrix_info->get_face_info(this->cell).cells_interior[0] :
+             &this->matrix_info->get_face_info(this->cell).cells_exterior[0]);
       for (unsigned int v = 0; v < n_vectorization_actual; ++v)
         {
-          Assert(cells[v] < dof_info->row_starts.size() - 1,
+          Assert(cells[v] < this->dof_info->row_starts.size() - 1,
                  ExcInternalError());
           const std::pair<unsigned int, unsigned int> *my_index_start =
-            &dof_info->row_starts[cells[v] * n_fe_components +
-                                  first_selected_component];
+            &this->dof_info->row_starts[cells[v] * n_fe_components +
+                                        first_selected_component];
 
           // check whether any of the SIMD lanes has constraints, i.e., the
           // constraint indicator which is the second entry of row_starts
@@ -4127,33 +4222,34 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
             has_constraints = true;
 
           dof_indices[v] =
-            dof_info->dof_indices.data() + my_index_start[0].first;
+            this->dof_info->dof_indices.data() + my_index_start[0].first;
         }
       for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
         dof_indices[v] = nullptr;
     }
   else
     {
-      AssertIndexRange((cell + 1) * n_lanes * n_fe_components,
-                       dof_info->row_starts.size());
+      AssertIndexRange((this->cell + 1) * n_lanes * n_fe_components,
+                       this->dof_info->row_starts.size());
       const unsigned int n_components_read =
         n_fe_components > 1 ? n_components : 1;
       for (unsigned int v = 0; v < n_vectorization_actual; ++v)
         {
           const std::pair<unsigned int, unsigned int> *my_index_start =
-            &dof_info->row_starts[(cell * n_lanes + v) * n_fe_components +
-                                  first_selected_component];
+            &this->dof_info
+               ->row_starts[(this->cell * n_lanes + v) * n_fe_components +
+                            first_selected_component];
           if (my_index_start[n_components_read].second !=
               my_index_start[0].second)
             has_constraints = true;
           Assert(my_index_start[n_components_read].first ==
                      my_index_start[0].first ||
-                   my_index_start[0].first < dof_info->dof_indices.size(),
+                   my_index_start[0].first < this->dof_info->dof_indices.size(),
                  ExcIndexRange(0,
                                my_index_start[0].first,
-                               dof_info->dof_indices.size()));
+                               this->dof_info->dof_indices.size()));
           dof_indices[v] =
-            dof_info->dof_indices.data() + my_index_start[0].first;
+            this->dof_info->dof_indices.data() + my_index_start[0].first;
         }
       for (unsigned int v = n_vectorization_actual; v < n_lanes; ++v)
         dof_indices[v] = nullptr;
@@ -4200,30 +4296,33 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         operation.process_empty(values_dofs[comp][i]);
   for (unsigned int v = 0; v < n_vectorization_actual; ++v)
     {
-      const unsigned int cell_index = is_face ? cells[v] : cell * n_lanes + v;
+      const unsigned int cell_index =
+        is_face ? cells[v] : this->cell * n_lanes + v;
       const unsigned int cell_dof_index =
         cell_index * n_fe_components + first_selected_component;
       const unsigned int n_components_read =
         n_fe_components > 1 ? n_components : 1;
       unsigned int index_indicators =
-        dof_info->row_starts[cell_dof_index].second;
+        this->dof_info->row_starts[cell_dof_index].second;
       unsigned int next_index_indicators =
-        dof_info->row_starts[cell_dof_index + 1].second;
+        this->dof_info->row_starts[cell_dof_index + 1].second;
 
       // For read_dof_values_plain, redirect the dof_indices field to the
       // unconstrained indices
       if (apply_constraints == false &&
-          dof_info->row_starts[cell_dof_index].second !=
-            dof_info->row_starts[cell_dof_index + n_components_read].second)
+          this->dof_info->row_starts[cell_dof_index].second !=
+            this->dof_info->row_starts[cell_dof_index + n_components_read]
+              .second)
         {
-          Assert(dof_info->row_starts_plain_indices[cell_index] !=
+          Assert(this->dof_info->row_starts_plain_indices[cell_index] !=
                    numbers::invalid_unsigned_int,
                  ExcNotInitialized());
           dof_indices[v] =
-            dof_info->plain_dof_indices.data() +
-            dof_info->component_dof_indices_offset[active_fe_index]
-                                                  [first_selected_component] +
-            dof_info->row_starts_plain_indices[cell_index];
+            this->dof_info->plain_dof_indices.data() +
+            this->dof_info
+              ->component_dof_indices_offset[this->active_fe_index]
+                                            [this->first_selected_component] +
+            this->dof_info->row_starts_plain_indices[cell_index];
           next_index_indicators = index_indicators;
         }
 
@@ -4233,7 +4332,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
           for (; index_indicators != next_index_indicators; ++index_indicators)
             {
               const std::pair<unsigned short, unsigned short> indicator =
-                dof_info->constraint_indicator[index_indicators];
+                this->dof_info->constraint_indicator[index_indicators];
               // run through values up to next constraint
               for (unsigned int j = 0; j < indicator.first; ++j)
                 for (unsigned int comp = 0; comp < n_components; ++comp)
@@ -4252,9 +4351,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                                           value[comp]);
 
               const Number *data_val =
-                matrix_info->constraint_pool_begin(indicator.second);
+                this->matrix_info->constraint_pool_begin(indicator.second);
               const Number *end_pool =
-                matrix_info->constraint_pool_end(indicator.second);
+                this->matrix_info->constraint_pool_end(indicator.second);
               for (; data_val != end_pool; ++data_val, ++dof_indices[v])
                 for (unsigned int comp = 0; comp < n_components; ++comp)
                   operation.process_constraint(*dof_indices[v],
@@ -4292,7 +4391,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                    ++index_indicators)
                 {
                   const std::pair<unsigned short, unsigned short> indicator =
-                    dof_info->constraint_indicator[index_indicators];
+                    this->dof_info->constraint_indicator[index_indicators];
 
                   // run through values up to next constraint
                   for (unsigned int j = 0; j < indicator.first; ++j)
@@ -4309,9 +4408,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                                             value);
 
                   const Number *data_val =
-                    matrix_info->constraint_pool_begin(indicator.second);
+                    this->matrix_info->constraint_pool_begin(indicator.second);
                   const Number *end_pool =
-                    matrix_info->constraint_pool_end(indicator.second);
+                    this->matrix_info->constraint_pool_end(indicator.second);
 
                   for (; data_val != end_pool; ++data_val, ++dof_indices[v])
                     operation.process_constraint(*dof_indices[v],
@@ -4338,7 +4437,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
               if (apply_constraints == true && comp + 1 < n_components)
                 next_index_indicators =
-                  dof_info->row_starts[cell_dof_index + comp + 2].second;
+                  this->dof_info->row_starts[cell_dof_index + comp + 2].second;
             }
         }
     }
@@ -4360,15 +4459,15 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   Assert(!local_dof_indices.empty(), ExcNotInitialized());
 
   unsigned int index =
-    first_selected_component * data->dofs_per_component_on_cell;
+    first_selected_component * this->data->dofs_per_component_on_cell;
   for (unsigned int comp = 0; comp < n_components; ++comp)
     {
-      for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+      for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
            ++i, ++index)
         {
           operation.process_empty(values_dofs[comp][i]);
           operation.process_dof_global(
-            local_dof_indices[data->lexicographic_numbering[index]],
+            local_dof_indices[this->data->lexicographic_numbering[index]],
             *src[0],
             values_dofs[comp][i][0]);
         }
@@ -4401,55 +4500,57 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                          internal::is_vectorizable<VectorType, Number>::value>
                                                                vector_selector;
   const internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex ind =
-    is_face ? dof_access_index :
+    is_face ? this->dof_access_index :
               internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
   const unsigned int n_lanes = mask.count();
 
   const std::vector<unsigned int> &dof_indices_cont =
-    dof_info->dof_indices_contiguous[ind];
+    this->dof_info->dof_indices_contiguous[ind];
 
   // Simple case: We have contiguous storage, so we can simply copy out the
   // data
-  if (dof_info->index_storage_variants[ind][cell] ==
+  if (this->dof_info->index_storage_variants[ind][this->cell] ==
         internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
           interleaved_contiguous &&
       n_lanes == VectorizedArrayType::size())
     {
       const unsigned int dof_index =
-        dof_indices_cont[cell * VectorizedArrayType::size()] +
-        dof_info->component_dof_indices_offset[active_fe_index]
-                                              [first_selected_component] *
+        dof_indices_cont[this->cell * VectorizedArrayType::size()] +
+        this->dof_info->component_dof_indices_offset[this->active_fe_index]
+                                                    [first_selected_component] *
           VectorizedArrayType::size();
       if (n_components == 1 || n_fe_components == 1)
         for (unsigned int comp = 0; comp < n_components; ++comp)
-          operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
-                                            dof_index,
-                                            *src[comp],
-                                            values_dofs[comp],
-                                            vector_selector);
+          operation.process_dofs_vectorized(
+            this->data->dofs_per_component_on_cell,
+            dof_index,
+            *src[comp],
+            values_dofs[comp],
+            vector_selector);
       else
-        operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
-                                            n_components,
-                                          dof_index,
-                                          *src[0],
-                                          values_dofs[0],
-                                          vector_selector);
+        operation.process_dofs_vectorized(
+          this->data->dofs_per_component_on_cell * n_components,
+          dof_index,
+          *src[0],
+          values_dofs[0],
+          vector_selector);
       return;
     }
 
   // More general case: Must go through the components one by one and apply
   // some transformations
   const unsigned int n_filled_lanes =
-    dof_info->n_vectorization_lanes_filled[ind][this->cell];
+    this->dof_info->n_vectorization_lanes_filled[ind][this->cell];
 
   unsigned int dof_indices[VectorizedArrayType::size()];
   for (unsigned int v = 0; v < n_filled_lanes; ++v)
     dof_indices[v] =
-      dof_indices_cont[cell * VectorizedArrayType::size() + v] +
-      dof_info->component_dof_indices_offset[active_fe_index]
-                                            [first_selected_component] *
-        dof_info->dof_indices_interleave_strides
-          [ind][cell * VectorizedArrayType::size() + v];
+      dof_indices_cont[this->cell * VectorizedArrayType::size() + v] +
+      this->dof_info
+          ->component_dof_indices_offset[this->active_fe_index]
+                                        [this->first_selected_component] *
+        this->dof_info->dof_indices_interleave_strides
+          [ind][this->cell * VectorizedArrayType::size() + v];
 
   for (unsigned int v = n_filled_lanes; v < VectorizedArrayType::size(); ++v)
     dof_indices[v] = numbers::invalid_unsigned_int;
@@ -4459,32 +4560,33 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   if (n_filled_lanes == VectorizedArrayType::size() &&
       n_lanes == VectorizedArrayType::size())
     {
-      if (dof_info->index_storage_variants[ind][cell] ==
+      if (this->dof_info->index_storage_variants[ind][this->cell] ==
           internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
             contiguous)
         {
           if (n_components == 1 || n_fe_components == 1)
             for (unsigned int comp = 0; comp < n_components; ++comp)
               operation.process_dofs_vectorized_transpose(
-                data->dofs_per_component_on_cell,
+                this->data->dofs_per_component_on_cell,
                 dof_indices,
                 *src[comp],
                 values_dofs[comp],
                 vector_selector);
           else
             operation.process_dofs_vectorized_transpose(
-              data->dofs_per_component_on_cell * n_components,
+              this->data->dofs_per_component_on_cell * n_components,
               dof_indices,
               *src[0],
               &values_dofs[0][0],
               vector_selector);
         }
-      else if (dof_info->index_storage_variants[ind][cell] ==
+      else if (this->dof_info->index_storage_variants[ind][this->cell] ==
                internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                  interleaved_contiguous_strided)
         {
           if (n_components == 1 || n_fe_components == 1)
-            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+            for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
+                 ++i)
               {
                 for (unsigned int comp = 0; comp < n_components; ++comp)
                   operation.process_dof_gather(dof_indices,
@@ -4495,13 +4597,14 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               }
           else
             for (unsigned int comp = 0; comp < n_components; ++comp)
-              for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+              for (unsigned int i = 0;
+                   i < this->data->dofs_per_component_on_cell;
                    ++i)
                 {
                   operation.process_dof_gather(
                     dof_indices,
                     *src[0],
-                    (comp * data->dofs_per_component_on_cell + i) *
+                    (comp * this->data->dofs_per_component_on_cell + i) *
                       VectorizedArrayType::size(),
                     values_dofs[comp][i],
                     vector_selector);
@@ -4509,15 +4612,16 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         }
       else
         {
-          Assert(dof_info->index_storage_variants[ind][cell] ==
+          Assert(this->dof_info->index_storage_variants[ind][this->cell] ==
                    internal::MatrixFreeFunctions::DoFInfo::
                      IndexStorageVariants::interleaved_contiguous_mixed_strides,
                  ExcNotImplemented());
           const unsigned int *offsets =
-            &dof_info->dof_indices_interleave_strides
-               [ind][VectorizedArrayType::size() * cell];
+            &this->dof_info->dof_indices_interleave_strides
+               [ind][VectorizedArrayType::size() * this->cell];
           if (n_components == 1 || n_fe_components == 1)
-            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+            for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
+                 ++i)
               {
                 for (unsigned int comp = 0; comp < n_components; ++comp)
                   operation.process_dof_gather(dof_indices,
@@ -4531,7 +4635,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
               }
           else
             for (unsigned int comp = 0; comp < n_components; ++comp)
-              for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+              for (unsigned int i = 0;
+                   i < this->data->dofs_per_component_on_cell;
                    ++i)
                 {
                   operation.process_dof_gather(dof_indices,
@@ -4548,9 +4653,10 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   else
     for (unsigned int comp = 0; comp < n_components; ++comp)
       {
-        for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+        for (unsigned int i = 0; i < this->data->dofs_per_component_on_cell;
+             ++i)
           operation.process_empty(values_dofs[comp][i]);
-        if (dof_info->index_storage_variants[ind][cell] ==
+        if (this->dof_info->index_storage_variants[ind][this->cell] ==
             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
               contiguous)
           {
@@ -4559,7 +4665,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
                   if (mask[v] == true)
                     for (unsigned int i = 0;
-                         i < data->dofs_per_component_on_cell;
+                         i < this->data->dofs_per_component_on_cell;
                          ++i)
                       operation.process_dof(dof_indices[v] + i,
                                             *src[comp],
@@ -4570,11 +4676,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
                   if (mask[v] == true)
                     for (unsigned int i = 0;
-                         i < data->dofs_per_component_on_cell;
+                         i < this->data->dofs_per_component_on_cell;
                          ++i)
                       operation.process_dof(
                         dof_indices[v] + i +
-                          comp * data->dofs_per_component_on_cell,
+                          comp * this->data->dofs_per_component_on_cell,
                         *src[0],
                         values_dofs[comp][i][v]);
               }
@@ -4582,8 +4688,8 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
         else
           {
             const unsigned int *offsets =
-              &dof_info->dof_indices_interleave_strides
-                 [ind][VectorizedArrayType::size() * cell];
+              &this->dof_info->dof_indices_interleave_strides
+                 [ind][VectorizedArrayType::size() * this->cell];
             for (unsigned int v = 0; v < n_filled_lanes; ++v)
               AssertIndexRange(offsets[v], VectorizedArrayType::size() + 1);
             if (n_components == 1 || n_fe_components == 1)
@@ -4591,7 +4697,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                 {
                   if (mask[v] == true)
                     for (unsigned int i = 0;
-                         i < data->dofs_per_component_on_cell;
+                         i < this->data->dofs_per_component_on_cell;
                          ++i)
                       operation.process_dof(dof_indices[v] + i * offsets[v],
                                             *src[comp],
@@ -4602,11 +4708,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
                 for (unsigned int v = 0; v < n_filled_lanes; ++v)
                   if (mask[v] == true)
                     for (unsigned int i = 0;
-                         i < data->dofs_per_component_on_cell;
+                         i < this->data->dofs_per_component_on_cell;
                          ++i)
                       operation.process_dof(
                         dof_indices[v] +
-                          (i + comp * data->dofs_per_component_on_cell) *
+                          (i + comp * this->data->dofs_per_component_on_cell) *
                             offsets[v],
                         *src[0],
                         values_dofs[comp][i][v]);
@@ -4757,34 +4863,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 /*------------------------------ access to data fields ----------------------*/
 
-template <int dim,
-          int n_components,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline const std::vector<unsigned int> &
-FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
-  get_internal_dof_numbering() const
-{
-  return data->lexicographic_numbering;
-}
-
-
-
-template <int dim,
-          int n_components,
-          typename Number,
-          bool is_face,
-          typename VectorizedArrayType>
-inline ArrayView<VectorizedArrayType>
-FEEvaluationBase<dim, n_components, Number, is_face, VectorizedArrayType>::
-  get_scratch_data() const
-{
-  return ArrayView<VectorizedArrayType>(
-    const_cast<VectorizedArrayType *>(scratch_data),
-    scratch_data_array->end() - scratch_data);
-}
-
 
 
 template <int dim,
@@ -4981,7 +5059,7 @@ inline DEAL_II_ALWAYS_INLINE
 #  endif
 
   AssertIndexRange(q_point, this->n_quadrature_points);
-  Assert(jacobian != nullptr, ExcNotInitialized());
+  Assert(this->jacobian != nullptr, ExcNotInitialized());
   const std::size_t nqp = this->n_quadrature_points;
   Tensor<1, n_components_, Tensor<1, dim, VectorizedArrayType>> grad_out;
 
@@ -4991,15 +5069,15 @@ inline DEAL_II_ALWAYS_INLINE
       for (unsigned int d = 0; d < dim; ++d)
         for (unsigned int comp = 0; comp < n_components; comp++)
           grad_out[comp][d] = gradients_quad[(comp * dim + d) * nqp + q_point] *
-                              jacobian[0][d][d];
+                              this->jacobian[0][d][d];
     }
   // cell with general/affine Jacobian
   else
     {
       const Tensor<2, dim, VectorizedArrayType> &jac =
-        jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
-                   q_point :
-                   0];
+        this->jacobian[this->cell_type > internal::MatrixFreeFunctions::affine ?
+                         q_point :
+                         0];
       for (unsigned int comp = 0; comp < n_components; comp++)
         for (unsigned int d = 0; d < dim; ++d)
           {
@@ -5030,7 +5108,7 @@ inline DEAL_II_ALWAYS_INLINE Tensor<1, n_components_, VectorizedArrayType>
          internal::ExcAccessToUninitializedField());
 #  endif
 
-  Assert(normal_x_jacobian != nullptr, ExcNotInitialized());
+  Assert(this->normal_x_jacobian != nullptr, ExcNotInitialized());
 
   const std::size_t                            nqp = this->n_quadrature_points;
   Tensor<1, n_components, VectorizedArrayType> grad_out;
@@ -5126,11 +5204,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 #  endif
   AssertIndexRange(q_point, this->n_quadrature_points);
 
-  Assert(jacobian != nullptr, ExcNotImplemented());
+  Assert(this->jacobian != nullptr, ExcNotImplemented());
   const Tensor<2, dim, VectorizedArrayType> &jac =
-    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
-               0 :
-               q_point];
+    this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
+                     0 :
+                     q_point];
 
   Tensor<1, n_components, Tensor<2, dim, VectorizedArrayType>> hessian_out;
 
@@ -5205,7 +5283,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   else
     {
       const auto &jac_grad =
-        mapping_data->jacobian_gradients
+        this->mapping_data->jacobian_gradients
           [1 - this->is_interior_face]
           [this->mapping_data->data_index_offsets[this->cell] + q_point];
       for (unsigned int comp = 0; comp < n_components; comp++)
@@ -5267,11 +5345,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 #  endif
   AssertIndexRange(q_point, this->n_quadrature_points);
 
-  Assert(jacobian != nullptr, ExcNotImplemented());
+  Assert(this->jacobian != nullptr, ExcNotImplemented());
   const Tensor<2, dim, VectorizedArrayType> &jac =
-    jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
-               0 :
-               q_point];
+    this->jacobian[this->cell_type <= internal::MatrixFreeFunctions::affine ?
+                     0 :
+                     q_point];
 
   const std::size_t      nqp  = this->n_quadrature_points;
   constexpr unsigned int hdim = (dim * (dim + 1)) / 2;
@@ -5312,7 +5390,7 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     {
       const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, VectorizedArrayType>>
         &jac_grad =
-          mapping_data->jacobian_gradients
+          this->mapping_data->jacobian_gradients
             [0][this->mapping_data->data_index_offsets[this->cell] + q_point];
       for (unsigned int comp = 0; comp < n_components; comp++)
         {
@@ -5412,13 +5490,14 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   const std::size_t nqp = this->n_quadrature_points;
   if (this->cell_type <= internal::MatrixFreeFunctions::affine)
     {
-      const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
+      const VectorizedArrayType JxW =
+        this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int comp = 0; comp < n_components; ++comp)
         values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
     }
   else
     {
-      const VectorizedArrayType JxW = J_value[q_point];
+      const VectorizedArrayType JxW = this->J_value[q_point];
       for (unsigned int comp = 0; comp < n_components; ++comp)
         values_quad[comp * nqp + q_point] = val_in[comp] * JxW;
     }
@@ -5448,10 +5527,11 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
   const std::size_t nqp = this->n_quadrature_points;
   if (!is_face && this->cell_type == internal::MatrixFreeFunctions::cartesian)
     {
-      const VectorizedArrayType JxW = J_value[0] * quadrature_weights[q_point];
+      const VectorizedArrayType JxW =
+        this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int d = 0; d < dim; ++d)
         {
-          const VectorizedArrayType factor = jacobian[0][d][d] * JxW;
+          const VectorizedArrayType factor = this->jacobian[0][d][d] * JxW;
           for (unsigned int comp = 0; comp < n_components; comp++)
             gradients_quad[(comp * dim + d) * nqp + q_point] =
               grad_in[comp][d] * factor;
@@ -5461,12 +5541,12 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
     {
       const Tensor<2, dim, VectorizedArrayType> jac =
         this->cell_type > internal::MatrixFreeFunctions::affine ?
-          jacobian[q_point] :
-          jacobian[0];
+          this->jacobian[q_point] :
+          this->jacobian[0];
       const VectorizedArrayType JxW =
         this->cell_type > internal::MatrixFreeFunctions::affine ?
-          J_value[q_point] :
-          J_value[0] * quadrature_weights[q_point];
+          this->J_value[q_point] :
+          this->J_value[0] * this->quadrature_weights[q_point];
       for (unsigned int comp = 0; comp < n_components; ++comp)
         for (unsigned int d = 0; d < dim; ++d)
           {
@@ -5591,22 +5671,19 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-template <int n_components_other>
 inline FEEvaluationAccess<dim,
                           n_components_,
                           Number,
                           is_face,
                           VectorizedArrayType>::
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other)
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other)
   : FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>(
       mapping,
       fe,
@@ -5692,18 +5769,15 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
 
 
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
 inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other)
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other)
   : FEEvaluationBase<dim, 1, Number, is_face, VectorizedArrayType>(
       mapping,
       fe,
@@ -5994,18 +6068,15 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
 
 
 template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
 inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
-  FEEvaluationAccess(const Mapping<dim> &      mapping,
-                     const FiniteElement<dim> &fe,
-                     const Quadrature<1> &     quadrature,
-                     const UpdateFlags         update_flags,
-                     const unsigned int        first_selected_component,
-                     const FEEvaluationBase<dim,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other)
+  FEEvaluationAccess(
+    const Mapping<dim> &      mapping,
+    const FiniteElement<dim> &fe,
+    const Quadrature<1> &     quadrature,
+    const UpdateFlags         update_flags,
+    const unsigned int        first_selected_component,
+    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
+      *other)
   : FEEvaluationBase<dim, dim, Number, is_face, VectorizedArrayType>(
       mapping,
       fe,
@@ -6392,18 +6463,14 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
 
 
 template <typename Number, bool is_face, typename VectorizedArrayType>
-template <int n_components_other>
 inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
-  FEEvaluationAccess(const Mapping<1> &      mapping,
-                     const FiniteElement<1> &fe,
-                     const Quadrature<1> &   quadrature,
-                     const UpdateFlags       update_flags,
-                     const unsigned int      first_selected_component,
-                     const FEEvaluationBase<1,
-                                            n_components_other,
-                                            Number,
-                                            is_face,
-                                            VectorizedArrayType> *other)
+  FEEvaluationAccess(
+    const Mapping<1> &      mapping,
+    const FiniteElement<1> &fe,
+    const Quadrature<1> &   quadrature,
+    const UpdateFlags       update_flags,
+    const unsigned int      first_selected_component,
+    const FEEvaluationBaseData<1, Number, is_face, VectorizedArrayType> *other)
   : FEEvaluationBase<1, 1, Number, is_face, VectorizedArrayType>(
       mapping,
       fe,
@@ -6708,9 +6775,7 @@ inline FEEvaluation<dim,
               quadrature,
               update_flags,
               first_selected_component,
-              static_cast<
-                FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
-                nullptr))
+              nullptr)
   , dofs_per_component(this->data->dofs_per_component_on_cell)
   , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
   , n_q_points(this->data->n_q_points)
@@ -6741,9 +6806,7 @@ inline FEEvaluation<dim,
               quadrature,
               update_flags,
               first_selected_component,
-              static_cast<
-                FEEvaluationBase<dim, 1, Number, false, VectorizedArrayType> *>(
-                nullptr))
+              nullptr)
   , dofs_per_component(this->data->dofs_per_component_on_cell)
   , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
   , n_q_points(this->data->n_q_points)
@@ -6759,20 +6822,16 @@ template <int dim,
           int n_components_,
           typename Number,
           typename VectorizedArrayType>
-template <int n_components_other>
 inline FEEvaluation<dim,
                     fe_degree,
                     n_q_points_1d,
                     n_components_,
                     Number,
                     VectorizedArrayType>::
-  FEEvaluation(const FiniteElement<dim> &                   fe,
-               const FEEvaluationBase<dim,
-                                      n_components_other,
-                                      Number,
-                                      false,
-                                      VectorizedArrayType> &other,
-               const unsigned int first_selected_component)
+  FEEvaluation(
+    const FiniteElement<dim> &                                           fe,
+    const FEEvaluationBaseData<dim, Number, false, VectorizedArrayType> &other,
+    const unsigned int first_selected_component)
   : BaseClass(other.mapped_geometry->get_fe_values().get_mapping(),
               fe,
               other.mapped_geometry->get_quadrature(),

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.