]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Address review comments
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Thu, 18 Jan 2024 21:18:29 +0000 (22:18 +0100)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 23 Jan 2024 09:33:33 +0000 (10:33 +0100)
include/deal.II/matrix_free/evaluation_template_face_factory.templates.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/non_matching/mapping_info.h
tests/non_matching/mapping_info_04.cc
tests/non_matching/mapping_info_05.cc

index 79519e39a84fd5be78c51a644291c50dac3f1df0..608d03e1b6a48a0a8511c161e9cc51385f0ba291 100644 (file)
@@ -47,6 +47,7 @@ namespace internal
   }
 
 
+
   template <int dim, typename Number>
   void
   FEFaceEvaluationFactory<dim, Number>::project_to_face(
index 1ed360ce38877ec5718033ec207d676a799e450e..dfe94b4dd46c0d8104d38b047673d89ce80bb889 100644 (file)
@@ -2750,20 +2750,18 @@ public:
   boundary_id() const;
 
   /**
-   * The number of degrees of freedom of a single component on the cell for
-   * the underlying evaluation object. Usually close to
-   * static_dofs_per_component, but the number depends on the actual element
-   * selected and is thus not static.
+   * Get the number of degrees of freedom of a single component which are
+   * projected onto a face.
    */
-  const unsigned int dofs_per_component_on_cell;
+  unsigned int
+  get_dofs_per_component_projected_to_face();
 
   /**
-   * The number of degrees of freedom on the cell accumulated over all
-   * components in the current evaluation object. Usually close to
-   * static_dofs_per_cell = static_dofs_per_component*n_components, but the
-   * number depends on the actual element selected and is thus not static.
+   * Get the number of degrees of freedom accumulated over all
+   * components which are projected onto a face.
    */
-  const unsigned int dofs_per_cell;
+  unsigned int
+  get_dofs_projected_to_face();
 
   /**
    * The number of degrees of freedom of a single component on the cell for
@@ -2771,7 +2769,7 @@ public:
    * static_dofs_per_component, but the number depends on the actual element
    * selected and is thus not static.
    */
-  const unsigned int dofs_per_component_on_face;
+  const unsigned int dofs_per_component;
 
   /**
    * The number of degrees of freedom on the cell accumulated over all
@@ -2779,7 +2777,7 @@ public:
    * static_dofs_per_cell = static_dofs_per_component*n_components, but the
    * number depends on the actual element selected and is thus not static.
    */
-  const unsigned int dofs_per_face;
+  const unsigned int dofs_per_cell;
 
   /**
    * The number of quadrature points in use. If the number of quadrature
@@ -8310,10 +8308,8 @@ inline FEFaceEvaluation<dim,
               active_fe_index,
               active_quad_index,
               face_type)
-  , dofs_per_component_on_cell(this->data->dofs_per_component_on_cell)
+  , dofs_per_component(this->data->dofs_per_component_on_cell)
   , dofs_per_cell(this->data->dofs_per_component_on_cell * n_components_)
-  , dofs_per_component_on_face(this->data->dofs_per_component_on_face)
-  , dofs_per_face(this->data->dofs_per_component_on_face * n_components_)
   , n_q_points(this->n_quadrature_points)
 {}
 
@@ -9414,6 +9410,45 @@ FEFaceEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+unsigned int
+FEFaceEvaluation<
+  dim,
+  fe_degree,
+  n_q_points_1d,
+  n_components_,
+  Number,
+  VectorizedArrayType>::get_dofs_per_component_projected_to_face()
+{
+  return this->data->dofs_per_component_on_face;
+}
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+unsigned int
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::get_dofs_projected_to_face()
+{
+  return this->data->dofs_per_component_on_face * n_components_;
+}
+
+
+
 /*------------------------- end FEFaceEvaluation ------------------------- */
 
 
index da845ef1bd97b2746c941bd277a74a47d3317923..8ff84419aa62762e5cb5dae63e9a6de9e1b0dbaf 100644 (file)
@@ -1111,1004 +1111,1309 @@ protected:
   const bool is_interior;
 };
 
-// ----------------------- template and inline function ----------------------
 
 
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(const Mapping<dim>       &mapping,
-                        const FiniteElement<dim> &fe,
-                        const UpdateFlags         update_flags,
-                        const unsigned int        first_selected_component)
-  : n_q_batches(numbers::invalid_unsigned_int)
-  , n_q_points(numbers::invalid_unsigned_int)
-  , n_q_points_scalar(numbers::invalid_unsigned_int)
-  , mapping(&mapping)
-  , fe(&fe)
-  , JxW_ptr(nullptr)
-  , update_flags(update_flags)
-  , mapping_info_on_the_fly(
-      std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
-        mapping,
-        update_flags))
-  , mapping_info(mapping_info_on_the_fly.get())
-  , current_cell_index(numbers::invalid_unsigned_int)
-  , current_face_number(numbers::invalid_unsigned_int)
-  , is_reinitialized(false)
-  , is_interior(true)
+/**
+ * This class provides an interface to the evaluation of interpolated solution
+ * values and gradients on cells on arbitrary reference point positions. These
+ * points can change from cell to cell, both with respect to their quantity as
+ * well to the location. The two typical use cases are evaluations on
+ * non-matching grids and particle simulations.
+ *
+ * The use of this class is similar to FEValues or FEEvaluation: The class is
+ * first initialized to a cell by calling `FEPointEvaluation::reinit(cell,
+ * unit_points)`, with the main difference to the other concepts that the
+ * underlying points in reference coordinates need to be passed along. Then,
+ * upon call to evaluate() or integrate(), the user can compute information at
+ * the give points. Eventually, the access functions get_value() or
+ * get_gradient() allow to query this information at a specific point index.
+ *
+ * The functionality is similar to creating an FEValues object with a
+ * Quadrature object on the `unit_points` on every cell separately and then
+ * calling FEValues::get_function_values or FEValues::get_function_gradients,
+ * and for some elements and mappings this is what actually happens
+ * internally. For specific combinations of Mapping and FiniteElement
+ * realizations, however, there is a much more efficient implementation that
+ * avoids the memory allocation and other expensive start-up cost of
+ * FEValues. Currently, the functionality is specialized for mappings derived
+ * from MappingQ and MappingCartesian and for finite elements with tensor
+ * product structure that work with the
+ * @ref matrixfree
+ * module. In those cases, the cost implied
+ * by this class is similar (or sometimes even somewhat lower) than using
+ * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`.
+ */
+template <int n_components_,
+          int dim,
+          int spacedim    = dim,
+          typename Number = double>
+class FEPointEvaluation
+  : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
 {
-  setup(first_selected_component);
-}
+public:
+  static constexpr unsigned int dimension    = dim;
+  static constexpr unsigned int n_components = n_components_;
 
+  using number_type = Number;
 
+  using ScalarNumber =
+    typename internal::VectorizedArrayTrait<Number>::value_type;
+  using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+    Number>::vectorized_value_type;
+  using ETT = typename internal::FEPointEvaluation::
+    EvaluatorTypeTraits<dim, n_components, Number>;
+  using value_type            = typename ETT::value_type;
+  using scalar_value_type     = typename ETT::scalar_value_type;
+  using vectorized_value_type = typename ETT::vectorized_value_type;
+  using gradient_type         = typename ETT::gradient_type;
+  using interface_vectorized_gradient_type =
+    typename ETT::interface_vectorized_gradient_type;
 
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(
+  /**
+   * Constructor.
+   *
+   * @param mapping The Mapping class describing the actual geometry of a cell
+   * passed to the evaluate() function.
+   *
+   * @param fe The FiniteElement object that is used for the evaluation, which
+   * is typically the same on all cells to be evaluated.
+   *
+   * @param update_flags Specify the quantities to be computed by the mapping
+   * during the call of reinit(). During evaluate() or integrate(), this data
+   * is queried to produce the desired result (e.g., the gradient of a finite
+   * element solution).
+   *
+   * @param first_selected_component For multi-component FiniteElement
+   * objects, this parameter allows to select a range of `n_components`
+   * components starting from this parameter.
+   */
+  FEPointEvaluation(const Mapping<dim>       &mapping,
+                    const FiniteElement<dim> &fe,
+                    const UpdateFlags         update_flags,
+                    const unsigned int        first_selected_component = 0);
+
+  /**
+   * Constructor to make the present class able to re-use the geometry
+   * data also used by other `FEPointEvaluation` objects.
+   *
+   * @param mapping_info The MappingInfo class describes the geometry-related
+   * data for evaluating finite-element solutions. This object enables to
+   * construct such an object on the outside, possibly re-using it between
+   * several objects or between several calls to the same cell and unit points.
+   *
+   * @param fe The FiniteElement object that is used for the evaluation, which
+   * is typically the same on all cells to be evaluated.
+   *
+   * @param first_selected_component For multi-component FiniteElement
+   * objects, this parameter allows to select a range of `n_components`
+   * components starting from this parameter.
+   */
+  FEPointEvaluation(
     NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
     const FiniteElement<dim>                        &fe,
-    const unsigned int                               first_selected_component,
-    const bool                                       is_interior)
-  : n_q_batches(numbers::invalid_unsigned_int)
-  , n_q_points(numbers::invalid_unsigned_int)
-  , n_q_points_scalar(numbers::invalid_unsigned_int)
-  , mapping(&mapping_info.get_mapping())
-  , fe(&fe)
-  , JxW_ptr(nullptr)
-  , update_flags(mapping_info.get_update_flags())
-  , mapping_info(&mapping_info)
-  , current_cell_index(numbers::invalid_unsigned_int)
-  , current_face_number(numbers::invalid_unsigned_int)
-  , is_reinitialized(false)
-  , is_interior(is_interior)
-{
-  setup(first_selected_component);
-  connection_is_reinitialized = mapping_info.connect_is_reinitialized(
-    [this]() { this->is_reinitialized = false; });
-}
+    const unsigned int first_selected_component = 0);
 
+  /**
+   * Set up the mapping information for the given cell, e.g., by computing the
+   * Jacobian of the mapping for the given points if gradients of the functions
+   * are requested.
+   *
+   * @param[in] cell An iterator to the current cell
+   *
+   * @param[in] unit_points List of points in the reference locations of the
+   * current cell where the FiniteElement object should be
+   * evaluated/integrated in the evaluate() and integrate() functions.
+   */
+  void
+  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+         const ArrayView<const Point<dim>> &unit_points);
 
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the single cell in the MappingInfo object.
+   */
+  void
+  reinit();
 
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(
-    FEPointEvaluationBase<n_components_, dim, spacedim, Number> &other) noexcept
-  : n_q_batches(other.n_q_batches)
-  , n_q_points(other.n_q_points)
-  , n_q_points_scalar(other.n_q_points_scalar)
-  , mapping(other.mapping)
-  , fe(other.fe)
-  , poly(other.poly)
-  , use_linear_path(other.use_linear_path)
-  , renumber(other.renumber)
-  , solution_renumbered(other.solution_renumbered)
-  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
-  , values(other.values)
-  , gradients(other.gradients)
-  , dofs_per_component(other.dofs_per_component)
-  , dofs_per_component_face(other.dofs_per_component_face)
-  , component_in_base_element(other.component_in_base_element)
-  , nonzero_shape_function_component(other.nonzero_shape_function_component)
-  , update_flags(other.update_flags)
-  , fe_values(other.fe_values)
-  , mapping_info_on_the_fly(
-      other.mapping_info_on_the_fly ?
-        std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
-          *mapping,
-          update_flags) :
-        nullptr)
-  , mapping_info(other.mapping_info)
-  , current_cell_index(other.current_cell_index)
-  , current_face_number(other.current_face_number)
-  , fast_path(other.fast_path)
-  , is_reinitialized(false)
-  , shapes(other.shapes)
-  , shapes_faces(other.shapes_faces)
-  , is_interior(other.is_interior)
-{
-  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
-    [this]() { this->is_reinitialized = false; });
-}
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the cell in the MappingInfo object.
+   */
+  void
+  reinit(const unsigned int cell_index);
 
 
+  /**
+   * This function interpolates the finite element solution, represented by
+   * `solution_values`, on the cell and `unit_points` passed to reinit().
+   *
+   * @param[in] solution_values This array is supposed to contain the unknown
+   * values on the element read out by
+   * `FEEvaluation::read_dof_values(global_vector)`.
+   *
+   * @param[in] evaluation_flags Flags specifying which quantities should be
+   * evaluated at the points.
+   */
+  template <std::size_t stride_view>
+  void
+  evaluate(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  FEPointEvaluationBase(
-    FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-      &&other) noexcept
-  : n_q_batches(other.n_q_batches)
-  , n_q_points(other.n_q_points)
-  , n_q_points_scalar(other.n_q_points_scalar)
-  , mapping(other.mapping)
-  , fe(other.fe)
-  , poly(other.poly)
-  , use_linear_path(other.use_linear_path)
-  , renumber(other.renumber)
-  , solution_renumbered(other.solution_renumbered)
-  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
-  , values(other.values)
-  , gradients(other.gradients)
-  , dofs_per_component(other.dofs_per_component)
-  , dofs_per_component_face(other.dofs_per_component_face)
-  , component_in_base_element(other.component_in_base_element)
-  , nonzero_shape_function_component(other.nonzero_shape_function_component)
-  , update_flags(other.update_flags)
-  , fe_values(other.fe_values)
-  , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
-  , mapping_info(other.mapping_info)
-  , current_cell_index(other.current_cell_index)
-  , current_face_number(other.current_face_number)
-  , fast_path(other.fast_path)
-  , is_reinitialized(false)
-  , shapes(other.shapes)
-  , shapes_faces(other.shapes_faces)
-  , is_interior(other.is_interior)
-{
-  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
-    [this]() { this->is_reinitialized = false; });
-}
-
-
+  /**
+   * This function interpolates the finite element solution, represented by
+   * `solution_values`, on the cell and `unit_points` passed to reinit().
+   *
+   * @param[in] solution_values This array is supposed to contain the unknown
+   * values on the element as returned by `cell->get_dof_values(global_vector,
+   * solution_values)`.
+   *
+   * @param[in] evaluation_flags Flags specifying which quantities should be
+   * evaluated at the points.
+   */
+  void
+  evaluate(const ArrayView<const ScalarNumber>    &solution_values,
+           const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
-template <int n_components_, int dim, int spacedim, typename Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  ~FEPointEvaluationBase()
-{
-  connection_is_reinitialized.disconnect();
-}
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used during
+   * `FEEvaluation::set_dof_values(global_vector)` or
+   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used to during
+   * `cell->set_dof_values(solution_values, global_vector)` or
+   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  void
+  integrate(const ArrayView<ScalarNumber>          &solution_values,
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points. This is
+   * similar to the integration of a bilinear form in terms of the test
+   * function, with the difference that this formula does not include a `JxW`
+   * factor (in contrast to the integrate function of this class). This allows
+   * the class to naturally embed point information (e.g. particles) into a
+   * finite element formulation.
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used during
+   * `FEEvaluation::set_dof_values(global_vector)` or
+   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  test_and_sum(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values = false);
 
-template <int n_components_, int dim, int spacedim, typename Number>
-void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::setup(
-  const unsigned int first_selected_component)
-{
-  AssertIndexRange(first_selected_component + n_components,
-                   fe->n_components() + 1);
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points. This is
+   * similar to the integration of a bilinear form in terms of the test
+   * function, with the difference that this formula does not include a `JxW`
+   * factor (in contrast to the integrate function of this class). This allows
+   * the class to naturally embed point information (e.g. particles) into a
+   * finite element formulation.
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used during
+   * `cell->set_dof_values(solution_values, global_vector)` or
+   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  void
+  test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
+               const EvaluationFlags::EvaluationFlags &integration_flags,
+               const bool                              sum_into_values = false);
 
-  shapes.reserve(100);
+  /**
+   * Return the normal vector. This class or the MappingInfo object passed to
+   * this function needs to be constructed with UpdateFlags containing
+   * `update_normal_vectors`.
+   */
+  Tensor<1, spacedim, Number>
+  normal_vector(const unsigned int point_index) const;
 
-  bool         same_base_element   = true;
-  unsigned int base_element_number = 0;
-  component_in_base_element        = 0;
-  unsigned int component           = 0;
-  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
-    if (component + fe->element_multiplicity(base_element_number) >
-        first_selected_component)
-      {
-        if (first_selected_component + n_components >
-            component + fe->element_multiplicity(base_element_number))
-          same_base_element = false;
-        component_in_base_element = first_selected_component - component;
-        break;
-      }
-    else
-      component += fe->element_multiplicity(base_element_number);
+private:
+  static constexpr std::size_t n_lanes_user_interface =
+    internal::VectorizedArrayTrait<Number>::width();
+  static constexpr std::size_t n_lanes_internal =
+    internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+  static constexpr std::size_t stride =
+    internal::VectorizedArrayTrait<Number>::stride();
 
-  if (internal::FEPointEvaluation::is_fast_path_supported(*mapping) &&
-      internal::FEPointEvaluation::is_fast_path_supported(
-        *fe, base_element_number) &&
-      same_base_element)
-    {
-      shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
-      renumber                = shape_info.lexicographic_numbering;
-      dofs_per_component      = shape_info.dofs_per_component_on_cell;
-      dofs_per_component_face = shape_info.dofs_per_component_on_face;
-      poly = internal::FEPointEvaluation::get_polynomial_space(
-        fe->base_element(base_element_number));
+  /**
+   * Resizes necessary data fields, reads in and renumbers solution values.
+   * Interpolates onto face if face path is selected.
+   */
+  template <bool is_linear, std::size_t stride_view>
+  void
+  prepare_evaluate_fast(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values);
 
-      bool is_lexicographic = true;
-      for (unsigned int i = 0; i < renumber.size(); ++i)
-        if (i != renumber[i])
-          is_lexicographic = false;
+  /**
+   * Evaluates the actual interpolation on the cell or face for a quadrature
+   * batch.
+   */
+  template <bool is_linear, std::size_t stride_view>
+  void
+  compute_evaluate_fast(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags,
+    const unsigned int                                       n_shapes,
+    const unsigned int                                       qb,
+    vectorized_value_type                                   &value,
+    interface_vectorized_gradient_type                      &gradient);
 
-      if (is_lexicographic)
-        renumber.clear();
+  /**
+   * Fast path of the evaluate function.
+   */
+  template <bool is_linear, std::size_t stride_view>
+  void
+  evaluate_fast(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
-      use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
-                         poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
-                         poly[1].value(1.) == 1.) &&
-                        (fe->n_components() == n_components);
+  /**
+   * Slow path of the evaluate function using FEValues.
+   */
+  template <std::size_t stride_view>
+  void
+  evaluate_slow(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
-      const unsigned int size_face = 3 * dofs_per_component_face * n_components;
-      const unsigned int size_cell = dofs_per_component * n_components;
-      scratch_data_scalar.resize(size_face + size_cell);
+  /**
+   * Integrates the product of the data passed in by submit_value() and
+   * submit_gradient() with the values or gradients of test functions on the
+   * cell or face for a given quadrature batch.
+   */
+  template <bool is_linear>
+  void
+  compute_integrate_fast(
+    const EvaluationFlags::EvaluationFlags  &integration_flags,
+    const unsigned int                       n_shapes,
+    const unsigned int                       qb,
+    const vectorized_value_type              value,
+    const interface_vectorized_gradient_type gradient,
+    vectorized_value_type                   *solution_values_vectorized_linear);
 
-      solution_renumbered.resize(dofs_per_component);
-      solution_renumbered_vectorized.resize(dofs_per_component);
+  /**
+   * Addition across the lanes of VectorizedArray as accumulated by the
+   * compute_integrate_fast_function(), writing the sum into the result vector.
+   * Applies face contributions to cell contributions for face path.
+   */
+  template <bool is_linear, std::size_t stride_view>
+  void
+  finish_integrate_fast(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    vectorized_value_type *solution_values_vectorized_linear,
+    const bool             sum_into_values);
 
-      fast_path = true;
-    }
-  else
-    {
-      nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
-      for (unsigned int d = 0; d < n_components; ++d)
-        {
-          const unsigned int component = first_selected_component + d;
-          for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
-            {
-              const bool is_primitive =
-                fe->is_primitive() || fe->is_primitive(i);
-              if (is_primitive)
-                nonzero_shape_function_component[i][d] =
-                  (component == fe->system_to_component_index(i).first);
-              else
-                nonzero_shape_function_component[i][d] =
-                  (fe->get_nonzero_components(i)[component] == true);
-            }
-        }
+  /**
+   * Fast path of the integrate function.
+   */
+  template <bool do_JxW, bool is_linear, std::size_t stride_view>
+  void
+  integrate_fast(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
-      fast_path = false;
-    }
-}
+  /**
+   * Slow path of the integrate function using FEValues.
+   */
+  template <bool do_JxW, std::size_t stride_view>
+  void
+  integrate_slow(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
+  /**
+   * Implementation of the integrate/test_and_sum function.
+   */
+  template <bool do_JxW, std::size_t stride_view>
+  void
+  do_integrate(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
+};
 
 
-template <int n_components_, int dim, int spacedim, typename Number>
-template <bool is_face, bool is_linear>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::do_reinit()
-{
-  const unsigned int geometry_index =
-    mapping_info->template compute_geometry_index_offset<is_face>(
-      current_cell_index, current_face_number);
 
-  cell_type = mapping_info->get_cell_type(geometry_index);
+/**
+ * This class provides an interface to the evaluation of interpolated solution
+ * values and gradients on faces on arbitrary reference point positions. These
+ * points can change from face to face, both with respect to their quantity as
+ * well to the location. A typical use case is evaluations on non-matching
+ * grids.
+ *
+ * The use of this class is similar to FEEvaluation: In the constructor, a
+ * reference to a NonMatching::MappingInfo object is passed, where the
+ * quadrature points in reference position is stored together with the mapping
+ * information. The class is then reinitialized to a cell by calling
+ * `FEFacePointEvaluation::reinit(face_index)` or
+ * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to
+ * evaluate() or integrate(), the user can compute information at the given
+ * points. Eventually, the access functions get_value() or get_gradient() allow
+ * to query this information at a specific point index.
+ */
+template <int n_components_,
+          int dim,
+          int spacedim    = dim,
+          typename Number = double>
+class FEFacePointEvaluation
+  : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+{
+public:
+  static constexpr unsigned int dimension    = dim;
+  static constexpr unsigned int n_components = n_components_;
 
-  const_cast<unsigned int &>(n_q_points_scalar) =
-    mapping_info->get_n_q_points_unvectorized(geometry_index);
+  using number_type = Number;
 
-  // round up n_q_points_scalar / n_lanes_internal
-  const_cast<unsigned int &>(n_q_batches) =
-    (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
+  using ScalarNumber =
+    typename internal::VectorizedArrayTrait<Number>::value_type;
+  using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+    Number>::vectorized_value_type;
+  using ETT = typename internal::FEPointEvaluation::
+    EvaluatorTypeTraits<dim, n_components, Number>;
+  using value_type            = typename ETT::value_type;
+  using scalar_value_type     = typename ETT::scalar_value_type;
+  using vectorized_value_type = typename ETT::vectorized_value_type;
+  using gradient_type         = typename ETT::gradient_type;
+  using interface_vectorized_gradient_type =
+    typename ETT::interface_vectorized_gradient_type;
 
-  const unsigned int n_q_points_before = n_q_points;
+  /**
+   * Constructor. Allows to select if interior or exterior face is selected.
+   */
+  FEFacePointEvaluation(
+    NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+    const FiniteElement<dim>                        &fe,
+    const bool                                       is_interior = true,
+    const unsigned int first_selected_component                  = 0);
 
-  const_cast<unsigned int &>(n_q_points) =
-    (stride == 1) ? n_q_batches : n_q_points_scalar;
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the face in the MappingInfo object. Used in element-centric loops (ECL).
+   */
+  void
+  reinit(const unsigned int cell_index, const unsigned int face_number);
 
-  if (n_q_points != n_q_points_before)
-    {
-      if (update_flags & update_values)
-        values.resize(n_q_points);
-      if (update_flags & update_gradients)
-        gradients.resize(n_q_points);
-    }
+  /**
+   * Reinitialize the evaluator to point to the correct precomputed mapping of
+   * the face in the MappingInfo object. Used in face-centric loops (FCL).
+   */
+  void
+  reinit(const unsigned int face_index);
 
-  if (n_q_points == 0)
-    {
-      is_reinitialized = true;
-      return;
-    }
+  /**
+   * This function interpolates the finite element solution, represented by
+   * `solution_values`, on the cell and `unit_points` passed to reinit().
+   *
+   * @param[in] solution_values This array is supposed to contain the unknown
+   * values on the element read out by
+   * `FEEvaluation::read_dof_values(global_vector)`.
+   *
+   * @param[in] evaluation_flags Flags specifying which quantities should be
+   * evaluated at the points.
+   */
+  template <std::size_t stride_view>
+  void
+  evaluate(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
-  // set unit point pointer
-  const unsigned int unit_point_offset =
-    mapping_info->compute_unit_point_index_offset(geometry_index);
+  /**
+   * This function interpolates the finite element solution, represented by
+   * `solution_values`, on the cell and `unit_points` passed to reinit().
+   *
+   * @param[in] solution_values This array is supposed to contain the unknown
+   * values on the element as returned by `cell->get_dof_values(global_vector,
+   * solution_values)`.
+   *
+   * @param[in] evaluation_flags Flags specifying which quantities should be
+   * evaluated at the points.
+   */
+  void
+  evaluate(const ArrayView<const ScalarNumber>    &solution_values,
+           const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
-  if (is_face)
-    unit_point_faces_ptr =
-      mapping_info->get_unit_point_faces(unit_point_offset);
-  else
-    unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used during
+   * `FEEvaluation::set_dof_values(global_vector)` or
+   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
-  // set data pointers
-  const unsigned int data_offset =
-    mapping_info->compute_data_index_offset(geometry_index);
-  const unsigned int compressed_data_offset =
-    mapping_info->compute_compressed_data_index_offset(geometry_index);
-#ifdef DEBUG
-  const UpdateFlags update_flags_mapping =
-    mapping_info->get_update_flags_mapping();
-  if (update_flags_mapping & UpdateFlags::update_quadrature_points)
-    real_point_ptr = mapping_info->get_real_point(data_offset);
-  if (update_flags_mapping & UpdateFlags::update_jacobians)
-    jacobian_ptr =
-      mapping_info->get_jacobian(compressed_data_offset, is_interior);
-  if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
-    inverse_jacobian_ptr =
-      mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
-  if (update_flags_mapping & UpdateFlags::update_normal_vectors)
-    normal_ptr = mapping_info->get_normal_vector(data_offset);
-  if (update_flags_mapping & UpdateFlags::update_JxW_values)
-    JxW_ptr = mapping_info->get_JxW(data_offset);
-#else
-  real_point_ptr = mapping_info->get_real_point(data_offset);
-  jacobian_ptr =
-    mapping_info->get_jacobian(compressed_data_offset, is_interior);
-  inverse_jacobian_ptr =
-    mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
-  normal_ptr = mapping_info->get_normal_vector(data_offset);
-  JxW_ptr    = mapping_info->get_JxW(data_offset);
-#endif
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used to during
+   * `cell->set_dof_values(solution_values, global_vector)` or
+   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  void
+  integrate(const ArrayView<ScalarNumber>          &solution_values,
+            const EvaluationFlags::EvaluationFlags &integration_flags,
+            const bool                              sum_into_values = false);
 
-  if (!is_linear && fast_path)
-    {
-      const std::size_t n_shapes = poly.size();
-
-      for (unsigned int qb = 0; qb < n_q_batches; ++qb)
-        if (is_face)
-          {
-            if (dim > 1)
-              {
-                shapes_faces.resize_fast(n_q_batches * n_shapes);
-                internal::compute_values_of_array(
-                  shapes_faces.data() + qb * n_shapes,
-                  poly,
-                  unit_point_faces_ptr[qb],
-                  update_flags & UpdateFlags::update_gradients ? 1 : 0);
-              }
-          }
-        else
-          {
-            shapes.resize_fast(n_q_batches * n_shapes);
-            internal::compute_values_of_array(
-              shapes.data() + qb * n_shapes,
-              poly,
-              unit_point_ptr[qb],
-              update_flags & UpdateFlags::update_gradients ? 1 : 0);
-          }
-    }
-
-  is_reinitialized = true;
-}
-
-
-
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluationBase<n_components_,
-                                            dim,
-                                            spacedim,
-                                            Number>::value_type &
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_value(
-  const unsigned int point_index) const
-{
-  AssertIndexRange(point_index, values.size());
-  return values[point_index];
-}
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used during
+   * `FEEvaluation::set_dof_values(global_vector)` or
+   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  test_and_sum(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values = false);
 
+  /**
+   * This function multiplies the quantities passed in by previous
+   * submit_value() or submit_gradient() calls by the value or gradient of the
+   * test functions, and performs summation over all given points multiplied be
+   * the Jacobian determinant times the quadrature weight (JxW).
+   *
+   * @param[out] solution_values This array will contain the result of the
+   * integral, which can be used to during
+   * `cell->set_dof_values(solution_values, global_vector)` or
+   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
+   * that for multi-component systems where only some of the components are
+   * selected by the present class, the entries in `solution_values` not touched
+   * by this class will be set to zero.
+   *
+   * @param[in] integration_flags Flags specifying which quantities should be
+   * integrated at the points.
+   *
+   * @param[in] sum_into_values Flag specifying if the integrated values
+   * should be summed into the solution values. For the default value
+   * `sum_into_values=false` every value of @p solution_values is zeroed out.
+   *
+   */
+  void
+  test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
+               const EvaluationFlags::EvaluationFlags &integration_flags,
+               const bool                              sum_into_values = false);
 
+  /**
+   * Evaluate values and gradients in face for the selected face (lane) of the
+   * batch. Default stride into the face dofs is width of
+   * VectorizedArray<selected_floating_point_type> which is the default
+   * vectorization over faces for FEFaceEvaluation.
+   */
+  template <int stride_face_dof = VectorizedArrayType::size()>
+  void
+  evaluate_in_face(const ScalarNumber                     *face_dof_values,
+                   const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluationBase<n_components_,
-                                            dim,
-                                            spacedim,
-                                            Number>::gradient_type &
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_gradient(
-  const unsigned int point_index) const
-{
-  AssertIndexRange(point_index, gradients.size());
-  return gradients[point_index];
-}
+  /**
+   * Integrate values and gradients in face for the selected face (lane) of the
+   * batch. Default stride into the face dofs is width of
+   * VectorizedArray<selected_floating_point_type> which is the default
+   * vectorization over faces for FEFaceEvaluation.
+   */
+  template <int stride_face_dof = VectorizedArrayType::size()>
+  void
+  integrate_in_face(ScalarNumber                           *face_dof_values,
+                    const EvaluationFlags::EvaluationFlags &integration_flags,
+                    const bool sum_into_values = false);
 
+  /**
+   * Return the normal vector. This class or the MappingInfo object passed to
+   * this function needs to be constructed with UpdateFlags containing
+   * `update_normal_vectors`.
+   */
+  Tensor<1, spacedim, Number>
+  normal_vector(const unsigned int point_index) const;
 
+private:
+  static constexpr std::size_t n_lanes_user_interface =
+    internal::VectorizedArrayTrait<Number>::width();
+  static constexpr std::size_t n_lanes_internal =
+    internal::VectorizedArrayTrait<VectorizedArrayType>::width();
+  static constexpr std::size_t stride =
+    internal::VectorizedArrayTrait<Number>::stride();
 
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_value(
-  const value_type  &value,
-  const unsigned int point_index)
-{
-  AssertIndexRange(point_index, n_q_points);
-  values[point_index] = value;
-}
+  template <bool is_linear, std::size_t stride_view>
+  void
+  do_evaluate(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
+  template <bool do_JxW, bool is_linear, std::size_t stride_view>
+  void
+  do_integrate(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    const bool                                         sum_into_values);
 
+  /**
+   * Actually does the evaluation templated on the chosen code path (linear or
+   * higher order).
+   */
+  template <bool is_linear, int stride_face_dof>
+  void
+  do_evaluate_in_face(const ScalarNumber                     *face_dof_values,
+                      const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
-template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_gradient(
-  const gradient_type &gradient,
-  const unsigned int   point_index)
-{
-  AssertIndexRange(point_index, n_q_points);
-  gradients[point_index] = gradient;
-}
+  /**
+   * Actually does the integration templated on the chosen code path (linear or
+   * higher order).
+   */
+  template <bool do_JxW, bool is_linear, int stride_face_dof>
+  void
+  do_integrate_in_face(
+    ScalarNumber                           *face_dof_values,
+    const EvaluationFlags::EvaluationFlags &integration_flags,
+    const bool                              sum_into_values);
+};
 
+// ----------------------- template and inline function ----------------------
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, dim, spacedim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::jacobian(
-  const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  FEPointEvaluationBase(const Mapping<dim>       &mapping,
+                        const FiniteElement<dim> &fe,
+                        const UpdateFlags         update_flags,
+                        const unsigned int        first_selected_component)
+  : n_q_batches(numbers::invalid_unsigned_int)
+  , n_q_points(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(numbers::invalid_unsigned_int)
+  , mapping(&mapping)
+  , fe(&fe)
+  , JxW_ptr(nullptr)
+  , update_flags(update_flags)
+  , mapping_info_on_the_fly(
+      std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
+        mapping,
+        update_flags))
+  , mapping_info(mapping_info_on_the_fly.get())
+  , current_cell_index(numbers::invalid_unsigned_int)
+  , current_face_number(numbers::invalid_unsigned_int)
+  , is_reinitialized(false)
+  , is_interior(true)
 {
-  AssertIndexRange(point_index, n_q_points);
-  Assert(jacobian_ptr != nullptr,
-         internal::FEPointEvaluation::
-           ExcFEPointEvaluationAccessToUninitializedMappingField(
-             "update_jacobians"));
-  return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions::
-                                     GeometryType::affine ?
-                        0 :
-                        point_index];
+  setup(first_selected_component);
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline DerivativeForm<1, spacedim, dim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::inverse_jacobian(
-  const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  FEPointEvaluationBase(
+    NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+    const FiniteElement<dim>                        &fe,
+    const unsigned int                               first_selected_component,
+    const bool                                       is_interior)
+  : n_q_batches(numbers::invalid_unsigned_int)
+  , n_q_points(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(numbers::invalid_unsigned_int)
+  , mapping(&mapping_info.get_mapping())
+  , fe(&fe)
+  , JxW_ptr(nullptr)
+  , update_flags(mapping_info.get_update_flags())
+  , mapping_info(&mapping_info)
+  , current_cell_index(numbers::invalid_unsigned_int)
+  , current_face_number(numbers::invalid_unsigned_int)
+  , is_reinitialized(false)
+  , is_interior(is_interior)
 {
-  AssertIndexRange(point_index, n_q_points);
-  Assert(inverse_jacobian_ptr != nullptr,
-         internal::FEPointEvaluation::
-           ExcFEPointEvaluationAccessToUninitializedMappingField(
-             "update_inverse_jacobians"));
-  return inverse_jacobian_ptr
-    [cell_type <=
-         ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ?
-       0 :
-       point_index];
+  setup(first_selected_component);
+  connection_is_reinitialized = mapping_info.connect_is_reinitialized(
+    [this]() { this->is_reinitialized = false; });
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline Number
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::JxW(
-  const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  FEPointEvaluationBase(
+    FEPointEvaluationBase<n_components_, dim, spacedim, Number> &other) noexcept
+  : n_q_batches(other.n_q_batches)
+  , n_q_points(other.n_q_points)
+  , n_q_points_scalar(other.n_q_points_scalar)
+  , mapping(other.mapping)
+  , fe(other.fe)
+  , poly(other.poly)
+  , use_linear_path(other.use_linear_path)
+  , renumber(other.renumber)
+  , solution_renumbered(other.solution_renumbered)
+  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
+  , values(other.values)
+  , gradients(other.gradients)
+  , dofs_per_component(other.dofs_per_component)
+  , dofs_per_component_face(other.dofs_per_component_face)
+  , component_in_base_element(other.component_in_base_element)
+  , nonzero_shape_function_component(other.nonzero_shape_function_component)
+  , update_flags(other.update_flags)
+  , fe_values(other.fe_values)
+  , mapping_info_on_the_fly(
+      other.mapping_info_on_the_fly ?
+        std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
+          *mapping,
+          update_flags) :
+        nullptr)
+  , mapping_info(other.mapping_info)
+  , current_cell_index(other.current_cell_index)
+  , current_face_number(other.current_face_number)
+  , fast_path(other.fast_path)
+  , is_reinitialized(false)
+  , shapes(other.shapes)
+  , shapes_faces(other.shapes_faces)
+  , is_interior(other.is_interior)
 {
-  AssertIndexRange(point_index, n_q_points);
-  Assert(JxW_ptr != nullptr,
-         internal::FEPointEvaluation::
-           ExcFEPointEvaluationAccessToUninitializedMappingField(
-             "update_JxW_values"));
-  return JxW_ptr[point_index];
+  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
+    [this]() { this->is_reinitialized = false; });
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline Point<spacedim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::real_point(
-  const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  FEPointEvaluationBase(
+    FEPointEvaluationBase<n_components_, dim, spacedim, Number>
+      &&other) noexcept
+  : n_q_batches(other.n_q_batches)
+  , n_q_points(other.n_q_points)
+  , n_q_points_scalar(other.n_q_points_scalar)
+  , mapping(other.mapping)
+  , fe(other.fe)
+  , poly(other.poly)
+  , use_linear_path(other.use_linear_path)
+  , renumber(other.renumber)
+  , solution_renumbered(other.solution_renumbered)
+  , solution_renumbered_vectorized(other.solution_renumbered_vectorized)
+  , values(other.values)
+  , gradients(other.gradients)
+  , dofs_per_component(other.dofs_per_component)
+  , dofs_per_component_face(other.dofs_per_component_face)
+  , component_in_base_element(other.component_in_base_element)
+  , nonzero_shape_function_component(other.nonzero_shape_function_component)
+  , update_flags(other.update_flags)
+  , fe_values(other.fe_values)
+  , mapping_info_on_the_fly(std::move(other.mapping_info_on_the_fly))
+  , mapping_info(other.mapping_info)
+  , current_cell_index(other.current_cell_index)
+  , current_face_number(other.current_face_number)
+  , fast_path(other.fast_path)
+  , is_reinitialized(false)
+  , shapes(other.shapes)
+  , shapes_faces(other.shapes_faces)
+  , is_interior(other.is_interior)
 {
-  AssertIndexRange(point_index, n_q_points);
-  Assert(real_point_ptr != nullptr,
-         internal::FEPointEvaluation::
-           ExcFEPointEvaluationAccessToUninitializedMappingField(
-             "update_quadrature_points"));
-  return real_point_ptr[point_index];
+  connection_is_reinitialized = mapping_info->connect_is_reinitialized(
+    [this]() { this->is_reinitialized = false; });
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline Point<dim, Number>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::unit_point(
-  const unsigned int point_index) const
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  ~FEPointEvaluationBase()
 {
-  AssertIndexRange(point_index, n_q_points);
-  Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
-  Point<dim, Number> unit_point;
-  for (unsigned int d = 0; d < dim; ++d)
-    unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
-      unit_point_ptr[point_index / stride][d], point_index % stride);
-  return unit_point;
+  connection_is_reinitialized.disconnect();
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
-FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
-  quadrature_point_indices() const
+void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::setup(
+  const unsigned int first_selected_component)
 {
-  return {0U, n_q_points};
-}
+  AssertIndexRange(first_selected_component + n_components,
+                   fe->n_components() + 1);
 
+  shapes.reserve(100);
 
+  bool         same_base_element   = true;
+  unsigned int base_element_number = 0;
+  component_in_base_element        = 0;
+  unsigned int component           = 0;
+  for (; base_element_number < fe->n_base_elements(); ++base_element_number)
+    if (component + fe->element_multiplicity(base_element_number) >
+        first_selected_component)
+      {
+        if (first_selected_component + n_components >
+            component + fe->element_multiplicity(base_element_number))
+          same_base_element = false;
+        component_in_base_element = first_selected_component - component;
+        break;
+      }
+    else
+      component += fe->element_multiplicity(base_element_number);
 
-/**
- * This class provides an interface to the evaluation of interpolated solution
- * values and gradients on cells on arbitrary reference point positions. These
- * points can change from cell to cell, both with respect to their quantity as
- * well to the location. The two typical use cases are evaluations on
- * non-matching grids and particle simulations.
- *
- * The use of this class is similar to FEValues or FEEvaluation: The class is
- * first initialized to a cell by calling `FEPointEvaluation::reinit(cell,
- * unit_points)`, with the main difference to the other concepts that the
- * underlying points in reference coordinates need to be passed along. Then,
- * upon call to evaluate() or integrate(), the user can compute information at
- * the give points. Eventually, the access functions get_value() or
- * get_gradient() allow to query this information at a specific point index.
- *
- * The functionality is similar to creating an FEValues object with a
- * Quadrature object on the `unit_points` on every cell separately and then
- * calling FEValues::get_function_values or FEValues::get_function_gradients,
- * and for some elements and mappings this is what actually happens
- * internally. For specific combinations of Mapping and FiniteElement
- * realizations, however, there is a much more efficient implementation that
- * avoids the memory allocation and other expensive start-up cost of
- * FEValues. Currently, the functionality is specialized for mappings derived
- * from MappingQ and MappingCartesian and for finite elements with tensor
- * product structure that work with the
- * @ref matrixfree
- * module. In those cases, the cost implied
- * by this class is similar (or sometimes even somewhat lower) than using
- * `FEValues::reinit(cell)` followed by `FEValues::get_function_gradients`.
- */
-template <int n_components_,
-          int dim,
-          int spacedim    = dim,
-          typename Number = double>
-class FEPointEvaluation
-  : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-{
-public:
-  static constexpr unsigned int dimension    = dim;
-  static constexpr unsigned int n_components = n_components_;
+  if (internal::FEPointEvaluation::is_fast_path_supported(*mapping) &&
+      internal::FEPointEvaluation::is_fast_path_supported(
+        *fe, base_element_number) &&
+      same_base_element)
+    {
+      shape_info.reinit(QMidpoint<1>(), *fe, base_element_number);
+      renumber                = shape_info.lexicographic_numbering;
+      dofs_per_component      = shape_info.dofs_per_component_on_cell;
+      dofs_per_component_face = shape_info.dofs_per_component_on_face;
+      poly = internal::FEPointEvaluation::get_polynomial_space(
+        fe->base_element(base_element_number));
 
-  using number_type = Number;
+      bool is_lexicographic = true;
+      for (unsigned int i = 0; i < renumber.size(); ++i)
+        if (i != renumber[i])
+          is_lexicographic = false;
 
-  using ScalarNumber =
-    typename internal::VectorizedArrayTrait<Number>::value_type;
-  using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
-    Number>::vectorized_value_type;
-  using ETT = typename internal::FEPointEvaluation::
-    EvaluatorTypeTraits<dim, n_components, Number>;
-  using value_type            = typename ETT::value_type;
-  using scalar_value_type     = typename ETT::scalar_value_type;
-  using vectorized_value_type = typename ETT::vectorized_value_type;
-  using gradient_type         = typename ETT::gradient_type;
-  using interface_vectorized_gradient_type =
-    typename ETT::interface_vectorized_gradient_type;
+      if (is_lexicographic)
+        renumber.clear();
 
-  /**
-   * Constructor.
-   *
-   * @param mapping The Mapping class describing the actual geometry of a cell
-   * passed to the evaluate() function.
-   *
-   * @param fe The FiniteElement object that is used for the evaluation, which
-   * is typically the same on all cells to be evaluated.
-   *
-   * @param update_flags Specify the quantities to be computed by the mapping
-   * during the call of reinit(). During evaluate() or integrate(), this data
-   * is queried to produce the desired result (e.g., the gradient of a finite
-   * element solution).
-   *
-   * @param first_selected_component For multi-component FiniteElement
-   * objects, this parameter allows to select a range of `n_components`
-   * components starting from this parameter.
-   */
-  FEPointEvaluation(const Mapping<dim>       &mapping,
-                    const FiniteElement<dim> &fe,
-                    const UpdateFlags         update_flags,
-                    const unsigned int        first_selected_component = 0)
-    : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
-        mapping,
-        fe,
-        update_flags,
-        first_selected_component)
-  {}
+      use_linear_path = (poly.size() == 2 && poly[0].value(0.) == 1. &&
+                         poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
+                         poly[1].value(1.) == 1.) &&
+                        (fe->n_components() == n_components);
 
-  /**
-   * Constructor to make the present class able to re-use the geometry
-   * data also used by other `FEPointEvaluation` objects.
-   *
-   * @param mapping_info The MappingInfo class describes the geometry-related
-   * data for evaluating finite-element solutions. This object enables to
-   * construct such an object on the outside, possibly re-using it between
-   * several objects or between several calls to the same cell and unit points.
-   *
-   * @param fe The FiniteElement object that is used for the evaluation, which
-   * is typically the same on all cells to be evaluated.
-   *
-   * @param first_selected_component For multi-component FiniteElement
-   * objects, this parameter allows to select a range of `n_components`
-   * components starting from this parameter.
-   */
-  FEPointEvaluation(
-    NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
-    const unsigned int first_selected_component = 0)
-    : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
-        mapping_info,
-        fe,
-        first_selected_component)
-  {}
+      const unsigned int size_face = 3 * dofs_per_component_face * n_components;
+      const unsigned int size_cell = dofs_per_component * n_components;
+      scratch_data_scalar.resize(size_face + size_cell);
 
-  /**
-   * Set up the mapping information for the given cell, e.g., by computing the
-   * Jacobian of the mapping for the given points if gradients of the functions
-   * are requested.
-   *
-   * @param[in] cell An iterator to the current cell
-   *
-   * @param[in] unit_points List of points in the reference locations of the
-   * current cell where the FiniteElement object should be
-   * evaluated/integrated in the evaluate() and integrate() functions.
-   */
-  void
-  reinit(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-         const ArrayView<const Point<dim>> &unit_points);
-
-  /**
-   * Reinitialize the evaluator to point to the correct precomputed mapping of
-   * the single cell in the MappingInfo object.
-   */
-  void
-  reinit();
+      solution_renumbered.resize(dofs_per_component);
+      solution_renumbered_vectorized.resize(dofs_per_component);
 
-  /**
-   * Reinitialize the evaluator to point to the correct precomputed mapping of
-   * the cell in the MappingInfo object.
-   */
-  void
-  reinit(const unsigned int cell_index);
+      fast_path = true;
+    }
+  else
+    {
+      nonzero_shape_function_component.resize(fe->n_dofs_per_cell());
+      for (unsigned int d = 0; d < n_components; ++d)
+        {
+          const unsigned int component = first_selected_component + d;
+          for (unsigned int i = 0; i < fe->n_dofs_per_cell(); ++i)
+            {
+              const bool is_primitive =
+                fe->is_primitive() || fe->is_primitive(i);
+              if (is_primitive)
+                nonzero_shape_function_component[i][d] =
+                  (component == fe->system_to_component_index(i).first);
+              else
+                nonzero_shape_function_component[i][d] =
+                  (fe->get_nonzero_components(i)[component] == true);
+            }
+        }
 
+      fast_path = false;
+    }
+}
 
-  /**
-   * This function interpolates the finite element solution, represented by
-   * `solution_values`, on the cell and `unit_points` passed to reinit().
-   *
-   * @param[in] solution_values This array is supposed to contain the unknown
-   * values on the element read out by
-   * `FEEvaluation::read_dof_values(global_vector)`.
-   *
-   * @param[in] evaluation_flags Flags specifying which quantities should be
-   * evaluated at the points.
-   */
-  template <std::size_t stride_view>
-  void
-  evaluate(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
-  /**
-   * This function interpolates the finite element solution, represented by
-   * `solution_values`, on the cell and `unit_points` passed to reinit().
-   *
-   * @param[in] solution_values This array is supposed to contain the unknown
-   * values on the element as returned by `cell->get_dof_values(global_vector,
-   * solution_values)`.
-   *
-   * @param[in] evaluation_flags Flags specifying which quantities should be
-   * evaluated at the points.
-   */
-  void
-  evaluate(const ArrayView<const ScalarNumber>    &solution_values,
-           const EvaluationFlags::EvaluationFlags &evaluation_flags);
 
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used during
-   * `FEEvaluation::set_dof_values(global_vector)` or
-   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  template <std::size_t stride_view>
-  void
-  integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags,
-            const bool                              sum_into_values = false);
+template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face, bool is_linear>
+inline void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::do_reinit()
+{
+  const unsigned int geometry_index =
+    mapping_info->template compute_geometry_index_offset<is_face>(
+      current_cell_index, current_face_number);
 
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used to during
-   * `cell->set_dof_values(solution_values, global_vector)` or
-   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  void
-  integrate(const ArrayView<ScalarNumber>          &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags,
-            const bool                              sum_into_values = false);
+  cell_type = mapping_info->get_cell_type(geometry_index);
 
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points. This is
-   * similar to the integration of a bilinear form in terms of the test
-   * function, with the difference that this formula does not include a `JxW`
-   * factor (in contrast to the integrate function of this class). This allows
-   * the class to naturally embed point information (e.g. particles) into a
-   * finite element formulation.
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used during
-   * `FEEvaluation::set_dof_values(global_vector)` or
-   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  template <std::size_t stride_view>
-  void
-  test_and_sum(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values = false);
+  const_cast<unsigned int &>(n_q_points_scalar) =
+    mapping_info->get_n_q_points_unvectorized(geometry_index);
 
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points. This is
-   * similar to the integration of a bilinear form in terms of the test
-   * function, with the difference that this formula does not include a `JxW`
-   * factor (in contrast to the integrate function of this class). This allows
-   * the class to naturally embed point information (e.g. particles) into a
-   * finite element formulation.
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used during
-   * `cell->set_dof_values(solution_values, global_vector)` or
-   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  void
-  test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
-               const EvaluationFlags::EvaluationFlags &integration_flags,
-               const bool                              sum_into_values = false);
+  // round up n_q_points_scalar / n_lanes_internal
+  const_cast<unsigned int &>(n_q_batches) =
+    (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
 
-  /**
-   * Return the normal vector. This class or the MappingInfo object passed to
-   * this function needs to be constructed with UpdateFlags containing
-   * `update_normal_vectors`.
-   */
-  Tensor<1, spacedim, Number>
-  normal_vector(const unsigned int point_index) const;
+  const unsigned int n_q_points_before = n_q_points;
 
-private:
-  static constexpr std::size_t n_lanes_user_interface =
-    internal::VectorizedArrayTrait<Number>::width();
-  static constexpr std::size_t n_lanes_internal =
-    internal::VectorizedArrayTrait<VectorizedArrayType>::width();
-  static constexpr std::size_t stride =
-    internal::VectorizedArrayTrait<Number>::stride();
+  const_cast<unsigned int &>(n_q_points) =
+    (stride == 1) ? n_q_batches : n_q_points_scalar;
 
-  /**
-   * Resizes necessary data fields, reads in and renumbers solution values.
-   * Interpolates onto face if face path is selected.
-   */
-  template <bool is_linear, std::size_t stride_view>
-  void
-  prepare_evaluate_fast(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values);
+  if (n_q_points != n_q_points_before)
+    {
+      if (update_flags & update_values)
+        values.resize(n_q_points);
+      if (update_flags & update_gradients)
+        gradients.resize(n_q_points);
+    }
 
-  /**
-   * Evaluates the actual interpolation on the cell or face for a quadrature
-   * batch.
-   */
-  template <bool is_linear, std::size_t stride_view>
-  void
-  compute_evaluate_fast(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags,
-    const unsigned int                                       n_shapes,
-    const unsigned int                                       qb,
-    vectorized_value_type                                   &value,
-    interface_vectorized_gradient_type                      &gradient);
+  if (n_q_points == 0)
+    {
+      is_reinitialized = true;
+      return;
+    }
 
-  /**
-   * Fast path of the evaluate function.
-   */
-  template <bool is_linear, std::size_t stride_view>
-  void
-  evaluate_fast(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
+  // set unit point pointer
+  const unsigned int unit_point_offset =
+    mapping_info->compute_unit_point_index_offset(geometry_index);
 
-  /**
-   * Slow path of the evaluate function using FEValues.
-   */
-  template <std::size_t stride_view>
-  void
-  evaluate_slow(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
-
-  /**
-   * Integrates the product of the data passed in by submit_value() and
-   * submit_gradient() with the values or gradients of test functions on the
-   * cell or face for a given quadrature batch.
-   */
-  template <bool is_linear>
-  void
-  compute_integrate_fast(
-    const EvaluationFlags::EvaluationFlags  &integration_flags,
-    const unsigned int                       n_shapes,
-    const unsigned int                       qb,
-    const vectorized_value_type              value,
-    const interface_vectorized_gradient_type gradient,
-    vectorized_value_type                   *solution_values_vectorized_linear);
+  if (is_face)
+    unit_point_faces_ptr =
+      mapping_info->get_unit_point_faces(unit_point_offset);
+  else
+    unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
 
-  /**
-   * Addition across the lanes of VectorizedArray as accumulated by the
-   * compute_integrate_fast_function(), writing the sum into the result vector.
-   * Applies face contributions to cell contributions for face path.
-   */
-  template <bool is_linear, std::size_t stride_view>
-  void
-  finish_integrate_fast(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    vectorized_value_type *solution_values_vectorized_linear,
-    const bool             sum_into_values);
+  // set data pointers
+  const unsigned int data_offset =
+    mapping_info->compute_data_index_offset(geometry_index);
+  const unsigned int compressed_data_offset =
+    mapping_info->compute_compressed_data_index_offset(geometry_index);
+#ifdef DEBUG
+  const UpdateFlags update_flags_mapping =
+    mapping_info->get_update_flags_mapping();
+  if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+    real_point_ptr = mapping_info->get_real_point(data_offset);
+  if (update_flags_mapping & UpdateFlags::update_jacobians)
+    jacobian_ptr =
+      mapping_info->get_jacobian(compressed_data_offset, is_interior);
+  if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+    inverse_jacobian_ptr =
+      mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
+  if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+    normal_ptr = mapping_info->get_normal_vector(data_offset);
+  if (update_flags_mapping & UpdateFlags::update_JxW_values)
+    JxW_ptr = mapping_info->get_JxW(data_offset);
+#else
+  real_point_ptr = mapping_info->get_real_point(data_offset);
+  jacobian_ptr =
+    mapping_info->get_jacobian(compressed_data_offset, is_interior);
+  inverse_jacobian_ptr =
+    mapping_info->get_inverse_jacobian(compressed_data_offset, is_interior);
+  normal_ptr = mapping_info->get_normal_vector(data_offset);
+  JxW_ptr    = mapping_info->get_JxW(data_offset);
+#endif
 
-  /**
-   * Fast path of the integrate function.
-   */
-  template <bool do_JxW, bool is_linear, std::size_t stride_view>
-  void
-  integrate_fast(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values);
+  if (!is_linear && fast_path)
+    {
+      const std::size_t n_shapes = poly.size();
 
-  /**
-   * Slow path of the integrate function using FEValues.
-   */
-  template <bool do_JxW, std::size_t stride_view>
-  void
-  integrate_slow(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values);
+      for (unsigned int qb = 0; qb < n_q_batches; ++qb)
+        if (is_face)
+          {
+            if (dim > 1)
+              {
+                shapes_faces.resize_fast(n_q_batches * n_shapes);
+                internal::compute_values_of_array(
+                  shapes_faces.data() + qb * n_shapes,
+                  poly,
+                  unit_point_faces_ptr[qb],
+                  update_flags & UpdateFlags::update_gradients ? 1 : 0);
+              }
+          }
+        else
+          {
+            shapes.resize_fast(n_q_batches * n_shapes);
+            internal::compute_values_of_array(
+              shapes.data() + qb * n_shapes,
+              poly,
+              unit_point_ptr[qb],
+              update_flags & UpdateFlags::update_gradients ? 1 : 0);
+          }
+    }
 
-  /**
-   * Implementation of the integrate/test_and_sum function.
-   */
-  template <bool do_JxW, std::size_t stride_view>
-  void
-  do_integrate(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values);
-};
+  is_reinitialized = true;
+}
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
+inline const typename FEPointEvaluationBase<n_components_,
+                                            dim,
+                                            spacedim,
+                                            Number>::value_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_value(
+  const unsigned int point_index) const
 {
-  this->current_cell_index  = numbers::invalid_unsigned_int;
-  this->current_face_number = numbers::invalid_unsigned_int;
-
-  if (this->use_linear_path)
-    this->template do_reinit<false, true>();
-  else
-    this->template do_reinit<false, false>();
+  AssertIndexRange(point_index, values.size());
+  return values[point_index];
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
-  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-  const ArrayView<const Point<dim>>                          &unit_points)
+inline const typename FEPointEvaluationBase<n_components_,
+                                            dim,
+                                            spacedim,
+                                            Number>::gradient_type &
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::get_gradient(
+  const unsigned int point_index) const
 {
-  // reinit is only allowed for mapping computation on the fly
-  AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
-              ExcNotImplemented());
+  AssertIndexRange(point_index, gradients.size());
+  return gradients[point_index];
+}
 
-  this->mapping_info->reinit(cell, unit_points);
 
-  if (!this->fast_path)
-    {
-      this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
-        *this->mapping,
-        *this->fe,
-        Quadrature<dim>(
-          std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
-        this->update_flags);
-      this->fe_values->reinit(cell);
-    }
 
-  if (this->use_linear_path)
-    this->template do_reinit<false, true>();
-  else
-    this->template do_reinit<false, false>();
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_value(
+  const value_type  &value,
+  const unsigned int point_index)
+{
+  AssertIndexRange(point_index, n_q_points);
+  values[point_index] = value;
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
 inline void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
-  const unsigned int cell_index)
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::submit_gradient(
+  const gradient_type &gradient,
+  const unsigned int   point_index)
 {
-  this->current_cell_index  = cell_index;
-  this->current_face_number = numbers::invalid_unsigned_int;
-
-  if (this->use_linear_path)
-    this->template do_reinit<false, true>();
-  else
-    this->template do_reinit<false, false>();
-
-  if (!this->fast_path)
-    {
-      std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
+  AssertIndexRange(point_index, n_q_points);
+  gradients[point_index] = gradient;
+}
 
-      for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
-        for (unsigned int d = 0; d < dim; ++d)
-          unit_points[v][d] =
-            this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
 
-      this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
-        *this->mapping,
-        *this->fe,
-        Quadrature<dim>(
-          std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
-        this->update_flags);
 
-      this->fe_values->reinit(
-        this->mapping_info->get_cell_iterator(this->current_cell_index));
-    }
+template <int n_components_, int dim, int spacedim, typename Number>
+inline DerivativeForm<1, dim, spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::jacobian(
+  const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, n_q_points);
+  Assert(jacobian_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_jacobians"));
+  return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions::
+                                     GeometryType::affine ?
+                        0 :
+                        point_index];
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-template <std::size_t stride_view>
-void
-FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
-  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
+inline DerivativeForm<1, spacedim, dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::inverse_jacobian(
+  const unsigned int point_index) const
 {
-  if (!this->is_reinitialized)
-    reinit();
-
-  if (this->n_q_points == 0)
-    return;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(inverse_jacobian_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_inverse_jacobians"));
+  return inverse_jacobian_ptr
+    [cell_type <=
+         ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ?
+       0 :
+       point_index];
+}
 
-  Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
 
-  if (!((evaluation_flags & EvaluationFlags::values) ||
-        (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
-    return;
 
-  AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
-  if (this->fast_path)
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Number
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::JxW(
+  const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, n_q_points);
+  Assert(JxW_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_JxW_values"));
+  return JxW_ptr[point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Point<spacedim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::real_point(
+  const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, n_q_points);
+  Assert(real_point_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_quadrature_points"));
+  return real_point_ptr[point_index];
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline Point<dim, Number>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::unit_point(
+  const unsigned int point_index) const
+{
+  AssertIndexRange(point_index, n_q_points);
+  Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
+  Point<dim, Number> unit_point;
+  for (unsigned int d = 0; d < dim; ++d)
+    unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+      unit_point_ptr[point_index / stride][d], point_index % stride);
+  return unit_point;
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+FEPointEvaluationBase<n_components_, dim, spacedim, Number>::
+  quadrature_point_indices() const
+{
+  return {0U, n_q_points};
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
+  NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+  const FiniteElement<dim>                        &fe,
+  const unsigned int                               first_selected_component)
+  : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+      mapping_info,
+      fe,
+      first_selected_component)
+{}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
+  const Mapping<dim>       &mapping,
+  const FiniteElement<dim> &fe,
+  const UpdateFlags         update_flags,
+  const unsigned int        first_selected_component)
+  : FEPointEvaluationBase<n_components_, dim, spacedim, Number>(
+      mapping,
+      fe,
+      update_flags,
+      first_selected_component)
+{}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
+{
+  this->current_cell_index  = numbers::invalid_unsigned_int;
+  this->current_face_number = numbers::invalid_unsigned_int;
+
+  if (this->use_linear_path)
+    this->template do_reinit<false, true>();
+  else
+    this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+  const ArrayView<const Point<dim>>                          &unit_points)
+{
+  // reinit is only allowed for mapping computation on the fly
+  AssertThrow(this->mapping_info_on_the_fly.get() != nullptr,
+              ExcNotImplemented());
+
+  this->mapping_info->reinit(cell, unit_points);
+
+  if (!this->fast_path)
+    {
+      this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
+        *this->mapping,
+        *this->fe,
+        Quadrature<dim>(
+          std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+        this->update_flags);
+      this->fe_values->reinit(cell);
+    }
+
+  if (this->use_linear_path)
+    this->template do_reinit<false, true>();
+  else
+    this->template do_reinit<false, false>();
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+inline void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
+  const unsigned int cell_index)
+{
+  this->current_cell_index  = cell_index;
+  this->current_face_number = numbers::invalid_unsigned_int;
+
+  if (this->use_linear_path)
+    this->template do_reinit<false, true>();
+  else
+    this->template do_reinit<false, false>();
+
+  if (!this->fast_path)
+    {
+      std::vector<Point<dim>> unit_points(this->n_q_points_scalar);
+
+      for (unsigned int v = 0; v < this->n_q_points_scalar; ++v)
+        for (unsigned int d = 0; d < dim; ++d)
+          unit_points[v][d] =
+            this->unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
+
+      this->fe_values = std::make_shared<FEValues<dim, spacedim>>(
+        *this->mapping,
+        *this->fe,
+        Quadrature<dim>(
+          std::vector<Point<dim>>(unit_points.begin(), unit_points.end())),
+        this->update_flags);
+
+      this->fe_values->reinit(
+        this->mapping_info->get_cell_iterator(this->current_cell_index));
+    }
+}
+
+
+
+template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
+void
+FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
+  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
+{
+  if (!this->is_reinitialized)
+    reinit();
+
+  if (this->n_q_points == 0)
+    return;
+
+  Assert(!(evaluation_flags & EvaluationFlags::hessians), ExcNotImplemented());
+
+  if (!((evaluation_flags & EvaluationFlags::values) ||
+        (evaluation_flags & EvaluationFlags::gradients))) // no evaluation flags
+    return;
+
+  AssertDimension(solution_values.size(), this->fe->dofs_per_cell);
+  if (this->fast_path)
     {
       if (this->use_linear_path)
         evaluate_fast<true>(solution_values, evaluation_flags);
@@ -2766,286 +3071,6 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::normal_vector(
 
 
 
-/**
- * This class provides an interface to the evaluation of interpolated solution
- * values and gradients on faces on arbitrary reference point positions. These
- * points can change from face to face, both with respect to their quantity as
- * well to the location. A typical use case is evaluations on non-matching
- * grids.
- *
- * The use of this class is similar to FEEvaluation: In the constructor, a
- * reference to a NonMatching::MappingInfo object is passed, where the
- * quadrature points in reference position is stored together with the mapping
- * information. The class is then reinitialized to a cell by calling
- * `FEFacePointEvaluation::reinit(face_index)` or
- * `FEFacePointEvaluation::reinit(cell_index, face_number)`. Then, upon call to
- * evaluate() or integrate(), the user can compute information at the given
- * points. Eventually, the access functions get_value() or get_gradient() allow
- * to query this information at a specific point index.
- */
-template <int n_components_,
-          int dim,
-          int spacedim    = dim,
-          typename Number = double>
-class FEFacePointEvaluation
-  : public FEPointEvaluationBase<n_components_, dim, spacedim, Number>
-{
-public:
-  static constexpr unsigned int dimension    = dim;
-  static constexpr unsigned int n_components = n_components_;
-
-  using number_type = Number;
-
-  using ScalarNumber =
-    typename internal::VectorizedArrayTrait<Number>::value_type;
-  using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
-    Number>::vectorized_value_type;
-  using ETT = typename internal::FEPointEvaluation::
-    EvaluatorTypeTraits<dim, n_components, Number>;
-  using value_type            = typename ETT::value_type;
-  using scalar_value_type     = typename ETT::scalar_value_type;
-  using vectorized_value_type = typename ETT::vectorized_value_type;
-  using gradient_type         = typename ETT::gradient_type;
-  using interface_vectorized_gradient_type =
-    typename ETT::interface_vectorized_gradient_type;
-
-  /**
-   * Constructor. Allows to select if interior or exterior face is selected.
-   */
-  FEFacePointEvaluation(
-    NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
-    const FiniteElement<dim>                        &fe,
-    const bool                                       is_interior = true,
-    const unsigned int first_selected_component                  = 0);
-
-  /**
-   * Reinitialize the evaluator to point to the correct precomputed mapping of
-   * the face in the MappingInfo object. Used in element-centric loops (ECL).
-   */
-  void
-  reinit(const unsigned int cell_index, const unsigned int face_number);
-
-  /**
-   * Reinitialize the evaluator to point to the correct precomputed mapping of
-   * the face in the MappingInfo object. Used in face-centric loops (FCL).
-   */
-  void
-  reinit(const unsigned int face_index);
-
-  /**
-   * This function interpolates the finite element solution, represented by
-   * `solution_values`, on the cell and `unit_points` passed to reinit().
-   *
-   * @param[in] solution_values This array is supposed to contain the unknown
-   * values on the element read out by
-   * `FEEvaluation::read_dof_values(global_vector)`.
-   *
-   * @param[in] evaluation_flags Flags specifying which quantities should be
-   * evaluated at the points.
-   */
-  template <std::size_t stride_view>
-  void
-  evaluate(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
-
-  /**
-   * This function interpolates the finite element solution, represented by
-   * `solution_values`, on the cell and `unit_points` passed to reinit().
-   *
-   * @param[in] solution_values This array is supposed to contain the unknown
-   * values on the element as returned by `cell->get_dof_values(global_vector,
-   * solution_values)`.
-   *
-   * @param[in] evaluation_flags Flags specifying which quantities should be
-   * evaluated at the points.
-   */
-  void
-  evaluate(const ArrayView<const ScalarNumber>    &solution_values,
-           const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used during
-   * `FEEvaluation::set_dof_values(global_vector)` or
-   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  template <std::size_t stride_view>
-  void
-  integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags,
-            const bool                              sum_into_values = false);
-
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used to during
-   * `cell->set_dof_values(solution_values, global_vector)` or
-   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  void
-  integrate(const ArrayView<ScalarNumber>          &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags,
-            const bool                              sum_into_values = false);
-
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used during
-   * `FEEvaluation::set_dof_values(global_vector)` or
-   * `FEEvaluation::distribute_local_to_global(global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  template <std::size_t stride_view>
-  void
-  test_and_sum(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values = false);
-
-  /**
-   * This function multiplies the quantities passed in by previous
-   * submit_value() or submit_gradient() calls by the value or gradient of the
-   * test functions, and performs summation over all given points multiplied be
-   * the Jacobian determinant times the quadrature weight (JxW).
-   *
-   * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used to during
-   * `cell->set_dof_values(solution_values, global_vector)` or
-   * `cell->distribute_local_to_global(solution_values, global_vector)`. Note
-   * that for multi-component systems where only some of the components are
-   * selected by the present class, the entries in `solution_values` not touched
-   * by this class will be set to zero.
-   *
-   * @param[in] integration_flags Flags specifying which quantities should be
-   * integrated at the points.
-   *
-   * @param[in] sum_into_values Flag specifying if the integrated values
-   * should be summed into the solution values. Defaults to false.
-   *
-   */
-  void
-  test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
-               const EvaluationFlags::EvaluationFlags &integration_flags,
-               const bool                              sum_into_values = false);
-
-  /**
-   * Evaluate values and gradients in face for the selected face (lane) of the
-   * batch. Default stride into the face dofs is width of
-   * VectorizedArray<selected_floating_point_type> which is the default
-   * vectorization over faces for FEFaceEvaluation.
-   */
-  template <int stride_face_dof = VectorizedArrayType::size()>
-  void
-  evaluate_in_face(const ScalarNumber                     *face_dof_values,
-                   const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
-  /**
-   * Integrate values and gradients in face for the selected face (lane) of the
-   * batch. Default stride into the face dofs is width of
-   * VectorizedArray<selected_floating_point_type> which is the default
-   * vectorization over faces for FEFaceEvaluation.
-   */
-  template <int stride_face_dof = VectorizedArrayType::size()>
-  void
-  integrate_in_face(ScalarNumber                           *face_dof_values,
-                    const EvaluationFlags::EvaluationFlags &integration_flags,
-                    const bool sum_into_values = false);
-
-  /**
-   * Return the normal vector. This class or the MappingInfo object passed to
-   * this function needs to be constructed with UpdateFlags containing
-   * `update_normal_vectors`.
-   */
-  Tensor<1, spacedim, Number>
-  normal_vector(const unsigned int point_index) const;
-
-private:
-  static constexpr std::size_t n_lanes_user_interface =
-    internal::VectorizedArrayTrait<Number>::width();
-  static constexpr std::size_t n_lanes_internal =
-    internal::VectorizedArrayTrait<VectorizedArrayType>::width();
-  static constexpr std::size_t stride =
-    internal::VectorizedArrayTrait<Number>::stride();
-
-  template <bool is_linear, std::size_t stride_view>
-  void
-  do_evaluate(
-    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
-
-  template <bool do_JxW, bool is_linear, std::size_t stride_view>
-  void
-  do_integrate(
-    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
-    const EvaluationFlags::EvaluationFlags            &integration_flags,
-    const bool                                         sum_into_values);
-
-  /**
-   * Actually does the evaluation templated on the chosen code path (linear or
-   * higher order).
-   */
-  template <bool is_linear, int stride_face_dof>
-  void
-  do_evaluate_in_face(const ScalarNumber                     *face_dof_values,
-                      const EvaluationFlags::EvaluationFlags &evaluation_flags);
-
-  /**
-   * Actually does the integration templated on the chosen code path (linear or
-   * higher order).
-   */
-  template <bool do_JxW, bool is_linear, int stride_face_dof>
-  void
-  do_integrate_in_face(
-    ScalarNumber                           *face_dof_values,
-    const EvaluationFlags::EvaluationFlags &integration_flags,
-    const bool                              sum_into_values);
-};
-
-
-
 template <int n_components_, int dim, int spacedim, typename Number>
 FEFacePointEvaluation<n_components_, dim, spacedim, Number>::
   FEFacePointEvaluation(
index 8663c0f4f6576300dba51492bc016f4bfb97634e..916e62848f84c3ef23881dd670e25a3735480d31 100644 (file)
@@ -1577,8 +1577,10 @@ namespace NonMatching
         face_number.emplace_back(f_m, f_p);
 
         Assert(
-          cell_m->combined_face_orientation(f_m) == 1 &&
-            cell_p->combined_face_orientation(f_p) == 1,
+          cell_m->combined_face_orientation(f_m) ==
+              ReferenceCell::default_combined_face_orientation() &&
+            cell_p->combined_face_orientation(f_p) ==
+              ReferenceCell::default_combined_face_orientation(),
           ExcMessage(
             "Non standard face orientation is currently not implemented."));
 
index bee712570ee0d1149fe799bb759862768fd4c3b0..3c4297aa5159b0c6a23a985078d5bac5bb0b1b0f 100644 (file)
@@ -20,8 +20,6 @@
 
 #include <deal.II/base/timer.h>
 
-#include <deal.II/distributed/tria.h>
-
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_renumbering.h>
 
@@ -84,7 +82,7 @@ test_dg_fcl(const unsigned int degree, const bool curved_mesh)
 
   const unsigned int n_q_points = degree + 1;
 
-  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  Triangulation<dim> tria;
 
   if (curved_mesh && dim > 1)
     GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1, 6);
index f7dc976cda23f794c368d498a25b7367766abe13..8ddcc184e972b540b93487a697f5292458a9f72b 100644 (file)
@@ -20,8 +20,6 @@
 
 #include <deal.II/base/timer.h>
 
-#include <deal.II/distributed/tria.h>
-
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_renumbering.h>
 
@@ -80,7 +78,7 @@ test_dg_ecl(const unsigned int degree, const bool curved_mesh)
 
   const unsigned int n_q_points = degree + 1;
 
-  parallel::distributed::Triangulation<dim> tria(MPI_COMM_WORLD);
+  Triangulation<dim> tria;
 
   if (curved_mesh && dim > 1)
     GridGenerator::hyper_shell(tria, Point<dim>(), 0.5, 1, 6);
@@ -318,7 +316,7 @@ test_dg_ecl(const unsigned int degree, const bool curved_mesh)
                     if (mask[v] == false)
                       {
                         for (unsigned int i = 0;
-                             i < 2 * fe_eval_m.dofs_per_face;
+                             i < 2 * fe_eval_m.get_dofs_projected_to_face();
                              ++i)
                           fe_eval_m.get_scratch_data().begin()[i][v] = 0.;
                         continue;

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.