]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make evaluation_kernels only include FEEvaluationData file
authorMartin Kronbichler <martin.kronbichler@it.uu.se>
Tue, 30 Nov 2021 21:37:52 +0000 (22:37 +0100)
committerMartin Kronbichler <martin.kronbichler@it.uu.se>
Mon, 6 Dec 2021 13:28:08 +0000 (14:28 +0100)
include/deal.II/base/array_view.h
include/deal.II/matrix_free/dof_info.h
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/evaluation_template_factory.templates.h
include/deal.II/matrix_free/evaluation_template_factory_hanging_nodes.templates.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/fe_evaluation_base_data.h [new file with mode: 0644]
include/deal.II/matrix_free/mapping_info.h
include/deal.II/matrix_free/mapping_info.templates.h
include/deal.II/matrix_free/mapping_info_storage.h [new file with mode: 0644]
include/deal.II/matrix_free/mapping_info_storage.templates.h [new file with mode: 0644]

index f38a7d14ae3e880dd493b03c5a49e79d029e9afd..6a4e2017dd2235d89f75e0f14aa87fec8a28f017 100644 (file)
@@ -414,7 +414,8 @@ inline void
 ArrayView<ElementType, MemorySpaceType>::reinit(value_type *starting_element,
                                                 const std::size_t n_elements)
 {
-  *this = ArrayView(starting_element, n_elements);
+  this->starting_element = starting_element;
+  this->n_elements       = n_elements;
 }
 
 
index 08178f7eaba6d51b2625c6d09bd8ed12d4116d48..7b0583ea9e9af29397fc1b16ca9f06dae4cac215 100644 (file)
@@ -31,7 +31,6 @@
 
 #include <deal.II/matrix_free/face_info.h>
 #include <deal.II/matrix_free/hanging_nodes_internal.h>
-#include <deal.II/matrix_free/mapping_info.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/task_info.h>
 #include <deal.II/matrix_free/vector_data_exchange.h>
@@ -49,7 +48,10 @@ namespace internal
   {
     template <int dim>
     class HangingNodes;
-  }
+
+    template <typename, typename>
+    struct FPArrayComparator;
+  } // namespace MatrixFreeFunctions
 } // namespace internal
 #endif
 
@@ -86,7 +88,7 @@ namespace internal
       std::pair<std::vector<Number>, types::global_dof_index> next_constraint;
       std::map<std::vector<Number>,
                types::global_dof_index,
-               FPArrayComparator<Number>>
+               FPArrayComparator<Number, VectorizedArray<Number>>>
         constraints;
     };
 
index 0532e46907e03ada13ace0c969fe6cfd0048c59c..dc84931d6343e682a11b4d8cc2d6447d99743677 100644 (file)
@@ -26,6 +26,7 @@
 
 #include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/evaluation_flags.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
 #include <deal.II/matrix_free/hanging_nodes_internal.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 DEAL_II_NAMESPACE_OPEN
 
 
-// forward declaration
-template <int, typename, bool, typename>
-class FEEvaluationBaseData;
-
-
-
 namespace internal
 {
   // Select evaluator type from element shape function type
index 5a98c3781b270a5dafd5ae1b289c881d81585589..3f6272b615d21e1f40eb55065c5d4354370c6850 100644 (file)
@@ -24,7 +24,6 @@
 #include <deal.II/matrix_free/evaluation_selector.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/evaluation_template_factory_internal.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index cf03c7c2d5987062e324f267b235e346b907475b..3c37bdd3391eeb2cff3d4c14da04ed5e0d0ac029 100644 (file)
@@ -23,7 +23,7 @@
 #include <deal.II/matrix_free/evaluation_kernels_hanging_nodes.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/evaluation_template_factory_internal.h>
-#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
 
 DEAL_II_NAMESPACE_OPEN
 
index 4fc5c6cc2a7d3d81d20e3963a1b02798dbcf9b09..d74ef0310780c3a78b9c4e8f5be7b6e1bbe5c832 100644 (file)
@@ -33,6 +33,7 @@
 #include <deal.II/matrix_free/evaluation_kernels.h>
 #include <deal.II/matrix_free/evaluation_selector.h>
 #include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/fe_evaluation_base_data.h>
 #include <deal.II/matrix_free/mapping_data_on_the_fly.h>
 #include <deal.II/matrix_free/matrix_free.h>
 #include <deal.II/matrix_free/shape_info.h>
@@ -45,461 +46,6 @@ DEAL_II_NAMESPACE_OPEN
 
 
 
-namespace internal
-{
-  DeclException0(ExcAccessToUninitializedField);
-
-  DeclException1(
-    ExcMatrixFreeAccessToUninitializedMappingField,
-    std::string,
-    << "You are requesting information from an FEEvaluation/FEFaceEvaluation "
-    << "object for which this kind of information has not been computed. What "
-    << "information these objects compute is determined by the update_* flags you "
-    << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, "
-    << "the operation you are attempting requires the <" << arg1
-    << "> flag to be set, but it was apparently not specified "
-    << "upon initialization.");
-} // namespace internal
-
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d            = fe_degree + 1,
-          int n_components_            = 1,
-          typename Number              = double,
-          typename VectorizedArrayType = VectorizedArray<Number>>
-class FEEvaluation;
-
-
-
-/**
- * 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 four template arguments:
- *
- * @tparam dim Dimension in which this class is to be used
- *
- * @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, typename Number, bool is_face, typename VectorizedArrayType>
-class FEEvaluationBaseData
-{
-  static_assert(
-    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
-    "Type of Number and of VectorizedArrayType do not match.");
-
-public:
-  static constexpr unsigned int dimension = dim;
-
-  /**
-   * Destructor.
-   */
-  ~FEEvaluationBaseData();
-
-  /**
-   * Return the index offset within the geometry fields for the cell the @p
-   * reinit() function has been called for. This index can be used to access
-   * an index into a field that has the same compression behavior as the
-   * Jacobian of the geometry, e.g., to store an effective coefficient tensors
-   * that combines a coefficient with the geometry for lower memory transfer
-   * as the available data fields.
-   */
-  unsigned int
-  get_mapping_data_index_offset() const;
-
-  /**
-   * Return the type of the cell the @p reinit() function has been called for.
-   * Valid values are @p cartesian for Cartesian cells (which allows for
-   * considerable data compression), @p affine for cells with affine mappings,
-   * and @p general for general cells without any compressed storage applied.
-   */
-  internal::MatrixFreeFunctions::GeometryType
-  get_cell_type() const;
-
-  /**
-   * Return a reference to the ShapeInfo object currently in use.
-   */
-  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
-  get_shape_info() const;
-
-  /**
-   * Return a reference to the DoFInfo object currently in use.
-   */
-  const internal::MatrixFreeFunctions::DoFInfo &
-  get_dof_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 $J^{-\mathrm T}$ of the
-   * Jacobian of the mapping between the unit to the real cell defined as
-   * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor
-   * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space
-   * coordinates and rows to real cell coordinates. Thus, the returned tensor
-   * represents a covariant transformation, which is used in the
-   * FEEvaluationBase::get_gradient() function to transform the unit cell
-   * gradients to gradients on the real cell by a multiplication $J^{-\mathrm
-   * T} \hat{\nabla} u_h$.
-   */
-  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_cell_batches() +
-   * MatrixFree::n_ghost_cell_batches() 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_cell_batches() +
-   * MatrixFree::n_ghost_cell_batches() 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;
-
-  /**
-   * Return the number of the quadrature formula of the present cell.
-   */
-  unsigned int
-  get_quadrature_index() const;
-
-  /**
-   * Return index of the current cell or face.
-   */
-  unsigned int
-  get_current_cell_index() const;
-
-  /**
-   * Return the active FE index for this class for efficient indexing in the
-   * hp-case.
-   */
-  unsigned int
-  get_active_fe_index() const;
-
-  /**
-   * Return the active quadrature index for this class for efficient indexing in
-   * the hp-case.
-   */
-  unsigned int
-  get_active_quadrature_index() const;
-
-  /**
-   * Return the underlying MatrixFree object.
-   */
-  const MatrixFree<dim, Number, VectorizedArrayType> &
-  get_matrix_free() 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,
-    const unsigned int active_fe_index,
-    const unsigned int active_quad_index,
-    const unsigned int face_type);
-
-  /**
-   * 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;
-
-  /**
-   * 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;
-
-  /**
-   * 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;
-
-  /**
-   * A pointer to the underlying quadrature formula specified at construction.
-   * Also contained in matrix_info, but it simplifies code if we store a
-   * reference to it.
-   */
-  const typename internal::MatrixFreeFunctions::MappingInfoStorage<
-    (is_face ? dim - 1 : dim),
-    dim,
-    Number,
-    VectorizedArrayType>::QuadratureDescriptor *descriptor;
-
-  /**
-   * The number of quadrature points in the current evaluation context.
-   */
-  const unsigned int n_quadrature_points;
-
-  /**
-   * 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. For cells it defines if
-   * the dof values should be read from the actual cell corresponding to the
-   * interior face or the neighboring cell corresponding to the exterior face.
-   */
-  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
@@ -540,8 +86,8 @@ protected:
 template <int dim,
           int n_components_,
           typename Number,
-          bool is_face                 = false,
-          typename VectorizedArrayType = VectorizedArray<Number>>
+          bool is_face,
+          typename VectorizedArrayType>
 class FEEvaluationBase
   : public FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>
 {
@@ -1143,12 +689,70 @@ public:
 
   //@}
 
+  /**
+   * Provides a unified interface to access data in a vector of
+   * VectorizedArray fields of length MatrixFree::n_cell_batches() +
+   * MatrixFree::n_ghost_cell_batches() 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_cell_batches() +
+   * MatrixFree::n_ghost_cell_batches() 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 first selected component.
    */
   unsigned int
   get_first_selected_component() const;
 
+  /**
+   * Return the underlying MatrixFree object.
+   */
+  const MatrixFree<dim, Number, VectorizedArrayType> &
+  get_matrix_free() const;
+
 protected:
   /**
    * Constructor. Made protected to prevent users from directly using this
@@ -1232,6 +836,11 @@ protected:
   FEEvaluationBase &
   operator=(const FEEvaluationBase &other);
 
+  /**
+   * Destructor.
+   */
+  ~FEEvaluationBase();
+
   /**
    * A unified function to read from and write into vectors based on the given
    * template operation. It can perform the operation for @p read_dof_values,
@@ -1286,6 +895,16 @@ protected:
   void
   apply_hanging_node_constraints() const;
 
+  /**
+   * This is the general array for all data fields.
+   */
+  AlignedVector<VectorizedArrayType> *scratch_data_array;
+
+  /**
+   * A pointer to the underlying data.
+   */
+  const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info;
+
   /**
    * This field stores the values for local degrees of freedom (e.g. after
    * reading out from a vector but before applying unit cell transformations
@@ -3470,83 +3089,174 @@ namespace internal
 #ifndef DOXYGEN
 
 
-/*----------------------- FEEvaluationBaseData ------------------------*/
+namespace internal
+{
+  template <typename MatrixFreeType>
+  inline unsigned int
+  get_active_fe_index(const MatrixFreeType &data,
+                      const unsigned int    dof_no,
+                      const unsigned int    fe_degree,
+                      const unsigned int    first_selected_component,
+                      const unsigned int    active_fe_index)
+  {
+    return fe_degree != numbers::invalid_unsigned_int ?
+             data.get_dof_info(dof_no).fe_index_from_degree(
+               first_selected_component, fe_degree) :
+             (active_fe_index != numbers::invalid_unsigned_int ?
+                active_fe_index :
+                0);
+  }
 
-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,
-    const unsigned int active_fe_index_in,
-    const unsigned int active_quad_index_in,
-    const unsigned int face_type)
-  : scratch_data_array(data_in.acquire_scratch_data())
-  , quad_no(quad_no_in)
-  , matrix_info(&data_in)
-  , dof_info(&data_in.get_dof_info(dof_no))
-  , mapping_data(
+
+
+  template <bool is_face,
+            int  dim,
+            typename Number,
+            typename VectorizedArrayType>
+  inline unsigned int
+  get_active_quad_index(const MatrixFree<dim, Number, VectorizedArrayType> &mf,
+                        const unsigned int quad_no,
+                        const unsigned int fe_degree,
+                        const unsigned int n_q_points,
+                        const unsigned int active_fe_index,
+                        const unsigned int active_quad_index)
+  {
+    const auto &mapping_storage = internal::MatrixFreeFunctions::
+      MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
+        mf.get_mapping_info(), quad_no);
+    if (fe_degree == numbers::invalid_unsigned_int)
+      return active_quad_index != numbers::invalid_unsigned_int ?
+               active_quad_index :
+               std::min<unsigned int>(active_fe_index,
+                                      mapping_storage.descriptor.size() - 1);
+    else
+      return mapping_storage.quad_index_from_n_q_points(n_q_points);
+  }
+} // namespace internal
+
+
+/*----------------------- 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,
+                   const unsigned int fe_degree,
+                   const unsigned int n_q_points,
+                   const bool         is_interior_face,
+                   const unsigned int active_fe_index,
+                   const unsigned int active_quad_index,
+                   const unsigned int face_type)
+  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+      data_in.get_shape_info(
+        dof_no,
+        quad_no,
+        data_in.get_dof_info(dof_no)
+          .component_to_base_index[first_selected_component],
+        internal::get_active_fe_index(data_in,
+                                      dof_no,
+                                      fe_degree,
+                                      first_selected_component,
+                                      active_fe_index),
+        internal::get_active_quad_index<is_face>(
+          data_in,
+          quad_no,
+          fe_degree,
+          n_q_points,
+          internal::get_active_fe_index(data_in,
+                                        dof_no,
+                                        fe_degree,
+                                        first_selected_component,
+                                        active_fe_index),
+          active_quad_index)),
+      data_in.get_dof_info(dof_no),
       internal::MatrixFreeFunctions::
         MappingInfoCellsOrFaces<dim, Number, is_face, VectorizedArrayType>::get(
           data_in.get_mapping_info(),
-          quad_no))
-  , active_fe_index(fe_degree != numbers::invalid_unsigned_int ?
-                      data_in.get_dof_info(dof_no).fe_index_from_degree(
-                        first_selected_component,
-                        fe_degree) :
-                      (active_fe_index_in != numbers::invalid_unsigned_int ?
-                         active_fe_index_in :
-                         0))
-  , active_quad_index(
-      fe_degree != numbers::invalid_unsigned_int ?
-        (mapping_data->quad_index_from_n_q_points(n_q_points)) :
-        (active_quad_index_in != numbers::invalid_unsigned_int ?
-           active_quad_index_in :
-           std::min<unsigned int>(active_fe_index,
-                                  mapping_data->descriptor.size() - 1)))
-  , descriptor(
-      &mapping_data->descriptor
-         [is_face ?
-            (active_quad_index * std::max<unsigned int>(1, dim - 1) +
-             (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
-            active_quad_index])
-  , n_quadrature_points(descriptor->n_q_points)
-  , data(&data_in.get_shape_info(
-      dof_no,
-      quad_no_in,
-      dof_info->component_to_base_index[first_selected_component],
-      active_fe_index,
-      active_quad_index))
-  , jacobian(nullptr)
-  , J_value(nullptr)
-  , normal_vectors(nullptr)
-  , normal_x_jacobian(nullptr)
-  , quadrature_weights(descriptor->quadrature_weights.begin())
-  , cell(numbers::invalid_unsigned_int)
-  , is_interior_face(is_interior_face)
-  , dof_access_index(
-      is_face ?
-        (is_interior_face ?
-           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
-           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
-        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
-  , cell_type(internal::MatrixFreeFunctions::general)
-{
-  Assert(matrix_info->mapping_initialized() == true, ExcNotInitialized());
-  AssertDimension(matrix_info->get_task_info().vectorization_length,
-                  VectorizedArrayType::size());
-  AssertDimension(n_quadrature_points, descriptor->n_q_points);
+          quad_no),
+      quad_no,
+      is_interior_face,
+      internal::get_active_fe_index(data_in,
+                                    dof_no,
+                                    fe_degree,
+                                    first_selected_component,
+                                    active_fe_index),
+      internal::get_active_quad_index<is_face>(
+        data_in,
+        quad_no,
+        fe_degree,
+        n_q_points,
+        internal::get_active_fe_index(data_in,
+                                      dof_no,
+                                      fe_degree,
+                                      first_selected_component,
+                                      active_fe_index),
+        active_quad_index),
+      face_type)
+  , scratch_data_array(data_in.acquire_scratch_data())
+  , matrix_info(&data_in)
+  , 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) + ")"));
+
+  // do not check for correct dimensions of data fields here, should be done
+  // in derived classes
 }
 
 
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  FEEvaluationBaseData(
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType>::
+  FEEvaluationBase(
     const Mapping<dim> &      mapping,
     const FiniteElement<dim> &fe,
     const Quadrature<1> &     quadrature,
@@ -3554,115 +3264,129 @@ inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
     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)
-  , active_fe_index(numbers::invalid_unsigned_int)
-  , active_quad_index(numbers::invalid_unsigned_int)
-  , descriptor(nullptr)
-  , n_quadrature_points(
-      Utilities::fixed_power < is_face ? dim - 1 : dim > (quadrature.size()))
+  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
+      other != nullptr &&
+          other->mapped_geometry->get_quadrature() == quadrature ?
+        other->mapped_geometry :
+        std::make_shared<
+          internal::MatrixFreeFunctions::
+            MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
+          mapping,
+          quadrature,
+          update_flags))
+  , scratch_data_array(new AlignedVector<VectorizedArrayType>())
   , matrix_info(nullptr)
-  , dof_info(nullptr)
-  , mapping_data(nullptr)
-  ,
-  // select the correct base element from the given FE component
-  data(new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
-    Quadrature<dim - is_face>(quadrature),
-    fe,
-    fe.component_to_base_index(first_selected_component).first))
-  , jacobian(nullptr)
-  , 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)
-{
-  Assert(other == nullptr || other->mapped_geometry.get() != nullptr,
-         ExcInternalError());
-  if (other != nullptr &&
-      other->mapped_geometry->get_quadrature() == quadrature)
-    mapped_geometry = other->mapped_geometry;
-  else
-    mapped_geometry =
-      std::make_shared<internal::MatrixFreeFunctions::
-                         MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
-        mapping, quadrature, update_flags);
-  cell = 0;
+  , 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)
+{
+  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;
+
+  Assert(this->data == nullptr, ExcInternalError());
+  this->data =
+    new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+      Quadrature<dim - is_face>(quadrature),
+      fe,
+      fe.component_to_base_index(first_selected_component).first);
 
-  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();
+  set_data_pointers();
 }
 
 
 
-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 ?
+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)
+  : 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)
-  , active_fe_index(other.active_fe_index)
-  , active_quad_index(other.active_quad_index)
-  , descriptor(other.descriptor == nullptr ? nullptr : other.descriptor)
-  , n_quadrature_points(other.n_quadrature_points)
   , matrix_info(other.matrix_info)
-  , dof_info(other.dof_info)
-  , mapping_data(other.mapping_data)
-  , data(other.matrix_info == nullptr ?
-           new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
-             *other.data) :
-           other.data)
-  , jacobian(nullptr)
-  , J_value(nullptr)
-  , normal_vectors(nullptr)
-  , normal_x_jacobian(nullptr)
-  , quadrature_weights(other.descriptor == nullptr ?
-                         nullptr :
-                         descriptor->quadrature_weights.begin())
-  , cell(numbers::invalid_unsigned_int)
-  , cell_type(internal::MatrixFreeFunctions::general)
-  , is_interior_face(other.is_interior_face)
-  , dof_access_index(other.dof_access_index)
-{
-  // Create deep copy of mapped geometry for use in parallel...
-  if (other.mapped_geometry.get() != nullptr)
+  , 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)
+{
+  if (other.matrix_info == nullptr)
     {
-      mapped_geometry = std::make_shared<
+      Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
+      this->data =
+        new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+          *other.data);
+
+      // Create deep copy of mapped geometry for use in parallel
+      this->mapped_geometry = std::make_shared<
         internal::MatrixFreeFunctions::
           MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
         other.mapped_geometry->get_fe_values().get_mapping(),
         other.mapped_geometry->get_quadrature(),
         other.mapped_geometry->get_fe_values().get_update_flags());
-      mapping_data = &mapped_geometry->get_data_storage();
-      cell         = 0;
+      this->mapping_data = &this->mapped_geometry->get_data_storage();
+      this->cell         = 0;
 
-      jacobian = mapped_geometry->get_data_storage().jacobians[0].begin();
-      J_value  = mapped_geometry->get_data_storage().JxW_values.begin();
+      this->jacobian =
+        this->mapped_geometry->get_data_storage().jacobians[0].begin();
+      this->J_value =
+        this->mapped_geometry->get_data_storage().JxW_values.begin();
     }
+
+  set_data_pointers();
 }
 
 
 
-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)
+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)
 {
-  AssertDimension(quad_no, other.quad_no);
-  AssertDimension(active_fe_index, other.active_fe_index);
-  AssertDimension(active_quad_index, other.active_quad_index);
-
   // release old memory
   if (matrix_info == nullptr)
     {
-      delete data;
+      delete this->data;
       delete scratch_data_array;
     }
   else
@@ -3670,182 +3394,156 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::operator=(
       matrix_info->release_scratch_data(scratch_data_array);
     }
 
-  matrix_info  = other.matrix_info;
-  dof_info     = other.dof_info;
-  descriptor   = other.descriptor;
-  mapping_data = other.mapping_data;
+  this->FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  operator=(other);
+
+  matrix_info = other.matrix_info;
+
   if (other.matrix_info == nullptr)
     {
-      data = new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
-        *other.data);
+      Assert(other.mapped_geometry.get() != nullptr, ExcInternalError());
+      this->data =
+        new internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>(
+          *other.data);
       scratch_data_array = new AlignedVector<VectorizedArrayType>();
-    }
-  else
-    {
-      data               = other.data;
-      scratch_data_array = matrix_info->acquire_scratch_data();
-    }
 
-  quadrature_weights =
-    (descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr);
-  cell             = numbers::invalid_unsigned_int;
-  cell_type        = internal::MatrixFreeFunctions::general;
-  is_interior_face = other.is_interior_face;
-  dof_access_index = other.dof_access_index;
-
-  // Create deep copy of mapped geometry for use in parallel...
-  if (other.mapped_geometry.get() != nullptr)
-    {
-      mapped_geometry = std::make_shared<
+      // Create deep copy of mapped geometry for use in parallel
+      this->mapped_geometry = std::make_shared<
         internal::MatrixFreeFunctions::
           MappingDataOnTheFly<dim, Number, VectorizedArrayType>>(
         other.mapped_geometry->get_fe_values().get_mapping(),
         other.mapped_geometry->get_quadrature(),
         other.mapped_geometry->get_fe_values().get_update_flags());
-      cell         = 0;
-      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();
-    }
-
-  return *this;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  ~FEEvaluationBaseData()
-{
-  if (matrix_info != nullptr)
-    {
-      try
-        {
-          matrix_info->release_scratch_data(scratch_data_array);
-        }
-      catch (...)
-        {}
-    }
-  else
-    {
-      delete scratch_data_array;
-      delete data;
-      data = nullptr;
+      this->cell         = 0;
+      this->mapping_data = &this->mapped_geometry->get_data_storage();
+      this->jacobian =
+        this->mapped_geometry->get_data_storage().jacobians[0].begin();
+      this->J_value =
+        this->mapped_geometry->get_data_storage().JxW_values.begin();
     }
-  scratch_data_array = nullptr;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_mapping_data_index_offset() const
-{
-  if (matrix_info == nullptr)
-    return 0;
   else
     {
-      AssertIndexRange(cell, this->mapping_data->data_index_offsets.size());
-      return this->mapping_data->data_index_offsets[cell];
+      scratch_data_array = matrix_info->acquire_scratch_data();
     }
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline internal::MatrixFreeFunctions::GeometryType
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
-  const
-{
-  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
-  return cell_type;
-}
 
+  AssertDimension(n_fe_components, other.n_fe_components);
+  AssertDimension(first_selected_component, other.first_selected_component);
 
+  set_data_pointers();
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_shape_info() const
-{
-  Assert(data != nullptr, ExcInternalError());
-  return *data;
+  return *this;
 }
 
 
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const internal::MatrixFreeFunctions::DoFInfo &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_dof_info()
-  const
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+inline FEEvaluationBase<dim,
+                        n_components_,
+                        Number,
+                        is_face,
+                        VectorizedArrayType>::~FEEvaluationBase()
 {
-  Assert(dof_info != nullptr,
-         ExcMessage(
-           "FEEvaluation was not initialized with a MatrixFree object!"));
-  return *dof_info;
+  if (matrix_info != nullptr)
+    {
+      try
+        {
+          matrix_info->release_scratch_data(scratch_data_array);
+        }
+      catch (...)
+        {}
+    }
+  else
+    {
+      delete scratch_data_array;
+      delete this->data;
+    }
 }
 
 
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_normal_vector(const unsigned int q_point) const
+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()
 {
-  AssertIndexRange(q_point, n_quadrature_points);
-  Assert(normal_vectors != nullptr,
-         internal::ExcMatrixFreeAccessToUninitializedMappingField(
-           "update_normal_vectors"));
-  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
-    return normal_vectors[0];
-  else
-    return normal_vectors[q_point];
-}
+  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 = this->n_quadrature_points;
 
+  const unsigned int size_scratch_data =
+    std::max(tensor_dofs_per_component + 1, dofs_per_component) *
+      n_components_ * 3 +
+    2 * n_quadrature_points;
+  const unsigned int allocated_size =
+    size_scratch_data + n_components_ * dofs_per_component +
+    (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+     n_quadrature_points);
+  this->scratch_data_array->resize_fast(allocated_size);
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
-  const unsigned int q_point) const
-{
-  AssertIndexRange(q_point, n_quadrature_points);
-  Assert(J_value != nullptr,
-         internal::ExcMatrixFreeAccessToUninitializedMappingField(
-           "update_values|update_gradients"));
-  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+  // set the pointers to the correct position in the data array
+  for (unsigned int c = 0; c < n_components_; ++c)
     {
-      Assert(this->quadrature_weights != nullptr, ExcInternalError());
-      return J_value[0] * this->quadrature_weights[q_point];
+      values_dofs[c] =
+        this->scratch_data_array->begin() + c * dofs_per_component;
     }
-  else
-    return J_value[q_point];
+  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);
+  gradients_from_hessians_quad =
+    this->scratch_data_array->begin() +
+    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
+  hessians_quad =
+    this->scratch_data_array->begin() +
+    n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
+
+  this->scratch_data.reinit(this->scratch_data_array->begin() +
+                              n_components_ * dofs_per_component +
+                              (n_components_ *
+                               ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
+                               n_quadrature_points),
+                            size_scratch_data);
 }
 
 
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  inverse_jacobian(const unsigned int q_point) const
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+inline const MatrixFree<dim, Number, VectorizedArrayType> &
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  get_matrix_free() const
 {
-  AssertIndexRange(q_point, n_quadrature_points);
-  Assert(this->jacobian != nullptr,
-         internal::ExcMatrixFreeAccessToUninitializedMappingField(
-           "update_gradients"));
-  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
-    return jacobian[0];
-  else
-    return jacobian[q_point];
+  Assert(matrix_info != nullptr,
+         ExcMessage(
+           "FEEvaluation was not initialized with a MatrixFree object!"));
+  return *matrix_info;
 }
 
 
 
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
 inline std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
-  const
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  get_cell_ids() const
 {
   Assert(this->matrix_info != nullptr, ExcNotInitialized());
 
@@ -3864,7 +3562,7 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
     {
       // cell or interior face face (element-centric loop)
       for (unsigned int i = 0; i < n_lanes; ++i)
-        cells[i] = cell * n_lanes + i;
+        cells[i] = this->cell * n_lanes + i;
     }
   else if (is_face &&
            this->dof_access_index ==
@@ -3905,9 +3603,9 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
     {
       // face-centric faces
       const unsigned int *cells_ =
-        is_interior_face ?
-          &this->matrix_info->get_face_info(cell).cells_interior[0] :
-          &this->matrix_info->get_face_info(cell).cells_exterior[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 i = 0; i < VectorizedArrayType::size(); ++i)
         if (cells_[i] != numbers::invalid_unsigned_int)
           cells[i] = cells_[i];
@@ -3920,6 +3618,7 @@ FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_ids()
 namespace internal
 {
   template <int dim,
+            int n_components,
             typename Number,
             bool is_face,
             typename VectorizedArrayType,
@@ -3928,7 +3627,11 @@ namespace internal
             typename FU>
   inline void
   process_cell_data(
-    const FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType> &phi,
+    const FEEvaluationBase<dim,
+                           n_components,
+                           Number,
+                           is_face,
+                           VectorizedArrayType> &       phi,
     const MatrixFree<dim, Number, VectorizedArrayType> *matrix_info,
     GlobalVectorType &                                  array,
     VectorizedArrayType2 &                              out,
@@ -3946,277 +3649,69 @@ namespace internal
     for (unsigned int i = 0; i < VectorizedArrayType::size(); ++i)
       if (cells[i] != numbers::invalid_unsigned_int)
         fu(out[i],
-           array[cells[i] / VectorizedArrayType::size()]
-                [cells[i] % VectorizedArrayType::size()]);
-  }
-} // namespace internal
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-std::array<unsigned int, VectorizedArrayType::size()>
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_cell_or_face_ids() const
-{
-  const unsigned int v_len = VectorizedArrayType::size();
-  std::array<unsigned int, VectorizedArrayType::size()> cells;
-
-  // initialize array
-  for (unsigned int i = 0; i < v_len; ++i)
-    cells[i] = numbers::invalid_unsigned_int;
-
-  if ((is_face &&
-       this->dof_access_index ==
-         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
-       this->is_interior_face == false) ||
-      (!is_face && !this->is_interior_face))
-    {
-      // cell-based face-loop: plus face
-      for (unsigned int i = 0; i < v_len; ++i)
-        {
-          // compute actual (non vectorized) cell ID
-          const unsigned int cell_this = this->cell * v_len + i;
-          // compute face ID
-          unsigned int fn =
-            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
-                                                                  this->face_no,
-                                                                  i);
-
-          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);
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_quadrature_index() const
-{
-  return this->quad_no;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_current_cell_index() const
-{
-  if (is_face && this->dof_access_index ==
-                   internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
-    return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
-  else
-    return this->cell;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_active_fe_index() const
-{
-  return active_fe_index;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline unsigned int
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_active_quadrature_index() const
-{
-  return active_quad_index;
-}
-
-
-
-template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
-inline const MatrixFree<dim, Number, VectorizedArrayType> &
-FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
-  get_matrix_free() const
-{
-  Assert(matrix_info != nullptr,
-         ExcMessage(
-           "FEEvaluation was not initialized with a MatrixFree object!"));
-  return *matrix_info;
-}
+           array[cells[i] / VectorizedArrayType::size()]
+                [cells[i] % VectorizedArrayType::size()]);
+  }
+} // namespace internal
 
 
-/*----------------------- 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,
-                   const unsigned int active_fe_index,
-                   const unsigned int active_quad_index,
-                   const unsigned int face_type)
-  : FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>(
-      data_in,
-      dof_no,
-      first_selected_component,
-      quad_no_in,
-      fe_degree,
-      n_q_points,
-      is_interior_face,
-      active_fe_index,
-      active_quad_index,
-      face_type)
-  , 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)
+std::array<unsigned int, VectorizedArrayType::size()>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  get_cell_or_face_ids() const
 {
-  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) + ")"));
+  const unsigned int v_len = VectorizedArrayType::size();
+  std::array<unsigned int, VectorizedArrayType::size()> cells;
 
-  // do not check for correct dimensions of data fields here, should be done
-  // in derived classes
+  // initialize array
+  for (unsigned int i = 0; i < v_len; ++i)
+    cells[i] = numbers::invalid_unsigned_int;
+
+  if ((is_face &&
+       this->dof_access_index ==
+         internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+       this->is_interior_face == false) ||
+      (!is_face && !this->is_interior_face))
+    {
+      // cell-based face-loop: plus face
+      for (unsigned int i = 0; i < v_len; ++i)
+        {
+          // compute actual (non vectorized) cell ID
+          const unsigned int cell_this = this->cell * v_len + i;
+          // compute face ID
+          unsigned int fn =
+            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+                                                                  this->face_no,
+                                                                  i);
+
+          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] = this->cell * v_len + i;
+    }
+
+  return cells;
 }
 
 
@@ -4226,48 +3721,16 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-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)
+inline VectorizedArrayType
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  read_cell_data(const AlignedVector<VectorizedArrayType> &array) const
 {
-  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;
+  VectorizedArrayType out = Number(1.);
+  internal::process_cell_data(
+    *this, this->matrix_info, array, out, [](auto &local, const auto &global) {
+      local = global;
+    });
+  return out;
 }
 
 
@@ -4277,27 +3740,15 @@ template <int dim,
           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)
-  : 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)
+inline void
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  set_cell_data(AlignedVector<VectorizedArrayType> &array,
+                const VectorizedArrayType &         in) const
 {
-  set_data_pointers();
+  internal::process_cell_data(
+    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+      global = local;
+    });
 }
 
 
@@ -4307,24 +3758,18 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
-inline FEEvaluationBase<dim,
-                        n_components_,
-                        Number,
-                        is_face,
-                        VectorizedArrayType> &
+template <typename T>
+inline std::array<T, VectorizedArrayType::size()>
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-operator=(const FEEvaluationBase<dim,
-                                 n_components_,
-                                 Number,
-                                 is_face,
-                                 VectorizedArrayType> &other)
+  read_cell_data(const AlignedVector<std::array<T, VectorizedArrayType::size()>>
+                   &array) const
 {
-  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;
+  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;
 }
 
 
@@ -4334,48 +3779,17 @@ template <int dim,
           typename Number,
           bool is_face,
           typename VectorizedArrayType>
+template <typename T>
 inline void
 FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
-  set_data_pointers()
+  set_cell_data(
+    AlignedVector<std::array<T, VectorizedArrayType::size()>> &array,
+    const std::array<T, VectorizedArrayType::size()> &         in) const
 {
-  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 = this->n_quadrature_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 + 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);
-  gradients_from_hessians_quad =
-    this->scratch_data_array->begin() +
-    n_components * (dofs_per_component + (dim + 1) * n_quadrature_points);
-  hessians_quad =
-    this->scratch_data_array->begin() +
-    n_components * (dofs_per_component + (2 * dim + 1) * n_quadrature_points);
-  this->scratch_data = this->scratch_data_array->begin() +
-                       n_components_ * dofs_per_component +
-                       (n_components_ * ((dim * (dim + 1)) / 2 + 2 * dim + 1) *
-                        n_quadrature_points);
+  internal::process_cell_data(
+    *this, this->matrix_info, array, in, [](const auto &local, auto &global) {
+      global = local;
+    });
 }
 
 
@@ -6571,7 +5985,7 @@ inline FEEvaluationAccess<dim,
     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 quad_no,
     const unsigned int fe_degree,
     const unsigned int n_q_points,
     const bool         is_interior_face,
@@ -6582,7 +5996,7 @@ inline FEEvaluationAccess<dim,
       data_in,
       dof_no,
       first_selected_component,
-      quad_no_in,
+      quad_no,
       fe_degree,
       n_q_points,
       is_interior_face,
@@ -6679,7 +6093,7 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
     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 quad_no,
     const unsigned int fe_degree,
     const unsigned int n_q_points,
     const bool         is_interior_face,
@@ -6690,7 +6104,7 @@ inline FEEvaluationAccess<dim, 1, Number, is_face, VectorizedArrayType>::
       data_in,
       dof_no,
       first_selected_component,
-      quad_no_in,
+      quad_no,
       fe_degree,
       n_q_points,
       is_interior_face,
@@ -7006,7 +6420,7 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
     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 quad_no,
     const unsigned int fe_degree,
     const unsigned int n_q_points,
     const bool         is_interior_face,
@@ -7017,7 +6431,7 @@ inline FEEvaluationAccess<dim, dim, Number, is_face, VectorizedArrayType>::
       data_in,
       dof_no,
       first_selected_component,
-      quad_no_in,
+      quad_no,
       fe_degree,
       n_q_points,
       is_interior_face,
@@ -7417,7 +6831,7 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
   FEEvaluationAccess(const MatrixFree<1, Number, VectorizedArrayType> &data_in,
                      const unsigned int                                dof_no,
                      const unsigned int first_selected_component,
-                     const unsigned int quad_no_in,
+                     const unsigned int quad_no,
                      const unsigned int fe_degree,
                      const unsigned int n_q_points,
                      const bool         is_interior_face,
@@ -7428,7 +6842,7 @@ inline FEEvaluationAccess<1, 1, Number, is_face, VectorizedArrayType>::
       data_in,
       dof_no,
       first_selected_component,
-      quad_no_in,
+      quad_no,
       fe_degree,
       n_q_points,
       is_interior_face,
@@ -8377,7 +7791,7 @@ FEEvaluation<dim,
                this->values_quad,
                this->gradients_quad,
                this->hessians_quad,
-               this->scratch_data);
+               this->scratch_data.begin());
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::evaluate(
       n_components,
@@ -8387,7 +7801,7 @@ FEEvaluation<dim,
       this->values_quad,
       this->gradients_quad,
       this->hessians_quad,
-      this->scratch_data);
+      this->scratch_data.begin());
 
   if (hessians_on_general_cells)
     for (unsigned int i = 0; i < n_components * dim * n_q_points; ++i)
@@ -8808,7 +8222,7 @@ FEEvaluation<dim,
                 this->values_quad,
                 this->gradients_quad,
                 this->hessians_quad,
-                this->scratch_data,
+                this->scratch_data.begin(),
                 false);
   else
     internal::FEEvaluationFactory<dim, Number, VectorizedArrayType>::integrate(
@@ -8819,7 +8233,7 @@ FEEvaluation<dim,
       this->values_quad,
       this->gradients_quad,
       this->hessians_quad,
-      this->scratch_data,
+      this->scratch_data.begin(),
       false);
 
 #  ifdef DEBUG
@@ -8888,7 +8302,7 @@ FEEvaluation<dim,
         this->values_quad,
         this->gradients_quad,
         this->hessians_quad,
-        this->scratch_data,
+        this->scratch_data.begin(),
         this->data,
         integration_flag) == false)
     {
@@ -9256,7 +8670,7 @@ FEFaceEvaluation<dim,
           this->begin_values(),
           this->begin_gradients(),
           this->begin_hessians(),
-          this->scratch_data,
+          this->scratch_data.begin(),
           evaluation_flag_actual & EvaluationFlags::values,
           evaluation_flag_actual & EvaluationFlags::gradients,
           evaluation_flag_actual & EvaluationFlags::hessians,
@@ -9277,7 +8691,7 @@ FEFaceEvaluation<dim,
             this->begin_values(),
             this->begin_gradients(),
             this->begin_hessians(),
-            this->scratch_data,
+            this->scratch_data.begin(),
             evaluation_flag_actual & EvaluationFlags::values,
             evaluation_flag_actual & EvaluationFlags::gradients,
             evaluation_flag_actual & EvaluationFlags::hessians,
@@ -9293,7 +8707,7 @@ FEFaceEvaluation<dim,
                    this->begin_values(),
                    this->begin_gradients(),
                    this->begin_hessians(),
-                   this->scratch_data,
+                   this->scratch_data.begin(),
                    evaluation_flag_actual & EvaluationFlags::values,
                    evaluation_flag_actual & EvaluationFlags::gradients,
                    evaluation_flag_actual & EvaluationFlags::hessians,
@@ -9451,7 +8865,7 @@ FEFaceEvaluation<dim,
         this->begin_values(),
         this->begin_gradients(),
         this->begin_hessians(),
-        this->scratch_data,
+        this->scratch_data.begin(),
         evaluation_flag_actual & EvaluationFlags::values,
         evaluation_flag_actual & EvaluationFlags::gradients,
         evaluation_flag_actual & EvaluationFlags::hessians,
@@ -9467,7 +8881,7 @@ FEFaceEvaluation<dim,
                 this->begin_values(),
                 this->begin_gradients(),
                 this->begin_hessians(),
-                this->scratch_data,
+                this->scratch_data.begin(),
                 evaluation_flag_actual & EvaluationFlags::values,
                 evaluation_flag_actual & EvaluationFlags::gradients,
                 evaluation_flag_actual & EvaluationFlags::hessians,
@@ -9560,7 +8974,7 @@ FEFaceEvaluation<dim,
           this->begin_values(),
           this->begin_gradients(),
           this->begin_hessians(),
-          this->scratch_data,
+          this->scratch_data.begin(),
           evaluation_flag & EvaluationFlags::values,
           evaluation_flag & EvaluationFlags::gradients,
           evaluation_flag & EvaluationFlags::hessians,
@@ -9602,7 +9016,7 @@ FEFaceEvaluation<dim,
               this->begin_values(),
               this->begin_gradients(),
               this->begin_hessians(),
-              this->scratch_data,
+              this->scratch_data.begin(),
               evaluation_flag & EvaluationFlags::values,
               evaluation_flag & EvaluationFlags::gradients,
               evaluation_flag & EvaluationFlags::hessians,
@@ -9627,7 +9041,7 @@ FEFaceEvaluation<dim,
                               this->begin_values(),
                               this->begin_gradients(),
                               this->begin_hessians(),
-                              this->scratch_data,
+                              this->scratch_data.begin(),
                               evaluation_flag & EvaluationFlags::values,
                               evaluation_flag & EvaluationFlags::gradients,
                               evaluation_flag & EvaluationFlags::hessians,
@@ -9744,7 +9158,7 @@ FEFaceEvaluation<dim,
             this->begin_values(),
             this->begin_gradients(),
             this->begin_hessians(),
-            this->scratch_data,
+            this->scratch_data.begin(),
             evaluation_flag & EvaluationFlags::values,
             evaluation_flag & EvaluationFlags::gradients,
             evaluation_flag & EvaluationFlags::hessians,
@@ -9777,7 +9191,7 @@ FEFaceEvaluation<dim,
                               this->begin_values(),
                               this->begin_gradients(),
                               this->begin_hessians(),
-                              this->scratch_data,
+                              this->scratch_data.begin(),
                               evaluation_flag & EvaluationFlags::values,
                               evaluation_flag & EvaluationFlags::gradients,
                               evaluation_flag & EvaluationFlags::hessians,
diff --git a/include/deal.II/matrix_free/fe_evaluation_base_data.h b/include/deal.II/matrix_free/fe_evaluation_base_data.h
new file mode 100644 (file)
index 0000000..dcce525
--- /dev/null
@@ -0,0 +1,763 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_fe_evaluation_base_data_h
+#define dealii_matrix_free_fe_evaluation_base_data_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/array_view.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/smartpointer.h>
+#include <deal.II/base/symmetric_tensor.h>
+#include <deal.II/base/template_constraints.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/dof_info.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
+#include <deal.II/matrix_free/shape_info.h>
+#include <deal.II/matrix_free/type_traits.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+
+namespace internal
+{
+  DeclException0(ExcAccessToUninitializedField);
+
+  DeclException1(
+    ExcMatrixFreeAccessToUninitializedMappingField,
+    std::string,
+    << "You are requesting information from an FEEvaluation/FEFaceEvaluation "
+    << "object for which this kind of information has not been computed. What "
+    << "information these objects compute is determined by the update_* flags you "
+    << "pass to MatrixFree::reinit() via MatrixFree::AdditionalData. Here, "
+    << "the operation you are attempting requires the <" << arg1
+    << "> flag to be set, but it was apparently not specified "
+    << "upon initialization.");
+} // namespace internal
+
+// forward declarations
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d            = fe_degree + 1,
+          int n_components_            = 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+class FEEvaluation;
+
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool     = false,
+          typename = VectorizedArray<Number>>
+class FEEvaluationBase;
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    template <int, typename, typename>
+    class MappingDataOnTheFly;
+  }
+} // namespace internal
+
+
+/**
+ * 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 four template arguments:
+ *
+ * @tparam dim Dimension in which this class is to be used
+ *
+ * @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, typename Number, bool is_face, typename VectorizedArrayType>
+class FEEvaluationBaseData
+{
+  static_assert(
+    std::is_same<Number, typename VectorizedArrayType::value_type>::value,
+    "Type of Number and of VectorizedArrayType do not match.");
+
+  using ShapeInfoType =
+    internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType>;
+  using MappingInfoStorageType =
+    internal::MatrixFreeFunctions::MappingInfoStorage<(is_face ? dim - 1 : dim),
+                                                      dim,
+                                                      Number,
+                                                      VectorizedArrayType>;
+  using DoFInfo = internal::MatrixFreeFunctions::DoFInfo;
+
+public:
+  static constexpr unsigned int dimension = dim;
+
+  /**
+   * Return the index offset within the geometry fields for the cell the @p
+   * reinit() function has been called for. This index can be used to access
+   * an index into a field that has the same compression behavior as the
+   * Jacobian of the geometry, e.g., to store an effective coefficient tensors
+   * that combines a coefficient with the geometry for lower memory transfer
+   * as the available data fields.
+   */
+  unsigned int
+  get_mapping_data_index_offset() const;
+
+  /**
+   * Return the type of the cell the @p reinit() function has been called for.
+   * Valid values are @p cartesian for Cartesian cells (which allows for
+   * considerable data compression), @p affine for cells with affine mappings,
+   * and @p general for general cells without any compressed storage applied.
+   */
+  internal::MatrixFreeFunctions::GeometryType
+  get_cell_type() const;
+
+  /**
+   * Return a reference to the ShapeInfo object currently in use.
+   */
+  const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+  get_shape_info() const;
+
+  /**
+   * Return a reference to the DoFInfo object currently in use.
+   */
+  const internal::MatrixFreeFunctions::DoFInfo &
+  get_dof_info() 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;
+
+  /**
+   * 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 $J^{-\mathrm T}$ of the
+   * Jacobian of the mapping between the unit to the real cell defined as
+   * $J_{ij} = d x_i / d\hat x_j$. The $(i,j)$ entry of the returned tensor
+   * contains $d\hat x_j/dx_i$, i.e., columns refer to reference space
+   * coordinates and rows to real cell coordinates. Thus, the returned tensor
+   * represents a covariant transformation, which is used in the
+   * FEEvaluationBase::get_gradient() function to transform the unit cell
+   * gradients to gradients on the real cell by a multiplication $J^{-\mathrm
+   * T} \hat{\nabla} u_h$.
+   */
+  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;
+
+  /**
+   * 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 the number of the quadrature formula of the present cell.
+   */
+  unsigned int
+  get_quadrature_index() const;
+
+  /**
+   * Return index of the current cell or face.
+   */
+  unsigned int
+  get_current_cell_index() const;
+
+  /**
+   * Return the active FE index for this class for efficient indexing in the
+   * hp-case.
+   */
+  unsigned int
+  get_active_fe_index() const;
+
+  /**
+   * Return the active quadrature index for this class for efficient indexing in
+   * the hp-case.
+   */
+  unsigned int
+  get_active_quadrature_index() 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 ShapeInfoType &         shape_info,
+                       const DoFInfo &               dof_info,
+                       const MappingInfoStorageType &mapping_data,
+                       const unsigned int            quad_no,
+                       const bool                    is_interior_face,
+                       const unsigned int            active_fe_index,
+                       const unsigned int            active_quad_index,
+                       const unsigned int            face_type);
+
+  /**
+   * Constructor that comes with reduced functionality and works similar as
+   * FEValues.
+   */
+  FEEvaluationBaseData(
+    const std::shared_ptr<
+      internal::MatrixFreeFunctions::
+        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapping_data);
+
+  /**
+   * Copy constructor.
+   */
+  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);
+
+  /**
+   * 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 ShapeInfoType *data;
+
+  /**
+   * 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 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 MappingInfoStorageType *mapping_data;
+
+  /**
+   * The number of the quadrature formula of the present cell among all
+   * quadrature formulas available in the MatrixFree objects pass to derived
+   * classes.
+   */
+  const unsigned int quad_no;
+
+  /**
+   * 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;
+
+  /**
+   * A pointer to the underlying quadrature formula specified at construction.
+   * Also contained in mapping_data, but it simplifies code if we store a
+   * reference to it.
+   */
+  const typename MappingInfoStorageType::QuadratureDescriptor *descriptor;
+
+  /**
+   * The number of quadrature points in the current evaluation context.
+   */
+  const unsigned int n_quadrature_points;
+
+  /**
+   * 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;
+
+  /**
+   * This is the user-visible part of FEEvaluationBase::scratch_data_array,
+   * only showing the part that can be consumed by various users. It is set
+   * during the allocation of the internal data structures in
+   * FEEvaluationBase.
+   */
+  mutable ArrayView<VectorizedArrayType> scratch_data;
+
+  /**
+   * 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. For cells it defines if
+   * the dof values should be read from the actual cell corresponding to the
+   * interior face or the neighboring cell corresponding to the exterior face.
+   */
+  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, as an alternative to the entry point with
+   * MatrixFree.
+   */
+  std::shared_ptr<internal::MatrixFreeFunctions::
+                    MappingDataOnTheFly<dim, Number, VectorizedArrayType>>
+    mapped_geometry;
+
+  // Make FEEvaluation and FEEvaluationBase objects friends for access to
+  // protected member mapped_geometry.
+  template <int, int, typename, bool, typename>
+  friend class FEEvaluationBase;
+
+  template <int, int, int, int, typename, typename>
+  friend class FEEvaluation;
+};
+
+
+
+/*----------------------- Inline functions ----------------------------------*/
+
+#ifndef DOXYGEN
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(const ShapeInfoType &         shape_info,
+                       const DoFInfo &               dof_info,
+                       const MappingInfoStorageType &mapping_data,
+                       const unsigned int            quad_no,
+                       const bool                    is_interior_face,
+                       const unsigned int            active_fe_index,
+                       const unsigned int            active_quad_index,
+                       const unsigned int            face_type)
+  : data(&shape_info)
+  , dof_info(&dof_info)
+  , mapping_data(&mapping_data)
+  , quad_no(quad_no)
+  , active_fe_index(active_fe_index)
+  , active_quad_index(active_quad_index)
+  , descriptor(
+      &mapping_data.descriptor
+         [is_face ?
+            (active_quad_index * std::max<unsigned int>(1, dim - 1) +
+             (face_type == numbers::invalid_unsigned_int ? 0 : face_type)) :
+            active_quad_index])
+  , n_quadrature_points(descriptor->n_q_points)
+  , jacobian(nullptr)
+  , J_value(nullptr)
+  , normal_vectors(nullptr)
+  , normal_x_jacobian(nullptr)
+  , quadrature_weights(descriptor->quadrature_weights.begin())
+  , cell(numbers::invalid_unsigned_int)
+  , is_interior_face(is_interior_face)
+  , dof_access_index(
+      is_face ?
+        (is_interior_face ?
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+  , face_no(0)
+  , face_orientation(0)
+  , subface_index(0)
+  , cell_type(internal::MatrixFreeFunctions::general)
+{}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(
+    const std::shared_ptr<
+      internal::MatrixFreeFunctions::
+        MappingDataOnTheFly<dim, Number, VectorizedArrayType>> &mapped_geometry)
+  : data(nullptr)
+  , dof_info(nullptr)
+  , mapping_data(nullptr)
+  , quad_no(numbers::invalid_unsigned_int)
+  , active_fe_index(numbers::invalid_unsigned_int)
+  , active_quad_index(numbers::invalid_unsigned_int)
+  , descriptor(nullptr)
+  , n_quadrature_points(
+      mapped_geometry->get_data_storage().descriptor[0].n_q_points)
+  , jacobian(nullptr)
+  , 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)
+  , mapped_geometry(mapped_geometry)
+{
+  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();
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  FEEvaluationBaseData(const FEEvaluationBaseData &other)
+  : data(other.data)
+  , dof_info(other.dof_info)
+  , mapping_data(other.mapping_data)
+  , quad_no(other.quad_no)
+  , active_fe_index(other.active_fe_index)
+  , active_quad_index(other.active_quad_index)
+  , descriptor(other.descriptor)
+  , n_quadrature_points(other.n_quadrature_points)
+  , jacobian(nullptr)
+  , J_value(nullptr)
+  , normal_vectors(nullptr)
+  , normal_x_jacobian(nullptr)
+  , quadrature_weights(other.quadrature_weights)
+  , cell(numbers::invalid_unsigned_int)
+  , is_interior_face(other.is_interior_face)
+  , dof_access_index(
+      is_face ?
+        (is_interior_face ?
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+  , face_no(0)
+  , face_orientation(0)
+  , subface_index(0)
+  , cell_type(internal::MatrixFreeFunctions::general)
+{}
+
+
+
+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 &other)
+{
+  AssertDimension(quad_no, other.quad_no);
+  AssertDimension(active_fe_index, other.active_fe_index);
+  AssertDimension(active_quad_index, other.active_quad_index);
+  AssertDimension(n_quadrature_points, descriptor->n_q_points);
+
+  data               = other.data;
+  dof_info           = other.dof_info;
+  mapping_data       = other.mapping_data;
+  descriptor         = other.descriptor;
+  jacobian           = nullptr;
+  J_value            = nullptr;
+  normal_vectors     = nullptr;
+  normal_x_jacobian  = nullptr;
+  quadrature_weights = other.quadrature_weights;
+  cell               = numbers::invalid_unsigned_int;
+  is_interior_face   = other.is_interior_face;
+  dof_access_index =
+    is_face ?
+      (is_interior_face ?
+         internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
+         internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
+      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
+  face_no          = 0;
+  face_orientation = 0;
+  subface_index    = 0;
+  cell_type        = internal::MatrixFreeFunctions::general;
+
+  return *this;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_mapping_data_index_offset() const
+{
+  Assert(mapping_data != nullptr, ExcInternalError());
+
+  if (dof_info == nullptr)
+    return 0;
+  else
+    {
+      AssertIndexRange(cell, mapping_data->data_index_offsets.size());
+      return mapping_data->data_index_offsets[cell];
+    }
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline internal::MatrixFreeFunctions::GeometryType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_cell_type()
+  const
+{
+  Assert(cell != numbers::invalid_unsigned_int, ExcNotInitialized());
+  return cell_type;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const internal::MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_shape_info() const
+{
+  Assert(data != nullptr, ExcInternalError());
+  return *data;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline const internal::MatrixFreeFunctions::DoFInfo &
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::get_dof_info()
+  const
+{
+  Assert(dof_info != nullptr,
+         ExcMessage(
+           "FEEvaluation was not initialized with a MatrixFree object!"));
+  return *dof_info;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_normal_vector(const unsigned int q_point) const
+{
+  AssertIndexRange(q_point, n_quadrature_points);
+  Assert(normal_vectors != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_normal_vectors"));
+  if (this->cell_type <= internal::MatrixFreeFunctions::flat_faces)
+    return normal_vectors[0];
+  else
+    return normal_vectors[q_point];
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE VectorizedArrayType
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::JxW(
+  const unsigned int q_point) const
+{
+  AssertIndexRange(q_point, n_quadrature_points);
+  Assert(J_value != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_values|update_gradients"));
+  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+    {
+      Assert(this->quadrature_weights != nullptr, ExcInternalError());
+      return J_value[0] * this->quadrature_weights[q_point];
+    }
+  else
+    return J_value[q_point];
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  inverse_jacobian(const unsigned int q_point) const
+{
+  AssertIndexRange(q_point, n_quadrature_points);
+  Assert(this->jacobian != nullptr,
+         internal::ExcMatrixFreeAccessToUninitializedMappingField(
+           "update_gradients"));
+  if (this->cell_type <= internal::MatrixFreeFunctions::affine)
+    return jacobian[0];
+  else
+    return jacobian[q_point];
+}
+
+
+
+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 unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_quadrature_index() const
+{
+  return this->quad_no;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_current_cell_index() const
+{
+  if (is_face && this->dof_access_index ==
+                   internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
+    return this->cell * GeometryInfo<dim>::faces_per_cell + this->face_no;
+  else
+    return this->cell;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_active_fe_index() const
+{
+  return active_fe_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline unsigned int
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_active_quadrature_index() const
+{
+  return active_quad_index;
+}
+
+
+
+template <int dim, typename Number, bool is_face, typename VectorizedArrayType>
+inline ArrayView<VectorizedArrayType>
+FEEvaluationBaseData<dim, Number, is_face, VectorizedArrayType>::
+  get_scratch_data() const
+{
+  return scratch_data;
+}
+
+
+
+#endif // ifndef DOXYGEN
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index a97a97252ba3e3e51fd17967c7e1029fede46d04..4fc62e48bea8c9ab113bcb9b02b72b58c0aac754 100644 (file)
@@ -34,6 +34,7 @@
 
 #include <deal.II/matrix_free/face_info.h>
 #include <deal.II/matrix_free/helper_functions.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
 
 #include <memory>
 
@@ -45,285 +46,6 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
-    /**
-     * An enum to identify various types of cells and faces. The most general
-     * type is what we typically compute in the FEValues context but for many
-     * geometries we can save significant storage.
-     *
-     * @ingroup matrixfree
-     */
-    enum GeometryType : unsigned char
-    {
-      /**
-       * The cell or face is Cartesian.
-       */
-      cartesian = 0,
-
-      /**
-       * The cell or face can be described with an affine mapping.
-       */
-      affine = 1,
-
-      /**
-       * The face is flat, i.e., the normal factor on a face is the same on
-       * all quadrature points. This type is not assigned for cells.
-       */
-      flat_faces = 2,
-
-      /**
-       * There is no special information available for compressing the
-       * representation of the object under consideration.
-       */
-      general = 3
-    };
-
-
-
-    /**
-     * Definition of a structure that stores all cached data related to the
-     * evaluated geometry from the mapping. In order to support hp-adaptivity
-     * and compressed storage (in particular for Jacobians, JxW values, and
-     * normals), storage length can be different for different rows. Thus, it
-     * allows to jump at the data of individual rows similar to compressed row
-     * storage in sparse matrices. We have two different start indices for
-     * fields with different sizes. The first category of offsets are the
-     * indices for Jacobians of the transformation from unit to real cell (we
-     * store the inverse Jacobian), second derivatives, JxW values, and normal
-     * vectors. We keep separate arrays for all these data structures because
-     * a user code might access only some of them. In such a case, one array
-     * will be gone through in a contiguous order with access to all entries,
-     * which makes it easy for the processor to prefetch data. Having all data
-     * in a single array would require some strides in the access pattern,
-     * which is much more complicated for the processor to predict (and indeed
-     * leads to prefetching of data that does not get used on Intel processors
-     * such as BroadwellEP).
-     *
-     * The second category of indices are the offsets for the quadrature
-     * points. Quadrature points can be compressed less than the other fields
-     * and thus need longer fields. Quadrature point indices are often used in
-     * other contexts such as evaluation of right hand sides.
-     *
-     * The third component is a descriptor of data from the unit cells, called
-     * QuadratureDescriptor, which contains the quadrature weights and
-     * permutations of how to go through quadrature points in case of face
-     * data. The latter comes in a vector for the support of hp-adaptivity,
-     * with several data fields for the individual quadrature formulas.
-     *
-     * @ingroup matrixfree
-     */
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    struct MappingInfoStorage
-    {
-      static_assert(
-        std::is_same<Number, typename VectorizedArrayType::value_type>::value,
-        "Type of Number and of VectorizedArrayType do not match.");
-
-      struct QuadratureDescriptor
-      {
-        /**
-         * Constructor. Does nothing.
-         */
-        QuadratureDescriptor();
-
-        /**
-         * Set up the lengths in the various members of this struct.
-         */
-        template <int dim_q>
-        void
-        initialize(const Quadrature<dim_q> &quadrature,
-                   const UpdateFlags update_flags_inner_faces = update_default);
-
-        /**
-         * Set up the lengths in the various members of this struct.
-         */
-        void
-        initialize(const Quadrature<1> &quadrature_1d,
-                   const UpdateFlags update_flags_inner_faces = update_default);
-
-        /**
-         * Returns the memory consumption in bytes.
-         */
-        std::size_t
-        memory_consumption() const;
-
-        /**
-         * Number of quadrature points applied on the given cell or face.
-         */
-        unsigned int n_q_points;
-
-        /**
-         * Original one-dimensional quadrature formula applied on the given
-         * cell or face.
-         */
-        Quadrature<1> quadrature_1d;
-
-        /**
-         * Quadrature formula applied on the given cell or face.
-         */
-        Quadrature<structdim> quadrature;
-
-        /**
-         * Quadrature weights separated by dimension for use in specific
-         * situations.
-         */
-        std::array<AlignedVector<Number>, structdim> tensor_quadrature_weights;
-
-        /**
-         * A cached vector of quadrature weights in the given number format
-         * (non-vectorized, as it is cheap to broadcast the value to all lanes
-         * when it is used in a vectorized context).
-         */
-        AlignedVector<Number> quadrature_weights;
-
-        /**
-         * For quadrature on faces, the evaluation of basis functions is not
-         * in the correct order if a face is not in the standard orientation
-         * to a given element. This data structure is used to re-order the
-         * data evaluated on quadrature points to represent the correct order.
-         */
-        dealii::Table<2, unsigned int> face_orientations;
-      };
-
-      /**
-       * A class describing the layout of the sections in the @p data_storage
-       * field and also includes some data that depends on the number of
-       * quadrature points in the hp-context such as the inner quadrature
-       * formula and re-indexing for faces that are not in the standard
-       * orientation.
-       */
-      std::vector<QuadratureDescriptor> descriptor;
-
-      /**
-       * Collection of quadrature formulae applied on the given face.
-       *
-       * @note Only filled for faces, since faces might be quadrilateral or
-       *   triangle shaped.
-       */
-      std::vector<dealii::hp::QCollection<structdim>> q_collection;
-
-      /**
-       * Stores the index offset into the arrays @p jxw_values, @p jacobians,
-       * @p normal_vectors and the second derivatives. Note that affine cells
-       * have shorter fields of length 1, where the others have lengths equal
-       * to the number of quadrature points of the given cell.
-       */
-      AlignedVector<unsigned int> data_index_offsets;
-
-      /**
-       * The storage of the Jacobian determinant (times the quadrature weight
-       * in case the transformation is non-affine) on quadrature
-       * points.
-       *
-       * Indexed by @p data_index_offsets.
-       */
-      AlignedVector<VectorizedArrayType> JxW_values;
-
-      /**
-       * Stores the normal vectors.
-       *
-       * Indexed by @p data_index_offsets.
-       */
-      AlignedVector<Tensor<1, spacedim, VectorizedArrayType>> normal_vectors;
-
-      /**
-       * The storage of covariant transformation on quadrature points, i.e.,
-       * the inverse and transposed Jacobians of the transformation from the
-       * unit to the real cell.
-       *
-       * Indexed by @p data_index_offsets.
-       *
-       * Contains two fields for access from both sides for interior faces,
-       * but the default case (cell integrals or boundary integrals) only
-       * fills the zeroth component and ignores the first one.
-       */
-      std::array<AlignedVector<Tensor<2, spacedim, VectorizedArrayType>>, 2>
-        jacobians;
-
-      /**
-       * The storage of the gradients of the inverse Jacobian
-       * transformation. Because of symmetry, only the upper diagonal and
-       * diagonal part are needed. The first index runs through the
-       * derivatives, starting with the diagonal and then continuing row-wise,
-       * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then
-       * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index
-       * is the spatial coordinate.
-       *
-       * Indexed by @p data_index_offsets.
-       *
-       * Contains two fields for access from both sides for interior faces,
-       * but the default case (cell integrals or boundary integrals) only
-       * fills the zeroth component and ignores the first one.
-       */
-      std::array<
-        AlignedVector<Tensor<1,
-                             spacedim *(spacedim + 1) / 2,
-                             Tensor<1, spacedim, VectorizedArrayType>>>,
-        2>
-        jacobian_gradients;
-
-      /**
-       * Stores the Jacobian transformations times the normal vector (this
-       * represents a shortcut that is accessed often and can thus get higher
-       * performance).
-       *
-       * Indexed by @p data_index_offsets.
-       */
-      std::array<AlignedVector<Tensor<1, spacedim, VectorizedArrayType>>, 2>
-        normals_times_jacobians;
-
-      /**
-       * Stores the index offset of a particular cell into the quadrature
-       * points array in real coordinates. Note that Cartesian cells have
-       * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells
-       * (length is @p n_q_points) or faces.
-       */
-      AlignedVector<unsigned int> quadrature_point_offsets;
-
-      /**
-       * Stores the quadrature points in real coordinates, including a
-       * compression scheme for Cartesian cells where we do not need to store
-       * the full data on all points.
-       *
-       * Indexed by @p quadrature_point_offsets.
-       */
-      AlignedVector<Point<spacedim, VectorizedArrayType>> quadrature_points;
-
-      /**
-       * Clears all data fields except the descriptor vector.
-       */
-      void
-      clear_data_fields();
-
-      /**
-       * Returns the quadrature index for a given number of quadrature
-       * points. If not in hp-mode or if the index is not found, this
-       * function always returns index 0. Hence, this function does not
-       * check whether the given degree is actually present.
-       */
-      unsigned int
-      quad_index_from_n_q_points(const unsigned int n_q_points) const;
-
-      /**
-       * Prints a detailed summary of memory consumption in the different
-       * structures of this class to the given output stream.
-       */
-      template <typename StreamType>
-      void
-      print_memory_consumption(StreamType &    out,
-                               const TaskInfo &task_info) const;
-
-      /**
-       * Returns the memory consumption in bytes.
-       */
-      std::size_t
-      memory_consumption() const;
-    };
-
-
-
     /**
      * The class that stores all geometry-dependent data related with cell
      * interiors for use in the matrix-free class.
@@ -555,12 +277,12 @@ namespace internal
     template <int dim, typename Number, typename VectorizedArrayType>
     struct MappingInfoCellsOrFaces<dim, Number, false, VectorizedArrayType>
     {
-      static const MappingInfoStorage<dim, dim, Number, VectorizedArrayType> *
+      static const MappingInfoStorage<dim, dim, Number, VectorizedArrayType> &
       get(const MappingInfo<dim, Number, VectorizedArrayType> &mapping_info,
           const unsigned int                                   quad_no)
       {
         AssertIndexRange(quad_no, mapping_info.cell_data.size());
-        return &mapping_info.cell_data[quad_no];
+        return mapping_info.cell_data[quad_no];
       }
     };
 
@@ -568,12 +290,12 @@ namespace internal
     struct MappingInfoCellsOrFaces<dim, Number, true, VectorizedArrayType>
     {
       static const MappingInfoStorage<dim - 1, dim, Number, VectorizedArrayType>
-        *
+        &
         get(const MappingInfo<dim, Number, VectorizedArrayType> &mapping_info,
             const unsigned int                                   quad_no)
       {
         AssertIndexRange(quad_no, mapping_info.face_data.size());
-        return &mapping_info.face_data[quad_no];
+        return mapping_info.face_data[quad_no];
       }
     };
 
@@ -651,22 +373,6 @@ namespace internal
 
     /* ------------------- inline functions ----------------------------- */
 
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    inline unsigned int
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      quad_index_from_n_q_points(const unsigned int n_q_points) const
-    {
-      for (unsigned int i = 0; i < descriptor.size(); ++i)
-        if (n_q_points == descriptor[i].n_q_points)
-          return i;
-      return 0;
-    }
-
-
-
     template <int dim, typename Number, typename VectorizedArrayType>
     inline GeometryType
     MappingInfo<dim, Number, VectorizedArrayType>::get_cell_type(
index 2c3072008f27b35e8fe281a7423c8fcb7e22e23d..d5c35a4972f68480bafdcf5b69f2fc1c701e43f2 100644 (file)
@@ -30,6 +30,7 @@
 
 #include <deal.II/matrix_free/evaluation_template_factory.h>
 #include <deal.II/matrix_free/mapping_info.h>
+#include <deal.II/matrix_free/mapping_info_storage.templates.h>
 #include <deal.II/matrix_free/util.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -39,231 +40,6 @@ namespace internal
 {
   namespace MatrixFreeFunctions
   {
-    /* ------------------------ MappingInfoStorage implementation ---------- */
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      QuadratureDescriptor::QuadratureDescriptor()
-      : n_q_points(numbers::invalid_unsigned_int)
-    {}
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    template <int dim_q>
-    void
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      QuadratureDescriptor::initialize(
-        const Quadrature<dim_q> &quadrature,
-        const UpdateFlags        update_flags_inner_faces)
-    {
-      Assert(structdim + 1 <= spacedim ||
-               update_flags_inner_faces == update_default,
-             ExcMessage("Volume cells do not allow for setting inner faces"));
-      this->quadrature = quadrature;
-      n_q_points       = quadrature.size();
-      quadrature_weights.resize(n_q_points);
-      for (unsigned int i = 0; i < n_q_points; ++i)
-        quadrature_weights[i] = quadrature.weight(i);
-
-      // note: quadrature_1d and tensor_quadrature_weights are not set up
-
-      // TODO: set up face_orientations
-      (void)update_flags_inner_faces;
-    }
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    void
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      QuadratureDescriptor::initialize(
-        const Quadrature<1> &quadrature_1d,
-        const UpdateFlags    update_flags_inner_faces)
-    {
-      Assert(structdim + 1 <= spacedim ||
-               update_flags_inner_faces == update_default,
-             ExcMessage("Volume cells do not allow for setting inner faces"));
-      this->quadrature_1d = quadrature_1d;
-      quadrature          = Quadrature<structdim>(quadrature_1d);
-      n_q_points          = quadrature.size();
-      quadrature_weights.resize(n_q_points);
-      for (unsigned int i = 0; i < n_q_points; ++i)
-        quadrature_weights[i] = quadrature.weight(i);
-
-      for (int d = 0; d < structdim; ++d)
-        {
-          tensor_quadrature_weights[d].resize(quadrature_1d.size());
-          for (unsigned int i = 0; i < quadrature_1d.size(); ++i)
-            tensor_quadrature_weights[d][i] = quadrature_1d.weight(i);
-        }
-
-      // face orientation for faces in 3D
-      if (structdim == spacedim - 1 && spacedim == 3 &&
-          update_flags_inner_faces != update_default)
-        {
-          const unsigned int n = quadrature_1d.size();
-          face_orientations.reinit(8, n * n);
-          for (unsigned int j = 0, i = 0; j < n; ++j)
-            for (unsigned int k = 0; k < n; ++k, ++i)
-              {
-                // face_orientation=true,  face_flip=false, face_rotation=false
-                face_orientations[0][i] = i;
-                // face_orientation=false, face_flip=false, face_rotation=false
-                face_orientations[1][i] = j + k * n;
-                // face_orientation=true,  face_flip=true,  face_rotation=false
-                face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
-                // face_orientation=false, face_flip=true,  face_rotation=false
-                face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
-                // face_orientation=true,  face_flip=false, face_rotation=true
-                face_orientations[4][i] = j + (n - 1 - k) * n;
-                // face_orientation=false, face_flip=false, face_rotation=true
-                face_orientations[5][i] = k + (n - 1 - j) * n;
-                // face_orientation=true,  face_flip=true,  face_rotation=true
-                face_orientations[6][i] = (n - 1 - j) + k * n;
-                // face_orientation=false, face_flip=true,  face_rotation=true
-                face_orientations[7][i] = (n - 1 - k) + j * n;
-              }
-        }
-    }
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    std::size_t
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      QuadratureDescriptor::memory_consumption() const
-    {
-      std::size_t memory = sizeof(this) + quadrature.memory_consumption() +
-                           quadrature_weights.memory_consumption() +
-                           face_orientations.memory_consumption();
-      for (int d = 0; d < structdim; ++d)
-        memory += tensor_quadrature_weights[d].memory_consumption();
-      return memory;
-    }
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    void
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      clear_data_fields()
-    {
-      data_index_offsets.clear();
-      JxW_values.clear();
-      normal_vectors.clear();
-      for (unsigned int i = 0; i < 2; ++i)
-        {
-          jacobians[i].clear();
-          jacobian_gradients[i].clear();
-          normals_times_jacobians[i].clear();
-        }
-      quadrature_point_offsets.clear();
-      quadrature_points.clear();
-    }
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    std::size_t
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      memory_consumption() const
-    {
-      return MemoryConsumption::memory_consumption(descriptor) +
-             MemoryConsumption::memory_consumption(data_index_offsets) +
-             MemoryConsumption::memory_consumption(JxW_values) +
-             MemoryConsumption::memory_consumption(normal_vectors) +
-             MemoryConsumption::memory_consumption(jacobians[0]) +
-             MemoryConsumption::memory_consumption(jacobians[1]) +
-             MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
-             MemoryConsumption::memory_consumption(jacobian_gradients[1]) +
-             MemoryConsumption::memory_consumption(normals_times_jacobians[0]) +
-             MemoryConsumption::memory_consumption(normals_times_jacobians[1]) +
-             MemoryConsumption::memory_consumption(quadrature_point_offsets) +
-             MemoryConsumption::memory_consumption(quadrature_points);
-    }
-
-
-
-    template <int structdim,
-              int spacedim,
-              typename Number,
-              typename VectorizedArrayType>
-    template <typename StreamType>
-    void
-    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
-      print_memory_consumption(StreamType &out, const TaskInfo &task_info) const
-    {
-      // print_memory_statistics involves global communication, so we can
-      // disable the check here only if no processor has any such data
-      const std::size_t size =
-        Utilities::MPI::sum(jacobians[0].size(), task_info.communicator);
-      if (size > 0)
-        {
-          out << "      Memory JxW data:               ";
-          task_info.print_memory_statistics(
-            out,
-            MemoryConsumption::memory_consumption(data_index_offsets) +
-              MemoryConsumption::memory_consumption(JxW_values));
-          out << "      Memory Jacobian data:          ";
-          task_info.print_memory_statistics(
-            out,
-            MemoryConsumption::memory_consumption(jacobians[0]) +
-              MemoryConsumption::memory_consumption(jacobians[1]));
-          out << "      Memory second derivative data: ";
-          task_info.print_memory_statistics(
-            out,
-            MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
-              MemoryConsumption::memory_consumption(jacobian_gradients[1]));
-        }
-      const std::size_t normal_size =
-        Utilities::MPI::sum(normal_vectors.size(), task_info.communicator);
-      if (normal_size > 0)
-        {
-          out << "      Memory normal vectors data:    ";
-          task_info.print_memory_statistics(
-            out,
-            MemoryConsumption::memory_consumption(normal_vectors) +
-              MemoryConsumption::memory_consumption(
-                normals_times_jacobians[0]) +
-              MemoryConsumption::memory_consumption(
-                normals_times_jacobians[1]));
-        }
-
-      const std::size_t quad_size =
-        Utilities::MPI::sum(quadrature_points.size(), task_info.communicator);
-      if (quad_size > 0)
-        {
-          out << "      Memory quadrature points:      ";
-          task_info.print_memory_statistics(
-            out,
-            MemoryConsumption::memory_consumption(quadrature_point_offsets) +
-              MemoryConsumption::memory_consumption(quadrature_points));
-        }
-    }
-
-
-
-    /* ------------------------ MappingInfo implementation ----------------- */
-
     template <int dim, typename Number, typename VectorizedArrayType>
     void
     MappingInfo<dim, Number, VectorizedArrayType>::clear()
diff --git a/include/deal.II/matrix_free/mapping_info_storage.h b/include/deal.II/matrix_free/mapping_info_storage.h
new file mode 100644 (file)
index 0000000..d121e96
--- /dev/null
@@ -0,0 +1,346 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2011 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_mapping_info_storage_h
+#define dealii_matrix_free_mapping_info_storage_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/aligned_vector.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/quadrature.h>
+#include <deal.II/base/tensor.h>
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/fe/fe_update_flags.h>
+
+#include <deal.II/hp/q_collection.h>
+
+#include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    // forward declaration
+    struct TaskInfo;
+
+    /**
+     * An enum to identify various types of cells and faces. The most general
+     * type is what we typically compute in the FEValues context but for many
+     * geometries we can save significant storage.
+     *
+     * @ingroup matrixfree
+     */
+    enum GeometryType : unsigned char
+    {
+      /**
+       * The cell or face is Cartesian.
+       */
+      cartesian = 0,
+
+      /**
+       * The cell or face can be described with an affine mapping.
+       */
+      affine = 1,
+
+      /**
+       * The face is flat, i.e., the normal factor on a face is the same on
+       * all quadrature points. This type is not assigned for cells.
+       */
+      flat_faces = 2,
+
+      /**
+       * There is no special information available for compressing the
+       * representation of the object under consideration.
+       */
+      general = 3
+    };
+
+
+
+    /**
+     * Definition of a structure that stores all cached data related to the
+     * evaluated geometry from the mapping. In order to support hp-adaptivity
+     * and compressed storage (in particular for Jacobians, JxW values, and
+     * normals), storage length can be different for different rows. Thus, it
+     * allows to jump at the data of individual rows similar to compressed row
+     * storage in sparse matrices. We have two different start indices for
+     * fields with different sizes. The first category of offsets are the
+     * indices for Jacobians of the transformation from unit to real cell (we
+     * store the inverse Jacobian), second derivatives, JxW values, and normal
+     * vectors. We keep separate arrays for all these data structures because
+     * a user code might access only some of them. In such a case, one array
+     * will be gone through in a contiguous order with access to all entries,
+     * which makes it easy for the processor to prefetch data. Having all data
+     * in a single array would require some strides in the access pattern,
+     * which is much more complicated for the processor to predict (and indeed
+     * leads to prefetching of data that does not get used on Intel processors
+     * such as BroadwellEP).
+     *
+     * The second category of indices are the offsets for the quadrature
+     * points. Quadrature points can be compressed less than the other fields
+     * and thus need longer fields. Quadrature point indices are often used in
+     * other contexts such as evaluation of right hand sides.
+     *
+     * The third component is a descriptor of data from the unit cells, called
+     * QuadratureDescriptor, which contains the quadrature weights and
+     * permutations of how to go through quadrature points in case of face
+     * data. The latter comes in a vector for the support of hp-adaptivity,
+     * with several data fields for the individual quadrature formulas.
+     *
+     * @ingroup matrixfree
+     */
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    struct MappingInfoStorage
+    {
+      static_assert(
+        std::is_same<Number, typename VectorizedArrayType::value_type>::value,
+        "Type of Number and of VectorizedArrayType do not match.");
+
+      struct QuadratureDescriptor
+      {
+        /**
+         * Constructor. Does nothing.
+         */
+        QuadratureDescriptor();
+
+        /**
+         * Set up the lengths in the various members of this struct.
+         */
+        template <int dim_q>
+        void
+        initialize(const Quadrature<dim_q> &quadrature,
+                   const UpdateFlags update_flags_inner_faces = update_default);
+
+        /**
+         * Set up the lengths in the various members of this struct.
+         */
+        void
+        initialize(const Quadrature<1> &quadrature_1d,
+                   const UpdateFlags update_flags_inner_faces = update_default);
+
+        /**
+         * Returns the memory consumption in bytes.
+         */
+        std::size_t
+        memory_consumption() const;
+
+        /**
+         * Number of quadrature points applied on the given cell or face.
+         */
+        unsigned int n_q_points;
+
+        /**
+         * Original one-dimensional quadrature formula applied on the given
+         * cell or face.
+         */
+        Quadrature<1> quadrature_1d;
+
+        /**
+         * Quadrature formula applied on the given cell or face.
+         */
+        Quadrature<structdim> quadrature;
+
+        /**
+         * Quadrature weights separated by dimension for use in specific
+         * situations.
+         */
+        std::array<AlignedVector<Number>, structdim> tensor_quadrature_weights;
+
+        /**
+         * A cached vector of quadrature weights in the given number format
+         * (non-vectorized, as it is cheap to broadcast the value to all lanes
+         * when it is used in a vectorized context).
+         */
+        AlignedVector<Number> quadrature_weights;
+
+        /**
+         * For quadrature on faces, the evaluation of basis functions is not
+         * in the correct order if a face is not in the standard orientation
+         * to a given element. This data structure is used to re-order the
+         * data evaluated on quadrature points to represent the correct order.
+         */
+        dealii::Table<2, unsigned int> face_orientations;
+      };
+
+      /**
+       * A class describing the layout of the sections in the @p data_storage
+       * field and also includes some data that depends on the number of
+       * quadrature points in the hp-context such as the inner quadrature
+       * formula and re-indexing for faces that are not in the standard
+       * orientation.
+       */
+      std::vector<QuadratureDescriptor> descriptor;
+
+      /**
+       * Collection of quadrature formulae applied on the given face.
+       *
+       * @note Only filled for faces, since faces might be quadrilateral or
+       *   triangle shaped.
+       */
+      std::vector<dealii::hp::QCollection<structdim>> q_collection;
+
+      /**
+       * Stores the index offset into the arrays @p jxw_values, @p jacobians,
+       * @p normal_vectors and the second derivatives. Note that affine cells
+       * have shorter fields of length 1, where the others have lengths equal
+       * to the number of quadrature points of the given cell.
+       */
+      AlignedVector<unsigned int> data_index_offsets;
+
+      /**
+       * The storage of the Jacobian determinant (times the quadrature weight
+       * in case the transformation is non-affine) on quadrature
+       * points.
+       *
+       * Indexed by @p data_index_offsets.
+       */
+      AlignedVector<VectorizedArrayType> JxW_values;
+
+      /**
+       * Stores the normal vectors.
+       *
+       * Indexed by @p data_index_offsets.
+       */
+      AlignedVector<Tensor<1, spacedim, VectorizedArrayType>> normal_vectors;
+
+      /**
+       * The storage of covariant transformation on quadrature points, i.e.,
+       * the inverse and transposed Jacobians of the transformation from the
+       * unit to the real cell.
+       *
+       * Indexed by @p data_index_offsets.
+       *
+       * Contains two fields for access from both sides for interior faces,
+       * but the default case (cell integrals or boundary integrals) only
+       * fills the zeroth component and ignores the first one.
+       */
+      std::array<AlignedVector<Tensor<2, spacedim, VectorizedArrayType>>, 2>
+        jacobians;
+
+      /**
+       * The storage of the gradients of the inverse Jacobian
+       * transformation. Because of symmetry, only the upper diagonal and
+       * diagonal part are needed. The first index runs through the
+       * derivatives, starting with the diagonal and then continuing row-wise,
+       * i.e., $\partial^2/\partial x_1 \partial x_2$ first, then
+       * $\partial^2/\partial x_1 \partial x_3$, and so on. The second index
+       * is the spatial coordinate.
+       *
+       * Indexed by @p data_index_offsets.
+       *
+       * Contains two fields for access from both sides for interior faces,
+       * but the default case (cell integrals or boundary integrals) only
+       * fills the zeroth component and ignores the first one.
+       */
+      std::array<
+        AlignedVector<Tensor<1,
+                             spacedim *(spacedim + 1) / 2,
+                             Tensor<1, spacedim, VectorizedArrayType>>>,
+        2>
+        jacobian_gradients;
+
+      /**
+       * Stores the Jacobian transformations times the normal vector (this
+       * represents a shortcut that is accessed often and can thus get higher
+       * performance).
+       *
+       * Indexed by @p data_index_offsets.
+       */
+      std::array<AlignedVector<Tensor<1, spacedim, VectorizedArrayType>>, 2>
+        normals_times_jacobians;
+
+      /**
+       * Stores the index offset of a particular cell into the quadrature
+       * points array in real coordinates. Note that Cartesian cells have
+       * shorter fields (length is @p n_q_points_1d) than non-Cartesian cells
+       * (length is @p n_q_points) or faces.
+       */
+      AlignedVector<unsigned int> quadrature_point_offsets;
+
+      /**
+       * Stores the quadrature points in real coordinates, including a
+       * compression scheme for Cartesian cells where we do not need to store
+       * the full data on all points.
+       *
+       * Indexed by @p quadrature_point_offsets.
+       */
+      AlignedVector<Point<spacedim, VectorizedArrayType>> quadrature_points;
+
+      /**
+       * Clears all data fields except the descriptor vector.
+       */
+      void
+      clear_data_fields();
+
+      /**
+       * Returns the quadrature index for a given number of quadrature
+       * points. If not in hp-mode or if the index is not found, this
+       * function always returns index 0. Hence, this function does not
+       * check whether the given degree is actually present.
+       */
+      unsigned int
+      quad_index_from_n_q_points(const unsigned int n_q_points) const;
+
+      /**
+       * Prints a detailed summary of memory consumption in the different
+       * structures of this class to the given output stream.
+       */
+      template <typename StreamType>
+      void
+      print_memory_consumption(StreamType &    out,
+                               const TaskInfo &task_info) const;
+
+      /**
+       * Returns the memory consumption in bytes.
+       */
+      std::size_t
+      memory_consumption() const;
+    };
+
+
+
+    /* ------------------- inline functions ----------------------------- */
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    inline unsigned int
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      quad_index_from_n_q_points(const unsigned int n_q_points) const
+    {
+      for (unsigned int i = 0; i < descriptor.size(); ++i)
+        if (n_q_points == descriptor[i].n_q_points)
+          return i;
+      return 0;
+    }
+
+  } // end of namespace MatrixFreeFunctions
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/matrix_free/mapping_info_storage.templates.h b/include/deal.II/matrix_free/mapping_info_storage.templates.h
new file mode 100644 (file)
index 0000000..69509d1
--- /dev/null
@@ -0,0 +1,267 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2011 - 2021 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_matrix_free_mapping_info_storage_templates_h
+#define dealii_matrix_free_mapping_info_storage_templates_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/memory_consumption.h>
+#include <deal.II/base/multithread_info.h>
+#include <deal.II/base/thread_management.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_nothing.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/matrix_free/evaluation_template_factory.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
+#include <deal.II/matrix_free/task_info.h>
+#include <deal.II/matrix_free/util.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace MatrixFreeFunctions
+  {
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      QuadratureDescriptor::QuadratureDescriptor()
+      : n_q_points(numbers::invalid_unsigned_int)
+    {}
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    template <int dim_q>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      QuadratureDescriptor::initialize(
+        const Quadrature<dim_q> &quadrature,
+        const UpdateFlags        update_flags_inner_faces)
+    {
+      Assert(structdim + 1 <= spacedim ||
+               update_flags_inner_faces == update_default,
+             ExcMessage("Volume cells do not allow for setting inner faces"));
+      this->quadrature = quadrature;
+      n_q_points       = quadrature.size();
+      quadrature_weights.resize(n_q_points);
+      for (unsigned int i = 0; i < n_q_points; ++i)
+        quadrature_weights[i] = quadrature.weight(i);
+
+      // note: quadrature_1d and tensor_quadrature_weights are not set up
+
+      // TODO: set up face_orientations
+      (void)update_flags_inner_faces;
+    }
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      QuadratureDescriptor::initialize(
+        const Quadrature<1> &quadrature_1d,
+        const UpdateFlags    update_flags_inner_faces)
+    {
+      Assert(structdim + 1 <= spacedim ||
+               update_flags_inner_faces == update_default,
+             ExcMessage("Volume cells do not allow for setting inner faces"));
+      this->quadrature_1d = quadrature_1d;
+      quadrature          = Quadrature<structdim>(quadrature_1d);
+      n_q_points          = quadrature.size();
+      quadrature_weights.resize(n_q_points);
+      for (unsigned int i = 0; i < n_q_points; ++i)
+        quadrature_weights[i] = quadrature.weight(i);
+
+      for (int d = 0; d < structdim; ++d)
+        {
+          tensor_quadrature_weights[d].resize(quadrature_1d.size());
+          for (unsigned int i = 0; i < quadrature_1d.size(); ++i)
+            tensor_quadrature_weights[d][i] = quadrature_1d.weight(i);
+        }
+
+      // face orientation for faces in 3D
+      if (structdim == spacedim - 1 && spacedim == 3 &&
+          update_flags_inner_faces != update_default)
+        {
+          const unsigned int n = quadrature_1d.size();
+          face_orientations.reinit(8, n * n);
+          for (unsigned int j = 0, i = 0; j < n; ++j)
+            for (unsigned int k = 0; k < n; ++k, ++i)
+              {
+                // face_orientation=true,  face_flip=false, face_rotation=false
+                face_orientations[0][i] = i;
+                // face_orientation=false, face_flip=false, face_rotation=false
+                face_orientations[1][i] = j + k * n;
+                // face_orientation=true,  face_flip=true,  face_rotation=false
+                face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
+                // face_orientation=false, face_flip=true,  face_rotation=false
+                face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
+                // face_orientation=true,  face_flip=false, face_rotation=true
+                face_orientations[4][i] = j + (n - 1 - k) * n;
+                // face_orientation=false, face_flip=false, face_rotation=true
+                face_orientations[5][i] = k + (n - 1 - j) * n;
+                // face_orientation=true,  face_flip=true,  face_rotation=true
+                face_orientations[6][i] = (n - 1 - j) + k * n;
+                // face_orientation=false, face_flip=true,  face_rotation=true
+                face_orientations[7][i] = (n - 1 - k) + j * n;
+              }
+        }
+    }
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    std::size_t
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      QuadratureDescriptor::memory_consumption() const
+    {
+      std::size_t memory = sizeof(this) + quadrature.memory_consumption() +
+                           quadrature_weights.memory_consumption() +
+                           face_orientations.memory_consumption();
+      for (int d = 0; d < structdim; ++d)
+        memory += tensor_quadrature_weights[d].memory_consumption();
+      return memory;
+    }
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      clear_data_fields()
+    {
+      data_index_offsets.clear();
+      JxW_values.clear();
+      normal_vectors.clear();
+      for (unsigned int i = 0; i < 2; ++i)
+        {
+          jacobians[i].clear();
+          jacobian_gradients[i].clear();
+          normals_times_jacobians[i].clear();
+        }
+      quadrature_point_offsets.clear();
+      quadrature_points.clear();
+    }
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    std::size_t
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      memory_consumption() const
+    {
+      return MemoryConsumption::memory_consumption(descriptor) +
+             MemoryConsumption::memory_consumption(data_index_offsets) +
+             MemoryConsumption::memory_consumption(JxW_values) +
+             MemoryConsumption::memory_consumption(normal_vectors) +
+             MemoryConsumption::memory_consumption(jacobians[0]) +
+             MemoryConsumption::memory_consumption(jacobians[1]) +
+             MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
+             MemoryConsumption::memory_consumption(jacobian_gradients[1]) +
+             MemoryConsumption::memory_consumption(normals_times_jacobians[0]) +
+             MemoryConsumption::memory_consumption(normals_times_jacobians[1]) +
+             MemoryConsumption::memory_consumption(quadrature_point_offsets) +
+             MemoryConsumption::memory_consumption(quadrature_points);
+    }
+
+
+
+    template <int structdim,
+              int spacedim,
+              typename Number,
+              typename VectorizedArrayType>
+    template <typename StreamType>
+    void
+    MappingInfoStorage<structdim, spacedim, Number, VectorizedArrayType>::
+      print_memory_consumption(StreamType &out, const TaskInfo &task_info) const
+    {
+      // print_memory_statistics involves global communication, so we can
+      // disable the check here only if no processor has any such data
+      const std::size_t size =
+        Utilities::MPI::sum(jacobians[0].size(), task_info.communicator);
+      if (size > 0)
+        {
+          out << "      Memory JxW data:               ";
+          task_info.print_memory_statistics(
+            out,
+            MemoryConsumption::memory_consumption(data_index_offsets) +
+              MemoryConsumption::memory_consumption(JxW_values));
+          out << "      Memory Jacobian data:          ";
+          task_info.print_memory_statistics(
+            out,
+            MemoryConsumption::memory_consumption(jacobians[0]) +
+              MemoryConsumption::memory_consumption(jacobians[1]));
+          out << "      Memory second derivative data: ";
+          task_info.print_memory_statistics(
+            out,
+            MemoryConsumption::memory_consumption(jacobian_gradients[0]) +
+              MemoryConsumption::memory_consumption(jacobian_gradients[1]));
+        }
+      const std::size_t normal_size =
+        Utilities::MPI::sum(normal_vectors.size(), task_info.communicator);
+      if (normal_size > 0)
+        {
+          out << "      Memory normal vectors data:    ";
+          task_info.print_memory_statistics(
+            out,
+            MemoryConsumption::memory_consumption(normal_vectors) +
+              MemoryConsumption::memory_consumption(
+                normals_times_jacobians[0]) +
+              MemoryConsumption::memory_consumption(
+                normals_times_jacobians[1]));
+        }
+
+      const std::size_t quad_size =
+        Utilities::MPI::sum(quadrature_points.size(), task_info.communicator);
+      if (quad_size > 0)
+        {
+          out << "      Memory quadrature points:      ";
+          task_info.print_memory_statistics(
+            out,
+            MemoryConsumption::memory_consumption(quadrature_point_offsets) +
+              MemoryConsumption::memory_consumption(quadrature_points));
+        }
+    }
+
+  } // namespace MatrixFreeFunctions
+} // end of namespace internal
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

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.