]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Optimize FEPointEvaluation/NM::MappingInfo 15971/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Mon, 28 Aug 2023 11:20:46 +0000 (13:20 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Tue, 19 Sep 2023 08:55:50 +0000 (10:55 +0200)
include/deal.II/base/derivative_form.h
include/deal.II/base/floating_point_comparator.h
include/deal.II/matrix_free/fe_point_evaluation.h
include/deal.II/matrix_free/tensor_product_kernels.h
include/deal.II/non_matching/mapping_info.h

index e3a73fa695077aaa14f1ff3a50b8ae68c52eff57..8c1e4710f689211d28a32339041d581441195453 100644 (file)
@@ -591,6 +591,82 @@ transpose(const DerivativeForm<1, dim, spacedim, Number> &DF)
 }
 
 
+
+/**
+ * Specialization of apply_transformation() for a diagonal DerivativeForm.
+ *
+ * @relatesalso DerivativeForm
+ */
+template <int spacedim, int dim, typename Number1, typename Number2>
+inline Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>
+apply_diagonal_transformation(
+  const DerivativeForm<1, dim, spacedim, Number1> &grad_F,
+  const Tensor<1, dim, Number2>                   &d_x)
+{
+  Assert(dim == spacedim,
+         ExcMessage("Only dim = spacedim allowed for diagonal transformation"));
+  Tensor<1, spacedim, typename ProductType<Number1, Number2>::type> dest;
+  for (unsigned int i = 0; i < spacedim; ++i)
+    dest[i] = grad_F[i][i] * d_x[i];
+  return dest;
+}
+
+
+/**
+ * Similar to the previous apply_diagonal_transformation(), specialized for the
+ * case `dim == spacedim` where we can return a rank-2 tensor instead of the
+ * more general `DerivativeForm`.
+ * Each row of the result corresponds to one of the rows of @p D_X transformed
+ * by @p grad_F, equivalent to $\mathrm{D\_X} \, \mathrm{grad\_F}^T$ in matrix
+ * notation.
+ *
+ * @relatesalso DerivativeForm
+ */
+template <int dim, typename Number1, typename Number2>
+inline Tensor<2, dim, typename ProductType<Number1, Number2>::type>
+apply_diagonal_transformation(
+  const DerivativeForm<1, dim, dim, Number1> &grad_F,
+  const Tensor<2, dim, Number2>              &D_X)
+{
+  Tensor<2, dim, typename ProductType<Number1, Number2>::type> dest;
+  for (unsigned int i = 0; i < dim; ++i)
+    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
+
+  return dest;
+}
+
+
+
+/**
+ * Similar to the previous apply_diagonal_transformation().
+ * Each row of the result corresponds to one of the rows of @p D_X transformed
+ * by @p grad_F.
+ *
+ * @relatesalso DerivativeForm
+ */
+template <int spacedim,
+          int dim,
+          int n_components,
+          typename Number1,
+          typename Number2>
+inline Tensor<1,
+              n_components,
+              Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
+apply_diagonal_transformation(
+  const DerivativeForm<1, dim, spacedim, Number1>        &grad_F,
+  const Tensor<1, n_components, Tensor<1, dim, Number2>> &D_X)
+{
+  Tensor<1,
+         n_components,
+         Tensor<1, spacedim, typename ProductType<Number1, Number2>::type>>
+    dest;
+  for (unsigned int i = 0; i < n_components; ++i)
+    dest[i] = apply_diagonal_transformation(grad_F, D_X[i]);
+
+  return dest;
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index 9c3e1214d9680913407580ea75ae85a089943e6e..fc76666ab3c1cab50900c03f791c77cf8a63da04 100644 (file)
@@ -18,6 +18,7 @@
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/derivative_form.h>
 #include <deal.II/base/table.h>
 #include <deal.II/base/tensor.h>
 #include <deal.II/base/vectorization.h>
@@ -105,6 +106,14 @@ struct FloatingPointComparator
   ComparisonResult
   compare(const Tensor<rank, dim, T> &t1, const Tensor<rank, dim, T> &t2) const;
 
+  /**
+   * Compare two derivative forms.
+   */
+  template <int rank, int dim, int spacedim, typename T>
+  ComparisonResult
+  compare(const DerivativeForm<rank, dim, spacedim, T> &t1,
+          const DerivativeForm<rank, dim, spacedim, T> &t2) const;
+
   /**
    * Compare two tables.
    */
@@ -217,6 +226,24 @@ FloatingPointComparator<Number>::compare(const Tensor<rank, dim, T> &t1,
 
 
 
+template <typename Number>
+template <int rank, int dim, int spacedim, typename T>
+typename FloatingPointComparator<Number>::ComparisonResult
+FloatingPointComparator<Number>::compare(
+  const DerivativeForm<rank, dim, spacedim, T> &t1,
+  const DerivativeForm<rank, dim, spacedim, T> &t2) const
+{
+  for (unsigned int i = 0; i < dim; ++i)
+    {
+      const ComparisonResult result = compare(t1[i], t2[i]);
+      if (result != ComparisonResult::equal)
+        return result;
+    }
+  return ComparisonResult::equal;
+}
+
+
+
 template <typename Number>
 template <typename T>
 typename FloatingPointComparator<Number>::ComparisonResult
index 9c7a56854236f303c97bc040c0288abd8c46e559..9189fd0b2e38e477a95508f0444ddfead4c1598b 100644 (file)
@@ -30,6 +30,7 @@
 
 #include <deal.II/matrix_free/evaluation_flags.h>
 #include <deal.II/matrix_free/evaluation_kernels_face.h>
+#include <deal.II/matrix_free/mapping_info_storage.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
@@ -85,13 +86,12 @@ namespace internal
         result[component] = vector_entry;
       }
 
-      static void
-      write_value(VectorizedArrayType         &vector_entry,
-                  const unsigned int           component,
-                  const vectorized_value_type &result)
+      static ScalarNumber
+      sum_value(const unsigned int           component,
+                const vectorized_value_type &result)
       {
         AssertIndexRange(component, n_components);
-        vector_entry = result[component];
+        return result[component].sum();
       }
 
       static void
@@ -237,12 +237,10 @@ namespace internal
         result = vector_entry;
       }
 
-      static void
-      write_value(VectorizedArrayType &vector_entry,
-                  const unsigned int,
-                  const vectorized_value_type &result)
+      static ScalarNumber
+      sum_value(const unsigned int, const vectorized_value_type &result)
       {
-        vector_entry = result;
+        return result.sum();
       }
 
       static void
@@ -389,15 +387,16 @@ namespace internal
                  const unsigned int component,
                  scalar_value_type &result)
       {
+        AssertIndexRange(component, dim);
         result[component] = vector_entry;
       }
 
-      static void
-      write_value(VectorizedArrayType         &vector_entry,
-                  const unsigned int           component,
-                  const vectorized_value_type &result)
+      static ScalarNumber
+      sum_value(const unsigned int           component,
+                const vectorized_value_type &result)
       {
-        vector_entry = result[component];
+        AssertIndexRange(component, dim);
+        return result[component].sum();
       }
 
       static void
@@ -543,12 +542,10 @@ namespace internal
         result = vector_entry;
       }
 
-      static void
-      write_value(VectorizedArrayType &vector_entry,
-                  const unsigned int,
-                  const vectorized_value_type &result)
+      static ScalarNumber
+      sum_value(const unsigned int, const vectorized_value_type &result)
       {
-        vector_entry = result;
+        return result.sum();
       }
 
       static void
@@ -834,6 +831,23 @@ public:
   void
   reinit(const unsigned int cell_index, const unsigned int face_number);
 
+  /**
+   * 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().
@@ -847,7 +861,36 @@ public:
    */
   void
   evaluate(const ArrayView<const ScalarNumber>    &solution_values,
-           const EvaluationFlags::EvaluationFlags &evaluation_flags);
+           const EvaluationFlags::EvaluationFlags &evaluation_flags)
+  {
+    evaluate(StridedArrayView<const ScalarNumber, 1>(solution_values.data(),
+                                                     solution_values.size()),
+             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.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  integrate(const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+            const EvaluationFlags::EvaluationFlags &integration_flags);
+
 
   /**
    * This function multiplies the quantities passed in by previous
@@ -869,7 +912,12 @@ public:
    */
   void
   integrate(const ArrayView<ScalarNumber>          &solution_values,
-            const EvaluationFlags::EvaluationFlags &integration_flags);
+            const EvaluationFlags::EvaluationFlags &integration_flags)
+  {
+    integrate(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+                                                solution_values.size()),
+              integration_flags);
+  }
 
   /**
    * This function multiplies the quantities passed in by previous
@@ -882,7 +930,35 @@ public:
    * finite element formulation.
    *
    * @param[out] solution_values This array will contain the result of the
-   * integral, which can be used to during
+   * 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.
+   *
+   */
+  template <std::size_t stride_view>
+  void
+  test_and_sum(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_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. 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
@@ -895,7 +971,12 @@ public:
    */
   void
   test_and_sum(const ArrayView<ScalarNumber>          &solution_values,
-               const EvaluationFlags::EvaluationFlags &integration_flags);
+               const EvaluationFlags::EvaluationFlags &integration_flags)
+  {
+    test_and_sum(StridedArrayView<ScalarNumber, 1>(solution_values.data(),
+                                                   solution_values.size()),
+                 integration_flags);
+  }
 
   /**
    * Return the value at quadrature point number @p point_index after a call to
@@ -930,18 +1011,6 @@ public:
   const gradient_type &
   get_gradient(const unsigned int point_index) const;
 
-  /**
-   * Return the gradient in unit coordinates at the point with index
-   * `point_index` after a call to FEPointEvaluation::evaluate() with
-   * EvaluationFlags::gradients set, or the gradient that has been stored there
-   * with a call to FEPointEvaluation::submit_gradient(). If the object is
-   * vector-valued, a vector-valued return argument is given. Note that when
-   * vectorization is enabled, values from several points are grouped
-   * together.
-   */
-  const gradient_type &
-  get_unit_gradient(const unsigned int point_index) const;
-
   /**
    * Write a contribution that is tested by the gradient to the field
    * containing the values on points with the given `point_index`. Access to
@@ -1032,6 +1101,7 @@ private:
    * precomputes the @p shapes vector, holding the evaluation of 1D basis
    * functions of tensor product polynomials, if necessary.
    */
+  template <bool is_face, bool is_linear>
   void
   do_reinit();
 
@@ -1039,87 +1109,109 @@ private:
    * Resizes necessary data fields, reads in and renumbers solution values.
    * Interpolates onto face if face path is selected.
    */
+  template <bool is_face_path, bool is_linear, std::size_t stride_view>
   void
   prepare_evaluate_fast(
-    const ArrayView<const ScalarNumber>    &solution_values,
-    const EvaluationFlags::EvaluationFlags &evaluation_flags);
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
   /**
    * Evaluates the actual interpolation on the cell or face for a quadrature
    * batch.
    */
+  template <bool is_face_path, bool is_linear, std::size_t stride_view>
   void
   compute_evaluate_fast(
-    const EvaluationFlags::EvaluationFlags &evaluation_flags,
-    const unsigned int                      n_shapes,
-    const unsigned int                      qb,
-    vectorized_value_type                  &value,
-    interface_vectorized_gradient_type     &gradient);
+    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);
 
   /**
    * Fast path of the evaluate function.
    */
+  template <bool is_face_path, bool is_linear, std::size_t stride_view>
   void
-  evaluate_fast(const ArrayView<const ScalarNumber>    &solution_values,
-                const EvaluationFlags::EvaluationFlags &evaluation_flags);
+  evaluate_fast(
+    const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags                  &evaluation_flags);
 
   /**
    * Slow path of the evaluate function using FEValues.
    */
+  template <std::size_t stride_view>
   void
-  evaluate_slow(const ArrayView<const ScalarNumber>    &solution_values,
-                const EvaluationFlags::EvaluationFlags &evaluation_flags);
+  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_face_path, 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);
+    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);
 
   /**
    * 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_face_path, bool is_linear, std::size_t stride_view>
   void
   finish_integrate_fast(
-    const ArrayView<ScalarNumber>          &solution_values,
-    const EvaluationFlags::EvaluationFlags &integration_flags);
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags,
+    vectorized_value_type *solution_values_vectorized_linear);
 
   /**
    * Fast path of the integrate function.
    */
-  template <bool do_JxW>
+  template <bool        do_JxW,
+            bool        is_face_path,
+            bool        is_linear,
+            std::size_t stride_view>
   void
-  integrate_fast(const ArrayView<ScalarNumber>          &solution_values,
-                 const EvaluationFlags::EvaluationFlags &integration_flags);
+  integrate_fast(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags);
 
   /**
    * Slow path of the integrate function using FEValues.
    */
-  template <bool do_JxW>
+  template <bool do_JxW, std::size_t stride_view>
   void
-  integrate_slow(const ArrayView<ScalarNumber>          &solution_values,
-                 const EvaluationFlags::EvaluationFlags &integration_flags);
+  integrate_slow(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags);
 
   /**
    * Implementation of the integrate/test_and_sum function.
    */
-  template <bool do_JxW>
+  template <bool do_JxW, std::size_t stride_view>
   void
-  do_integrate(const ArrayView<ScalarNumber>          &solution_values,
-               const EvaluationFlags::EvaluationFlags &integration_flags);
+  do_integrate(
+    const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+    const EvaluationFlags::EvaluationFlags            &integration_flags);
 
   /**
    * Number of quadrature batches of the current cell/face.
    */
+  const unsigned int n_q_batches;
+
+  /**
+   * Number of quadrature points/batches of the current cell/face.
+   */
   const unsigned int n_q_points;
 
   /**
@@ -1144,9 +1236,9 @@ private:
   std::vector<Polynomials::Polynomial<double>> poly;
 
   /**
-   * Store whether the polynomials are linear with nodes at 0 and 1.
+   * Store whether the linear path should be used.
    */
-  bool polynomials_are_hat_functions;
+  bool use_linear_path;
 
   /**
    * Renumbering between the unknowns of unknowns implied by the FiniteElement
@@ -1181,11 +1273,6 @@ private:
    */
   std::vector<value_type> values;
 
-  /**
-   * Temporary array to store the gradients in unit coordinates at the points.
-   */
-  std::vector<gradient_type> unit_gradients;
-
   /**
    * Temporary array to store the gradients in real coordinates at the points.
    */
@@ -1234,6 +1321,11 @@ private:
    */
   const Number *JxW_ptr;
 
+  /**
+   * Cell type describing the geometry of the cell and compression of jacobians.
+   */
+  internal::MatrixFreeFunctions::GeometryType cell_type;
+
   /**
    * Number of unknowns per component, i.e., number of unique basis functions,
    * for the chosen FiniteElement (or base element).
@@ -1340,7 +1432,8 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   const FiniteElement<dim> &fe,
   const UpdateFlags         update_flags,
   const unsigned int        first_selected_component)
-  : n_q_points(numbers::invalid_unsigned_int)
+  : 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)
@@ -1366,7 +1459,8 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
   const FiniteElement<dim>                        &fe,
   const unsigned int                               first_selected_component)
-  : n_q_points(numbers::invalid_unsigned_int)
+  : 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)
@@ -1388,17 +1482,17 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   FEPointEvaluation<n_components_, dim, spacedim, Number> &other) noexcept
-  : n_q_points(other.n_q_points)
+  : 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)
-  , polynomials_are_hat_functions(other.polynomials_are_hat_functions)
+  , 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)
-  , unit_gradients(other.unit_gradients)
   , gradients(other.gradients)
   , dofs_per_component(other.dofs_per_component)
   , dofs_per_component_face(other.dofs_per_component_face)
@@ -1430,17 +1524,17 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
 template <int n_components_, int dim, int spacedim, typename Number>
 FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   FEPointEvaluation<n_components_, dim, spacedim, Number> &&other) noexcept
-  : n_q_points(other.n_q_points)
+  : 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)
-  , polynomials_are_hat_functions(other.polynomials_are_hat_functions)
+  , 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)
-  , unit_gradients(other.unit_gradients)
   , gradients(other.gradients)
   , dofs_per_component(other.dofs_per_component)
   , dofs_per_component_face(other.dofs_per_component_face)
@@ -1519,15 +1613,18 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::setup(
       if (is_lexicographic)
         renumber.clear();
 
-      polynomials_are_hat_functions =
-        (poly.size() == 2 && poly[0].value(0.) == 1. &&
-         poly[0].value(1.) == 0. && poly[1].value(0.) == 0. &&
-         poly[1].value(1.) == 1.);
+      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);
 
       const unsigned int size_face = 2 * dofs_per_component_face;
       const unsigned int size_cell = dofs_per_component;
       scratch_data_scalar.resize(size_face + size_cell);
 
+      solution_renumbered.resize(dofs_per_component);
+      solution_renumbered_vectorized.resize(dofs_per_component);
+
       fast_path = true;
     }
   else
@@ -1556,7 +1653,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::setup(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-void
+inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
   const typename Triangulation<dim, spacedim>::cell_iterator &cell,
   const ArrayView<const Point<dim>>                          &unit_points)
@@ -1577,48 +1674,50 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
       fe_values->reinit(cell);
     }
 
-  do_reinit();
+  if (use_linear_path)
+    do_reinit<false, true>();
+  else
+    do_reinit<false, false>();
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-void
+inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit()
 {
   current_cell_index  = numbers::invalid_unsigned_int;
   current_face_number = numbers::invalid_unsigned_int;
 
-  do_reinit();
+  if (use_linear_path)
+    do_reinit<false, true>();
+  else
+    do_reinit<false, false>();
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-void
+inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
   const unsigned int cell_index)
 {
   current_cell_index  = cell_index;
   current_face_number = numbers::invalid_unsigned_int;
 
-  do_reinit();
+  if (use_linear_path)
+    do_reinit<false, true>();
+  else
+    do_reinit<false, false>();
 
   if (!fast_path)
     {
-      const auto unit_points_vectorized = mapping_info->get_unit_point(
-        mapping_info->compute_unit_point_index_offset(current_cell_index,
-                                                      current_face_number));
-      const unsigned int n_q_points_unvectorized =
-        mapping_info->get_n_q_points_unvectorized(current_cell_index,
-                                                  current_face_number);
-
-      std::vector<Point<dim>> unit_points(n_q_points_unvectorized);
+      std::vector<Point<dim>> unit_points(n_q_points_scalar);
 
-      for (unsigned int v = 0; v < n_q_points_unvectorized; ++v)
+      for (unsigned int v = 0; v < n_q_points_scalar; ++v)
         for (unsigned int d = 0; d < dim; ++d)
-          unit_points[v][d] = unit_points_vectorized[v / n_lanes_internal][d]
-                                                    [v % n_lanes_internal];
+          unit_points[v][d] =
+            unit_point_ptr[v / n_lanes_internal][d][v % n_lanes_internal];
 
       fe_values = std::make_shared<FEValues<dim, spacedim>>(
         *mapping,
@@ -1634,7 +1733,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-void
+inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
   const unsigned int cell_index,
   const unsigned int face_number)
@@ -1642,27 +1741,44 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::reinit(
   current_cell_index  = cell_index;
   current_face_number = face_number;
 
-  do_reinit();
+  if (use_linear_path)
+    do_reinit<true, true>();
+  else
+    do_reinit<true, false>();
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face, bool is_linear>
 inline void
 FEPointEvaluation<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);
+
   const_cast<unsigned int &>(n_q_points_scalar) =
-    mapping_info->get_n_q_points_unvectorized(current_cell_index,
-                                              current_face_number);
+    mapping_info->get_n_q_points_unvectorized(geometry_index);
+
+  // 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;
+
+  const unsigned int n_q_points_before = n_q_points;
 
-  // round up n_points_scalar / n_lanes_user_interface
   const_cast<unsigned int &>(n_q_points) =
-    (n_q_points_scalar + n_lanes_user_interface - 1) / n_lanes_user_interface;
+    (stride == 1) ? n_q_batches : n_q_points_scalar;
 
-  if (update_flags & update_values)
-    values.resize(n_q_points, numbers::signaling_nan<value_type>());
-  if (update_flags & update_gradients)
-    gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
+  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);
+    }
 
   if (n_q_points == 0)
     {
@@ -1672,12 +1788,12 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_reinit()
 
   // use face path if mapping_info in face state and number of quadrature points
   // is large enough
-  use_face_path = mapping_info->is_face_state() && n_q_points_scalar >= 6;
+  use_face_path =
+    is_face && (mapping_info->is_face_state() && n_q_points_scalar >= 6);
 
   // set unit point pointer
   const unsigned int unit_point_offset =
-    mapping_info->compute_unit_point_index_offset(current_cell_index,
-                                                  current_face_number);
+    mapping_info->compute_unit_point_index_offset(geometry_index);
 
   if (use_face_path)
     unit_point_faces_ptr =
@@ -1686,35 +1802,43 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_reinit()
     unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
 
   // 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();
-  const unsigned int data_offset =
-    mapping_info->compute_data_index_offset(current_cell_index,
-                                            current_face_number);
   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(data_offset);
+    jacobian_ptr = mapping_info->get_jacobian(compressed_data_offset);
   if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
-    inverse_jacobian_ptr = mapping_info->get_inverse_jacobian(data_offset);
-  if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+    inverse_jacobian_ptr =
+      mapping_info->get_inverse_jacobian(compressed_data_offset);
+  if (is_face && 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);
+  inverse_jacobian_ptr =
+    mapping_info->get_inverse_jacobian(compressed_data_offset);
+  normal_ptr = mapping_info->get_normal_vector(data_offset);
+  JxW_ptr    = mapping_info->get_JxW(data_offset);
+#endif
 
-  if (fast_path && !polynomials_are_hat_functions)
+  if (!is_linear && fast_path)
     {
-      // round up n_q_points_scalar / n_lanes_internal
-      const std::size_t n_batches =
-        (n_q_points_scalar + n_lanes_internal - 1) / n_lanes_internal;
       const std::size_t n_shapes = poly.size();
 
-      for (unsigned int qb = 0; qb < n_batches; ++qb)
+      for (unsigned int qb = 0; qb < n_q_batches; ++qb)
         if (use_face_path)
           {
             if (dim > 1)
               {
-                shapes_faces.resize_fast(n_batches * n_shapes);
+                shapes_faces.resize_fast(n_q_batches * n_shapes);
                 internal::compute_values_of_array(
                   shapes_faces.data() + qb * n_shapes,
                   poly,
@@ -1724,7 +1848,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_reinit()
           }
         else
           {
-            shapes.resize_fast(n_batches * n_shapes);
+            shapes.resize_fast(n_q_batches * n_shapes);
             internal::compute_values_of_array(
               shapes.data() + qb * n_shapes,
               poly,
@@ -1739,45 +1863,38 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_reinit()
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face_path, bool is_linear, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::prepare_evaluate_fast(
-  const ArrayView<const ScalarNumber>    &solution_values,
-  const EvaluationFlags::EvaluationFlags &evaluation_flags)
+  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
 {
-  if (use_face_path)
-    {
-      if (solution_renumbered.size() != 2 * dofs_per_component_face)
-        solution_renumbered.resize(2 * dofs_per_component_face);
-    }
-  else
-    {
-      if (solution_renumbered.size() != dofs_per_component)
-        solution_renumbered.resize(dofs_per_component);
-    }
+  const unsigned int dofs_per_comp =
+    is_linear ? Utilities::pow(2, dim) : dofs_per_component;
+
   for (unsigned int comp = 0; comp < n_components; ++comp)
     {
       const std::size_t offset =
-        (component_in_base_element + comp) * dofs_per_component;
+        (component_in_base_element + comp) * dofs_per_comp;
 
-      if (use_face_path)
+      if (is_face_path)
         {
           const ScalarNumber *input;
-          if (renumber.empty())
+          if (is_linear || renumber.empty())
             {
-              for (unsigned int i = 0; i < dofs_per_component; ++i)
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
                 scratch_data_scalar[i] = solution_values[i + offset];
               input = scratch_data_scalar.data();
             }
           else
             {
               const unsigned int *renumber_ptr = renumber.data() + offset;
-              for (unsigned int i = 0; i < dofs_per_component; ++i)
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
                 scratch_data_scalar[i] = solution_values[renumber_ptr[i]];
               input = scratch_data_scalar.data();
             }
 
-          ScalarNumber *output =
-            scratch_data_scalar.begin() + dofs_per_component;
+          ScalarNumber *output = scratch_data_scalar.begin() + dofs_per_comp;
 
           internal::FEFaceNormalEvaluationImpl<dim, -1, ScalarNumber>::
             template interpolate<true, false>(1,
@@ -1787,14 +1904,16 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::prepare_evaluate_fast(
                                               output,
                                               current_face_number);
 
-          for (unsigned int i = 0; i < 2 * dofs_per_component_face; ++i)
+          const unsigned int dofs_per_comp_face =
+            is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
+          for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
             ETT::read_value(output[i], comp, solution_renumbered[i]);
         }
       else
         {
-          if (renumber.empty())
+          if ((is_linear && n_components == 1) || renumber.empty())
             {
-              for (unsigned int i = 0; i < dofs_per_component; ++i)
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
                 ETT::read_value(solution_values[i + offset],
                                 comp,
                                 solution_renumbered[i]);
@@ -1802,43 +1921,39 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::prepare_evaluate_fast(
           else
             {
               const unsigned int *renumber_ptr = renumber.data() + offset;
-              for (unsigned int i = 0; i < dofs_per_component; ++i)
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
                 ETT::read_value(solution_values[renumber_ptr[i]],
                                 comp,
                                 solution_renumbered[i]);
             }
         }
     }
-
-  // unit gradients are currently only implemented with the fast tensor
-  // path
-  unit_gradients.resize(n_q_points, numbers::signaling_nan<gradient_type>());
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face_path, bool is_linear, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
-  const EvaluationFlags::EvaluationFlags &evaluation_flags,
-  const unsigned int                      n_shapes,
-  const unsigned int                      qb,
-  vectorized_value_type                  &value,
-  interface_vectorized_gradient_type     &gradient)
+  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 (use_face_path)
+  if (is_face_path)
     {
       if (evaluation_flags & EvaluationFlags::gradients)
         {
           const std::array<vectorized_value_type, dim + 1> interpolated_value =
-            polynomials_are_hat_functions ?
+            is_linear ?
               internal::evaluate_tensor_product_value_and_gradient_linear<
                 dim - 1,
                 scalar_value_type,
                 VectorizedArrayType,
-                2>(n_shapes,
-                   solution_renumbered.data(),
-                   unit_point_faces_ptr[qb]) :
+                2>(solution_renumbered.data(), unit_point_faces_ptr[qb]) :
               internal::evaluate_tensor_product_value_and_gradient_shapes<
                 dim - 1,
                 scalar_value_type,
@@ -1889,30 +2004,42 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
         }
       else
         {
-          value = polynomials_are_hat_functions ?
-                    internal::evaluate_tensor_product_value_linear<
-                      dim - 1,
-                      scalar_value_type,
-                      VectorizedArrayType>(n_shapes,
-                                           solution_renumbered.data(),
-                                           unit_point_faces_ptr[qb]) :
-                    internal::evaluate_tensor_product_value_shapes<
-                      dim - 1,
-                      scalar_value_type,
-                      VectorizedArrayType,
-                      false>(shapes_faces.data() + qb * n_shapes,
-                             n_shapes,
-                             solution_renumbered.data());
+          value = is_linear ? internal::evaluate_tensor_product_value_linear<
+                                dim - 1,
+                                scalar_value_type,
+                                VectorizedArrayType>(solution_renumbered.data(),
+                                                     unit_point_faces_ptr[qb]) :
+                              internal::evaluate_tensor_product_value_shapes<
+                                dim - 1,
+                                scalar_value_type,
+                                VectorizedArrayType,
+                                false>(shapes_faces.data() + qb * n_shapes,
+                                       n_shapes,
+                                       solution_renumbered.data());
         }
     }
   else
     {
       if (evaluation_flags & EvaluationFlags::gradients)
         {
-          const std::array<vectorized_value_type, dim + 1> result =
-            polynomials_are_hat_functions ?
-              internal::evaluate_tensor_product_value_and_gradient_linear(
-                n_shapes, solution_renumbered.data(), unit_point_ptr[qb]) :
+          std::array<vectorized_value_type, dim + 1> result;
+          if constexpr (is_linear)
+            {
+              if constexpr (n_components == 1)
+                result =
+                  internal::evaluate_tensor_product_value_and_gradient_linear<
+                    dim,
+                    scalar_value_type,
+                    VectorizedArrayType,
+                    1,
+                    stride_view>(solution_values.data(), unit_point_ptr[qb]);
+              else
+                result =
+                  internal::evaluate_tensor_product_value_and_gradient_linear(
+                    solution_renumbered.data(), unit_point_ptr[qb]);
+            }
+          else
+            result =
               internal::evaluate_tensor_product_value_and_gradient_shapes<
                 dim,
                 scalar_value_type,
@@ -1930,17 +2057,26 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
         }
       else
         {
-          value =
-            polynomials_are_hat_functions ?
-              internal::evaluate_tensor_product_value_linear(
-                n_shapes, solution_renumbered.data(), unit_point_ptr[qb]) :
-              internal::evaluate_tensor_product_value_shapes<
-                dim,
-                scalar_value_type,
-                VectorizedArrayType,
-                false>(shapes.data() + qb * n_shapes,
-                       n_shapes,
-                       solution_renumbered.data());
+          if constexpr (is_linear)
+            {
+              if constexpr (n_components == 1)
+                value = internal::evaluate_tensor_product_value_linear<
+                  dim,
+                  scalar_value_type,
+                  VectorizedArrayType,
+                  stride_view>(solution_values.data(), unit_point_ptr[qb]);
+              else
+                value = internal::evaluate_tensor_product_value_linear(
+                  solution_renumbered.data(), unit_point_ptr[qb]);
+            }
+          else
+            value = internal::evaluate_tensor_product_value_shapes<
+              dim,
+              scalar_value_type,
+              VectorizedArrayType,
+              false>(shapes.data() + qb * n_shapes,
+                     n_shapes,
+                     solution_renumbered.data());
         }
     }
 }
@@ -1948,28 +2084,33 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_evaluate_fast(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face_path, bool is_linear, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
-  const ArrayView<const ScalarNumber>    &solution_values,
-  const EvaluationFlags::EvaluationFlags &evaluation_flags)
+  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
 {
-  prepare_evaluate_fast(solution_values, evaluation_flags);
-
-  // loop over quadrature batches qb / points q
-  const unsigned int                 n_shapes = poly.size();
-  vectorized_value_type              value;
-  interface_vectorized_gradient_type gradient;
-  for (unsigned int qb = 0, q = 0; q < n_q_points_scalar;
-       ++qb, q += n_lanes_internal)
+  if (!(is_linear && n_components == 1) || is_face_path)
+    prepare_evaluate_fast<is_face_path, is_linear>(solution_values,
+                                                   evaluation_flags);
+
+  // loop over quadrature batches qb
+  const unsigned int n_shapes = is_linear ? 2 : poly.size();
+
+  for (unsigned int qb = 0; qb < n_q_batches; ++qb)
     {
-      compute_evaluate_fast(evaluation_flags, n_shapes, qb, value, gradient);
+      vectorized_value_type              value;
+      interface_vectorized_gradient_type gradient;
+
+      compute_evaluate_fast<is_face_path, is_linear>(
+        solution_values, evaluation_flags, n_shapes, qb, value, gradient);
 
       if (evaluation_flags & EvaluationFlags::values)
         {
-          for (unsigned int v = 0;
-               v < stride && (stride == 1 || q + v < n_q_points_scalar);
-               ++v)
-            ETT::set_value(value, v, values[qb * stride + v]);
+          for (unsigned int v = 0, offset = qb * stride;
+               v < stride && (stride == 1 || offset < n_q_points_scalar);
+               ++v, ++offset)
+            ETT::set_value(value, v, values[offset]);
         }
       if (evaluation_flags & EvaluationFlags::gradients)
         {
@@ -1977,15 +2118,25 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
                    update_flags & update_inverse_jacobians,
                  ExcNotInitialized());
 
-          for (unsigned int v = 0;
-               v < stride && (stride == 1 || q + v < n_q_points_scalar);
-               ++v)
+          for (unsigned int v = 0, offset = qb * stride;
+               v < stride && (stride == 1 || offset < n_q_points_scalar);
+               ++v, ++offset)
             {
-              const unsigned int offset = qb * stride + v;
-              ETT::set_gradient(gradient, v, unit_gradients[offset]);
+              gradient_type unit_gradient;
+              ETT::set_gradient(gradient, v, unit_gradient);
               gradients[offset] =
-                apply_transformation(inverse_jacobian_ptr[offset].transpose(),
-                                     unit_gradients[offset]);
+                cell_type <=
+                    internal::MatrixFreeFunctions::GeometryType::cartesian ?
+                  apply_diagonal_transformation(inverse_jacobian_ptr[0],
+                                                unit_gradient) :
+                  apply_transformation(
+                    inverse_jacobian_ptr
+                      [cell_type <=
+                           internal::MatrixFreeFunctions::GeometryType::affine ?
+                         0 :
+                         offset]
+                        .transpose(),
+                    unit_gradient);
             }
         }
     }
@@ -1994,10 +2145,11 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_fast(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_slow(
-  const ArrayView<const ScalarNumber>    &solution_values,
-  const EvaluationFlags::EvaluationFlags &evaluation_flags)
+  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
 {
   // slow path with FEValues
   Assert(fe_values.get() != nullptr,
@@ -2076,10 +2228,11 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate_slow(
 
 
 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 ArrayView<const ScalarNumber>    &solution_values,
-  const EvaluationFlags::EvaluationFlags &evaluation_flags)
+  const StridedArrayView<const ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags                  &evaluation_flags)
 {
   if (!is_reinitialized)
     reinit();
@@ -2095,7 +2248,22 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
 
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
   if (fast_path)
-    evaluate_fast(solution_values, evaluation_flags);
+    {
+      if (use_face_path)
+        {
+          if (use_linear_path)
+            evaluate_fast<true, true>(solution_values, evaluation_flags);
+          else
+            evaluate_fast<true, false>(solution_values, evaluation_flags);
+        }
+      else
+        {
+          if (use_linear_path)
+            evaluate_fast<false, true>(solution_values, evaluation_flags);
+          else
+            evaluate_fast<false, false>(solution_values, evaluation_flags);
+        }
+    }
   else
     evaluate_slow(solution_values, evaluation_flags);
 }
@@ -2103,15 +2271,17 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::evaluate(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face_path, bool is_linear>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::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)
+  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 (use_face_path)
+  if (is_face_path)
     {
       if (integration_flags & EvaluationFlags::gradients)
         {
@@ -2158,6 +2328,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_integrate_fast(
             Assert(false, ExcInternalError());
 
           internal::integrate_tensor_product_value_and_gradient<
+            is_linear,
             dim - 1,
             VectorizedArrayType,
             vectorized_value_type,
@@ -2165,47 +2336,51 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_integrate_fast(
                n_shapes,
                value_face.data(),
                gradient_in_face,
-               solution_renumbered_vectorized.data(),
+               is_linear ? solution_values_vectorized_linear :
+                           solution_renumbered_vectorized.data(),
                unit_point_faces_ptr[qb],
-               polynomials_are_hat_functions,
                qb != 0);
         }
       else
-        internal::integrate_tensor_product_value<dim - 1,
+        internal::integrate_tensor_product_value<is_linear,
+                                                 dim - 1,
                                                  VectorizedArrayType,
                                                  vectorized_value_type>(
           shapes_faces.data() + qb * n_shapes,
           n_shapes,
           value,
-          solution_renumbered_vectorized.data(),
+          is_linear ? solution_values_vectorized_linear :
+                      solution_renumbered_vectorized.data(),
           unit_point_faces_ptr[qb],
-          polynomials_are_hat_functions,
           qb != 0);
     }
   else
     {
       if (integration_flags & EvaluationFlags::gradients)
         internal::integrate_tensor_product_value_and_gradient<
+          is_linear,
           dim,
           VectorizedArrayType,
           vectorized_value_type>(shapes.data() + qb * n_shapes,
                                  n_shapes,
                                  &value,
                                  gradient,
-                                 solution_renumbered_vectorized.data(),
+                                 is_linear ?
+                                   solution_values_vectorized_linear :
+                                   solution_renumbered_vectorized.data(),
                                  unit_point_ptr[qb],
-                                 polynomials_are_hat_functions,
                                  qb != 0);
       else
-        internal::integrate_tensor_product_value<dim,
+        internal::integrate_tensor_product_value<is_linear,
+                                                 dim,
                                                  VectorizedArrayType,
                                                  vectorized_value_type>(
           shapes.data() + qb * n_shapes,
           n_shapes,
           value,
-          solution_renumbered_vectorized.data(),
+          is_linear ? solution_values_vectorized_linear :
+                      solution_renumbered_vectorized.data(),
           unit_point_ptr[qb],
-          polynomials_are_hat_functions,
           qb != 0);
     }
 }
@@ -2213,31 +2388,40 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::compute_integrate_fast(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <bool is_face_path, bool is_linear, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags,
+  vectorized_value_type *solution_values_vectorized_linear)
 {
-  std::fill(solution_values.begin(), solution_values.end(), ScalarNumber());
+  if (fe->n_components() > n_components)
+    for (unsigned int i = 0; i < solution_values.size(); ++i)
+      solution_values[i] = 0;
+
+  const unsigned int dofs_per_comp =
+    is_linear ? Utilities::pow(2, dim) : dofs_per_component;
+
   for (unsigned int comp = 0; comp < n_components; ++comp)
     {
       const std::size_t offset =
-        (component_in_base_element + comp) * dofs_per_component;
+        (component_in_base_element + comp) * dofs_per_comp;
 
-      if (use_face_path)
+      if (is_face_path)
         {
-          const unsigned int size_input = 2 * dofs_per_component_face;
+          const unsigned int dofs_per_comp_face =
+            is_linear ? Utilities::pow(2, dim - 1) : dofs_per_component_face;
+
+          const unsigned int size_input = 2 * dofs_per_comp_face;
           ScalarNumber      *input      = scratch_data_scalar.begin();
           ScalarNumber      *output     = input + size_input;
 
-          for (unsigned int i = 0; i < 2 * dofs_per_component_face; ++i)
-            {
-              VectorizedArrayType vectorized_input;
-              ETT::write_value(vectorized_input,
-                               comp,
+          for (unsigned int i = 0; i < 2 * dofs_per_comp_face; ++i)
+            input[i] =
+              ETT::sum_value(comp,
+                             is_linear ?
+                               *(solution_values_vectorized_linear + i) :
                                solution_renumbered_vectorized[i]);
-              input[i] = vectorized_input.sum();
-            }
 
           internal::FEFaceNormalEvaluationImpl<dim, -1, ScalarNumber>::
             template interpolate<false, false>(1,
@@ -2247,33 +2431,36 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
                                                output,
                                                current_face_number);
 
-          if (renumber.empty())
-            for (unsigned int i = 0; i < dofs_per_component; ++i)
-              solution_values[i + offset] = output[i];
+          if (is_linear || renumber.empty())
+            {
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
+                solution_values[i + offset] = output[i];
+            }
           else
-            for (unsigned int i = 0; i < dofs_per_component; ++i)
-              solution_values[renumber[i + offset]] = output[i];
+            {
+              const unsigned int *renumber_ptr = renumber.data() + offset;
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
+                solution_values[renumber_ptr[i]] = output[i];
+            }
         }
       else
         {
-          if (renumber.empty())
-            for (unsigned int i = 0; i < dofs_per_component; ++i)
-              {
-                VectorizedArrayType result;
-                ETT::write_value(result,
-                                 comp,
-                                 solution_renumbered_vectorized[i]);
-                solution_values[i + offset] = result.sum();
-              }
+          if (is_linear || renumber.empty())
+            {
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
+                solution_values[i + offset] =
+                  ETT::sum_value(comp,
+                                 is_linear ?
+                                   *(solution_values_vectorized_linear + i) :
+                                   solution_renumbered_vectorized[i]);
+            }
           else
-            for (unsigned int i = 0; i < dofs_per_component; ++i)
-              {
-                VectorizedArrayType result;
-                ETT::write_value(result,
-                                 comp,
-                                 solution_renumbered_vectorized[i]);
-                solution_values[renumber[i + offset]] = result.sum();
-              }
+            {
+              const unsigned int *renumber_ptr = renumber.data() + offset;
+              for (unsigned int i = 0; i < dofs_per_comp; ++i)
+                solution_values[renumber_ptr[i]] =
+                  ETT::sum_value(comp, solution_renumbered_vectorized[i]);
+            }
         }
     }
 }
@@ -2281,107 +2468,103 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::finish_integrate_fast(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW>
+template <bool        do_JxW,
+          bool        is_face_path,
+          bool        is_linear,
+          std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_fast(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags)
 {
-  // fast path with tensor product integration
-  if (use_face_path)
-    {
-      if (solution_renumbered_vectorized.size() != 2 * dofs_per_component_face)
-        solution_renumbered_vectorized.resize(2 * dofs_per_component_face);
-    }
-  else
-    {
-      if (solution_renumbered_vectorized.size() != dofs_per_component)
-        solution_renumbered_vectorized.resize(dofs_per_component);
-    }
+  // zero out lanes of incomplete last quadrature point batch
+  if constexpr (stride == 1)
+    if (const unsigned int n_filled_lanes =
+          n_q_points_scalar & (n_lanes_internal - 1);
+        n_filled_lanes > 0)
+      {
+        if (integration_flags & EvaluationFlags::values)
+          for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+            ETT::set_zero_value(values.back(), v);
+        if (integration_flags & EvaluationFlags::gradients)
+          for (unsigned int v = n_filled_lanes; v < n_lanes_internal; ++v)
+            ETT::set_zero_gradient(gradients.back(), v);
+      }
 
-  // loop over quadrature batches qb / points q
-  const unsigned int n_shapes = poly.size();
-  for (unsigned int qb = 0, q = 0; q < n_q_points_scalar;
-       ++qb, q += n_lanes_internal)
-    {
-      const bool incomplete_last_batch =
-        q + n_lanes_user_interface > n_q_points_scalar;
+  std::array<vectorized_value_type,
+             is_linear ? (is_face_path ? 2 * Utilities::pow(2, dim - 1) :
+                                         Utilities::pow(2, dim)) :
+                         0>
+    solution_values_vectorized_linear;
 
+  // loop over quadrature batches qb
+  const unsigned int n_shapes = is_linear ? 2 : poly.size();
+
+  const bool cartesian_cell =
+    cell_type <= internal::MatrixFreeFunctions::GeometryType::cartesian;
+  const bool affine_cell =
+    cell_type <= internal::MatrixFreeFunctions::GeometryType::affine;
+  for (unsigned int qb = 0; qb < n_q_batches; ++qb)
+    {
       vectorized_value_type                 value = {};
       Tensor<1, dim, vectorized_value_type> gradient;
 
       if (integration_flags & EvaluationFlags::values)
-        {
-          // zero out lanes of incomplete last quadrature point batch
-          if (incomplete_last_batch)
-            {
-              const unsigned int n_filled_lanes_last_batch =
-                n_q_points_scalar % n_lanes_internal;
-              for (unsigned int v = n_filled_lanes_last_batch;
-                   v < n_lanes_internal;
-                   ++v)
-                ETT::set_zero_value(values[qb], v);
-            }
+        for (unsigned int v = 0, offset = qb * stride;
+             v < stride && (stride == 1 || offset < n_q_points_scalar);
+             ++v, ++offset)
+          ETT::get_value(value,
+                         v,
+                         do_JxW ? values[offset] * JxW_ptr[offset] :
+                                  values[offset]);
 
-          for (unsigned int v = 0;
-               v < stride && (stride == 1 || q + v < n_q_points_scalar);
-               ++v)
-            {
-              const unsigned int offset = qb * stride + v;
-              if (do_JxW)
-                values[offset] *= JxW_ptr[offset];
-              ETT::get_value(value, v, values[offset]);
-            }
-        }
       if (integration_flags & EvaluationFlags::gradients)
-        {
-          // zero out lanes of incomplete last quadrature point batch
-          if (incomplete_last_batch)
-            {
-              const unsigned int n_filled_lanes_last_batch =
-                n_q_points_scalar % n_lanes_internal;
-              for (unsigned int v = n_filled_lanes_last_batch;
-                   v < n_lanes_internal;
-                   ++v)
-                ETT::set_zero_gradient(gradients[qb], v);
-            }
-
-          for (unsigned int v = 0;
-               v < stride && (stride == 1 || q + v < n_q_points_scalar);
-               ++v)
-            {
-              const unsigned int offset = qb * stride + v;
-              if (do_JxW)
-                gradients[offset] *= JxW_ptr[offset];
-              ETT::get_gradient(
-                gradient,
-                v,
-                apply_transformation(inverse_jacobian_ptr[offset],
-                                     gradients[offset]));
-            }
-        }
+        for (unsigned int v = 0, offset = qb * stride;
+             v < stride && (stride == 1 || offset < n_q_points_scalar);
+             ++v, ++offset)
+          {
+            const auto grad_w =
+              do_JxW ? gradients[offset] * JxW_ptr[offset] : gradients[offset];
+            ETT::get_gradient(
+              gradient,
+              v,
+              cartesian_cell ?
+                apply_diagonal_transformation(inverse_jacobian_ptr[0], grad_w) :
+                apply_transformation(
+                  inverse_jacobian_ptr[affine_cell ? 0 : offset], grad_w));
+          }
 
-      compute_integrate_fast(integration_flags, n_shapes, qb, value, gradient);
+      compute_integrate_fast<is_face_path, is_linear>(
+        integration_flags,
+        n_shapes,
+        qb,
+        value,
+        gradient,
+        solution_values_vectorized_linear.data());
     }
 
   // add between the lanes and write into the result
-  finish_integrate_fast(solution_values, integration_flags);
+  finish_integrate_fast<is_face_path, is_linear>(
+    solution_values,
+    integration_flags,
+    solution_values_vectorized_linear.data());
 }
 
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW>
+template <bool do_JxW, std::size_t stride_view>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags)
 {
   // slow path with FEValues
   Assert(fe_values.get() != nullptr,
          ExcMessage(
            "Not initialized. Please call FEPointEvaluation::reinit()!"));
-  std::fill(solution_values.begin(), solution_values.end(), 0.0);
+  for (unsigned int i = 0; i < solution_values.size(); ++i)
+    solution_values[i] = 0;
 
   const std::size_t n_points = fe_values->get_quadrature().size();
 
@@ -2447,28 +2630,24 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate_slow(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
-template <bool do_JxW>
+template <bool do_JxW, std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags)
 {
   if (!is_reinitialized)
     reinit();
 
-  if (n_q_points == 0) // no evaluation points provided
-    {
-      std::fill(solution_values.begin(), solution_values.end(), 0.0);
-      return;
-    }
-
   Assert(!(integration_flags & EvaluationFlags::hessians), ExcNotImplemented());
 
-  if (!((integration_flags & EvaluationFlags::values) ||
+  if (n_q_points == 0 || // no evaluation points provided
+      !((integration_flags & EvaluationFlags::values) ||
         (integration_flags &
          EvaluationFlags::gradients))) // no integration flags
     {
-      std::fill(solution_values.begin(), solution_values.end(), 0.0);
+      for (unsigned int i = 0; i < solution_values.size(); ++i)
+        solution_values[i] = 0;
       return;
     }
 
@@ -2479,7 +2658,26 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
 
   AssertDimension(solution_values.size(), fe->dofs_per_cell);
   if (fast_path)
-    integrate_fast<do_JxW>(solution_values, integration_flags);
+    {
+      if (use_face_path)
+        {
+          if (use_linear_path)
+            integrate_fast<do_JxW, true, true>(solution_values,
+                                               integration_flags);
+          else
+            integrate_fast<do_JxW, true, false>(solution_values,
+                                                integration_flags);
+        }
+      else
+        {
+          if (use_linear_path)
+            integrate_fast<do_JxW, false, true>(solution_values,
+                                                integration_flags);
+          else
+            integrate_fast<do_JxW, false, false>(solution_values,
+                                                 integration_flags);
+        }
+    }
   else
     integrate_slow<do_JxW>(solution_values, integration_flags);
 }
@@ -2487,10 +2685,11 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::do_integrate(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags)
 {
   do_integrate<true>(solution_values, integration_flags);
 }
@@ -2498,10 +2697,11 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::integrate(
 
 
 template <int n_components_, int dim, int spacedim, typename Number>
+template <std::size_t stride_view>
 void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::test_and_sum(
-  const ArrayView<ScalarNumber>          &solution_values,
-  const EvaluationFlags::EvaluationFlags &integration_flags)
+  const StridedArrayView<ScalarNumber, stride_view> &solution_values,
+  const EvaluationFlags::EvaluationFlags            &integration_flags)
 {
   do_integrate<false>(solution_values, integration_flags);
 }
@@ -2532,22 +2732,6 @@ inline const typename FEPointEvaluation<n_components_, dim, spacedim, Number>::
 
 
 
-template <int n_components_, int dim, int spacedim, typename Number>
-inline const typename FEPointEvaluation<n_components_, dim, spacedim, Number>::
-  gradient_type &
-  FEPointEvaluation<n_components_, dim, spacedim, Number>::get_unit_gradient(
-    const unsigned int point_index) const
-{
-  Assert(fast_path,
-         ExcMessage("Unit gradients are currently only implemented for tensor "
-                    "product finite elements combined with MappingQ "
-                    "mappings"));
-  AssertIndexRange(point_index, unit_gradients.size());
-  return unit_gradients[point_index];
-}
-
-
-
 template <int n_components_, int dim, int spacedim, typename Number>
 inline void
 FEPointEvaluation<n_components_, dim, spacedim, Number>::submit_value(
@@ -2582,7 +2766,10 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::jacobian(
          internal::FEPointEvaluation::
            ExcFEPointEvaluationAccessToUninitializedMappingField(
              "update_jacobians"));
-  return jacobian_ptr[point_index];
+  return jacobian_ptr[cell_type <= ::dealii::internal::MatrixFreeFunctions::
+                                     GeometryType::affine ?
+                        0 :
+                        point_index];
 }
 
 
@@ -2597,7 +2784,11 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::inverse_jacobian(
          internal::FEPointEvaluation::
            ExcFEPointEvaluationAccessToUninitializedMappingField(
              "update_inverse_jacobians"));
-  return inverse_jacobian_ptr[point_index];
+  return inverse_jacobian_ptr
+    [cell_type <=
+         ::dealii::internal::MatrixFreeFunctions::GeometryType::affine ?
+       0 :
+       point_index];
 }
 
 
index 837510b98638dca3901ad6f740773b0b98e97ac3..e6f1fa7b95b90183f1e16ae58404c3b2790ae5b0 100644 (file)
@@ -2272,14 +2272,16 @@ namespace internal
    * Specializes @p evaluate_tensor_product_value_and_gradient() for linear
    * polynomials which massively reduces the necessary instructions.
    */
-  template <int dim, typename Number, typename Number2, int n_values = 1>
+  template <int dim,
+            typename Number,
+            typename Number2,
+            int n_values = 1,
+            int stride   = 1>
   inline std::array<typename ProductTypeNoPoint<Number, Number2>::type,
                     dim + n_values>
   evaluate_tensor_product_value_and_gradient_linear(
-    const unsigned int               n_shapes,
-    const Number                    *values,
-    const Point<dim, Number2>       &p,
-    const std::vector<unsigned int> &renumber = {})
+    const Number              *values,
+    const Point<dim, Number2> &p)
   {
     static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
     static_assert(1 <= n_values && n_values <= 2,
@@ -2287,15 +2289,12 @@ namespace internal
 
     using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
 
+    static_assert(
+      n_values == 1 || stride == 1,
+      "Either n_values or stride has to be one for correct data access!");
     // If n_values > 1, we want to interpolate from a second array,
     // placed in the same array immediately after the main data. This
     // is used to interpolate normal derivatives onto faces.
-    const Number *values_2 =
-      n_values > 1 ? values + Utilities::fixed_power<dim>(n_shapes) : nullptr;
-
-    AssertDimension(n_shapes, 2);
-    for (unsigned int i = 0; i < renumber.size(); ++i)
-      AssertDimension(renumber[i], i);
 
     std::array<Number3, dim + n_values> result;
     if (dim == 0)
@@ -2303,23 +2302,23 @@ namespace internal
         // we only need the value on faces of a 1d element
         result[0] = values[0];
         if (n_values > 1)
-          result[1] = values_2[0];
+          result[1] = values[1];
       }
     else if (dim == 1)
       {
         // gradient
-        result[0] = Number3(values[1] - values[0]);
+        result[0] = Number3(values[stride] - values[0]);
         // values
         result[1] = Number3(values[0]) + p[0] * result[0];
         if (n_values > 1)
-          result[2] = Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]);
+          result[2] = Number3(values[2]) + p[0] * (values[3] - values[2]);
       }
     else if (dim == 2)
       {
-        const Number3 val10 = Number3(values[1] - values[0]);
-        const Number3 val32 = Number3(values[3] - values[2]);
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
         const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2]) + p[0] * val32;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
 
         // gradient
         result[0] = val10 + p[1] * (val32 - val10);
@@ -2331,25 +2330,25 @@ namespace internal
         if (n_values > 1)
           {
             const Number3 tmp0_2 =
-              Number3(values_2[0]) + p[0] * (values_2[1] - values_2[0]);
+              Number3(values[4]) + p[0] * (values[5] - values[4]);
             const Number3 tmp1_2 =
-              Number3(values_2[2]) + p[0] * (values_2[3] - values_2[0]);
+              Number3(values[6]) + p[0] * (values[7] - values[6]);
             result[3] = tmp0_2 + p[1] * (tmp1_2 - tmp0_2);
           }
       }
     else if (dim == 3)
       {
-        const Number3 val10 = Number3(values[1] - values[0]);
-        const Number3 val32 = Number3(values[3] - values[2]);
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
         const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2]) + p[0] * val32;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
         const Number3 tmp10 = tmp1 - tmp0;
         const Number3 tmpy0 = tmp0 + p[1] * tmp10;
 
-        const Number3 val54 = Number3(values[5] - values[4]);
-        const Number3 val76 = Number3(values[7] - values[6]);
-        const Number3 tmp2  = Number3(values[4]) + p[0] * val54;
-        const Number3 tmp3  = Number3(values[6]) + p[0] * val76;
+        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
+        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
+        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
+        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
         const Number3 tmp32 = tmp3 - tmp2;
         const Number3 tmpy1 = tmp2 + p[1] * tmp32;
 
@@ -2419,8 +2418,8 @@ namespace internal
     std::array<Number3, dim + 1> result;
     if (d_linear)
       {
-        result = evaluate_tensor_product_value_and_gradient_linear(
-          poly.size(), values.data(), p, renumber);
+        result =
+          evaluate_tensor_product_value_and_gradient_linear(values.data(), p);
       }
     else
       {
@@ -2548,23 +2547,15 @@ namespace internal
 
 
 
-  template <int dim, typename Number, typename Number2>
+  template <int dim, typename Number, typename Number2, int stride = 1>
   inline typename ProductTypeNoPoint<Number, Number2>::type
-  evaluate_tensor_product_value_linear(
-    const unsigned int               n_shapes,
-    const Number                    *values,
-    const Point<dim, Number2>       &p,
-    const std::vector<unsigned int> &renumber = {})
+  evaluate_tensor_product_value_linear(const Number              *values,
+                                       const Point<dim, Number2> &p)
   {
-    (void)n_shapes;
     static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
 
     using Number3 = typename ProductTypeNoPoint<Number, Number2>::type;
 
-    AssertDimension(n_shapes, 2);
-    for (unsigned int i = 0; i < renumber.size(); ++i)
-      AssertDimension(renumber[i], i);
-
     if (dim == 0)
       {
         // we only need the value on faces of a 1d element
@@ -2572,28 +2563,28 @@ namespace internal
       }
     else if (dim == 1)
       {
-        return Number3(values[0]) + p[0] * Number3(values[1] - values[0]);
+        return Number3(values[0]) + p[0] * Number3(values[stride] - values[0]);
       }
     else if (dim == 2)
       {
-        const Number3 val10 = Number3(values[1] - values[0]);
-        const Number3 val32 = Number3(values[3] - values[2]);
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
         const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2]) + p[0] * val32;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
         return tmp0 + p[1] * (tmp1 - tmp0);
       }
     else if (dim == 3)
       {
-        const Number3 val10 = Number3(values[1] - values[0]);
-        const Number3 val32 = Number3(values[3] - values[2]);
+        const Number3 val10 = Number3(values[stride] - values[0]);
+        const Number3 val32 = Number3(values[3 * stride] - values[2 * stride]);
         const Number3 tmp0  = Number3(values[0]) + p[0] * val10;
-        const Number3 tmp1  = Number3(values[2]) + p[0] * val32;
+        const Number3 tmp1  = Number3(values[2 * stride]) + p[0] * val32;
         const Number3 tmpy0 = tmp0 + p[1] * (tmp1 - tmp0);
 
-        const Number3 val54 = Number3(values[5] - values[4]);
-        const Number3 val76 = Number3(values[7] - values[6]);
-        const Number3 tmp2  = Number3(values[4]) + p[0] * val54;
-        const Number3 tmp3  = Number3(values[6]) + p[0] * val76;
+        const Number3 val54 = Number3(values[5 * stride] - values[4 * stride]);
+        const Number3 val76 = Number3(values[7 * stride] - values[6 * stride]);
+        const Number3 tmp2  = Number3(values[4 * stride]) + p[0] * val54;
+        const Number3 tmp3  = Number3(values[6 * stride]) + p[0] * val76;
         const Number3 tmpy1 = tmp2 + p[1] * (tmp3 - tmp2);
 
         return tmpy0 + p[2] * (tmpy1 - tmpy0);
@@ -2617,10 +2608,7 @@ namespace internal
     typename ProductTypeNoPoint<Number, Number2>::type result;
     if (d_linear)
       {
-        result = evaluate_tensor_product_value_linear(poly.size(),
-                                                      values.data(),
-                                                      p,
-                                                      renumber);
+        result = evaluate_tensor_product_value_linear(values.data(), p);
       }
     else
       {
@@ -3049,19 +3037,15 @@ namespace internal
             int  n_values = 1>
   inline void
   integrate_add_tensor_product_value_and_gradient_linear(
-    const unsigned int             n_shapes,
     const Number2                 *value,
     const Tensor<1, dim, Number2> &gradient,
     Number2                       *values,
     const Point<dim, Number>      &p)
   {
-    (void)n_shapes;
     static_assert(0 <= dim && dim <= 3, "Only dim=0,1,2,3 implemented");
     static_assert(1 <= n_values && n_values <= 2,
                   "Only n_values=1,2 implemented");
 
-    AssertDimension(n_shapes, 2);
-
     // Note that 'add' is a template argument, so the compiler will remove
     // these checks
     if (dim == 0)
@@ -3210,7 +3194,11 @@ namespace internal
    * function depending on if values should be added to or set and if
    * polynomials are linear.
    */
-  template <int dim, typename Number, typename Number2, int n_values = 1>
+  template <bool is_linear,
+            int  dim,
+            typename Number,
+            typename Number2,
+            int n_values = 1>
   inline void
   integrate_tensor_product_value_and_gradient(
     const dealii::ndarray<Number, 2, dim> *shapes,
@@ -3219,7 +3207,6 @@ namespace internal
     const Tensor<1, dim, Number2>         &gradient,
     Number2                               *values,
     const Point<dim, Number>              &p,
-    const bool                             is_linear,
     const bool                             do_add)
   {
     if (do_add)
@@ -3230,7 +3217,7 @@ namespace internal
             Number,
             Number2,
             true,
-            n_values>(n_shapes, value, gradient, values, p);
+            n_values>(value, gradient, values, p);
         else
           internal::integrate_add_tensor_product_value_and_gradient_shapes<
             dim,
@@ -3247,7 +3234,7 @@ namespace internal
             Number,
             Number2,
             false,
-            n_values>(n_shapes, value, gradient, values, p);
+            n_values>(value, gradient, values, p);
         else
           internal::integrate_add_tensor_product_value_and_gradient_shapes<
             dim,
@@ -3382,16 +3369,12 @@ namespace internal
    */
   template <int dim, typename Number, typename Number2, bool add>
   inline void
-  integrate_add_tensor_product_value_linear(const unsigned int        n_shapes,
-                                            const Number2            &value,
+  integrate_add_tensor_product_value_linear(const Number2            &value,
                                             Number2                  *values,
                                             const Point<dim, Number> &p)
   {
-    (void)n_shapes;
     static_assert(dim >= 0 && dim <= 3, "Only dim=0,1,2,3 implemented");
 
-    AssertDimension(n_shapes, 2);
-
     if (dim == 0)
       {
         if (add)
@@ -3481,14 +3464,13 @@ namespace internal
    * function depending on if values should be added to or set and if
    * polynomials are linear.
    */
-  template <int dim, typename Number, typename Number2>
+  template <bool is_linear, int dim, typename Number, typename Number2>
   inline void
   integrate_tensor_product_value(const dealii::ndarray<Number, 2, dim> *shapes,
                                  const unsigned int        n_shapes,
                                  const Number2            &value,
                                  Number2                  *values,
                                  const Point<dim, Number> &p,
-                                 const bool                is_linear,
                                  const bool                do_add)
   {
     if (do_add)
@@ -3497,8 +3479,7 @@ namespace internal
           internal::integrate_add_tensor_product_value_linear<dim,
                                                               Number,
                                                               Number2,
-                                                              true>(n_shapes,
-                                                                    value,
+                                                              true>(value,
                                                                     values,
                                                                     p);
         else
@@ -3516,8 +3497,7 @@ namespace internal
           internal::integrate_add_tensor_product_value_linear<dim,
                                                               Number,
                                                               Number2,
-                                                              false>(n_shapes,
-                                                                     value,
+                                                              false>(value,
                                                                      values,
                                                                      p);
         else
index eb094cc43284cfebe2901f722c0ed4538b1ac188..9e674a4f715f248113d4407672469ebe3a3fe0d4 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/config.h>
 
+#include "deal.II/base/floating_point_comparator.h"
 #include <deal.II/base/aligned_vector.h>
 #include <deal.II/base/exceptions.h>
 #include <deal.II/base/vectorization.h>
@@ -31,6 +32,8 @@
 #include <deal.II/fe/mapping_q.h>
 #include <deal.II/fe/mapping_related_data.h>
 
+#include <deal.II/matrix_free/mapping_info_storage.h>
+
 #include <memory>
 
 
@@ -179,6 +182,49 @@ namespace NonMatching
           }
       }
     };
+
+    template <int dim, int spacedim = dim>
+    dealii::internal::MatrixFreeFunctions::GeometryType
+    compute_geometry_type(
+      const double diameter,
+      const std::vector<DerivativeForm<1, dim, spacedim, double>>
+        &inverse_jacobians)
+    {
+      const auto   jac_0 = inverse_jacobians[0];
+      const double zero_tolerance_double =
+        1. / diameter * std::numeric_limits<double>::epsilon() * 1024.;
+      bool jacobian_constant = true;
+      for (unsigned int q = 1; q < inverse_jacobians.size(); ++q)
+        {
+          const DerivativeForm<1, dim, spacedim> &jac = inverse_jacobians[q];
+          for (unsigned int d = 0; d < dim; ++d)
+            for (unsigned int e = 0; e < spacedim; ++e)
+              if (std::fabs(jac_0[d][e] - jac[d][e]) > zero_tolerance_double)
+                jacobian_constant = false;
+          if (!jacobian_constant)
+            break;
+        }
+
+      // check whether the Jacobian is diagonal to machine
+      // accuracy
+      bool cell_cartesian = jacobian_constant;
+      for (unsigned int d = 0; d < dim; ++d)
+        for (unsigned int e = 0; e < dim; ++e)
+          if (d != e)
+            if (std::fabs(jac_0[d][e]) > zero_tolerance_double)
+              {
+                cell_cartesian = false;
+                break;
+              }
+
+      // return cell type
+      if (cell_cartesian)
+        return dealii::internal::MatrixFreeFunctions::GeometryType::cartesian;
+      else if (jacobian_constant)
+        return dealii::internal::MatrixFreeFunctions::GeometryType::affine;
+      else
+        return dealii::internal::MatrixFreeFunctions::GeometryType::general;
+    }
   } // namespace internal
 
   /**
@@ -412,26 +458,44 @@ namespace NonMatching
     boost::signals2::connection
     connect_is_reinitialized(const std::function<void()> &set_is_reinitialized);
 
+    /**
+     * Compute the geometry index offset of the current cell/face.
+     */
+    template <bool is_face>
+    unsigned int
+    compute_geometry_index_offset(const unsigned int cell_index,
+                                  const unsigned int face_number) const;
+
     /**
      * Compute the unit points index offset for the current cell/face.
      */
     unsigned int
-    compute_unit_point_index_offset(const unsigned int cell_index,
-                                    const unsigned int face_number) const;
+    compute_unit_point_index_offset(const unsigned int geometry_index) const;
 
     /**
      * Compute the data index offset for the current cell/face.
      */
     unsigned int
-    compute_data_index_offset(const unsigned int cell_index,
-                              const unsigned int face_number) const;
+    compute_data_index_offset(const unsigned int geometry_index) const;
+
+    /**
+     * Compute the data index offset for the current cell/face.
+     */
+    unsigned int
+    compute_compressed_data_index_offset(
+      const unsigned int geometry_index) const;
 
     /**
      * Get number of unvectorized quadrature points.
      */
     unsigned int
-    get_n_q_points_unvectorized(const unsigned int cell_index,
-                                const unsigned int face_number) const;
+    get_n_q_points_unvectorized(const unsigned int geometry_index) const;
+
+    /**
+     * Get cell geometry type.
+     */
+    dealii::internal::MatrixFreeFunctions::GeometryType
+    get_cell_type(const unsigned int geometry_index) const;
 
     /**
      * Return cell iterator.
@@ -503,7 +567,9 @@ namespace NonMatching
                        const unsigned int         n_q_points,
                        const unsigned int         n_q_points_unvectorized,
                        const MappingData         &mapping_data,
-                       const std::vector<double> &weights);
+                       const std::vector<double> &weights,
+                       const unsigned int compressed_unit_point_index_offset,
+                       const bool         affine_cell);
 
     /**
      * Compute the compressed cell index.
@@ -512,11 +578,18 @@ namespace NonMatching
     compute_compressed_cell_index(const unsigned int cell_index) const;
 
     /**
-     * Compute the geometry index offset of the current cell/face.
+     * Compute the mapping information for cells/surface.
      */
-    unsigned int
-    compute_geometry_index_offset(const unsigned int cell_index,
-                                  const unsigned int face_number) const;
+    template <typename ContainerType, typename QuadratureType>
+    void
+    do_reinit_cells(
+      const ContainerType               &cell_iterator_range,
+      const std::vector<QuadratureType> &quadrature_vector,
+      const unsigned int                 n_unfiltered_cells,
+      const std::function<
+        void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+             const QuadratureType &quadrature,
+             MappingData          &mapping_data)> &compute_mapping_data);
 
     /**
      * Enum class for reinitialized states.
@@ -552,7 +625,7 @@ namespace NonMatching
     /**
      * Offset to point to the first unit point of a cell/face.
      */
-    AlignedVector<unsigned int> unit_points_index;
+    std::vector<unsigned int> unit_points_index;
 
     /**
      * A pointer to the internal data of the underlying mapping.
@@ -581,10 +654,22 @@ namespace NonMatching
     const AdditionalData additional_data;
 
     /**
-     * Stores the index offset into the arrays @p JxW_values, @p jacobians,
-     * @p inverse_jacobians and @p normal_vectors.
+     * Stores whether a cell is Cartesian (cell type 0), has constant
+     * transform data (Jacobians) (cell type 1), or is general (cell type
+     * 3). Type 2 is only used for faces and no cells are assigned this
+     * value.
+     */
+    std::vector<dealii::internal::MatrixFreeFunctions::GeometryType> cell_type;
+
+    /**
+     * Stores the index offset into the arrays @p JxW_values and @p normal_vectors.
+     */
+    std::vector<unsigned int> data_index_offsets;
+
+    /**
+     * Stores the index offset into the arrays @p jacobians and @p inverse_jacobians.
      */
-    AlignedVector<unsigned int> data_index_offsets;
+    std::vector<unsigned int> compressed_data_index_offsets;
 
     /**
      * The storage of the Jacobian determinant times the quadrature weight on
@@ -605,7 +690,7 @@ namespace NonMatching
      * The storage of contravariant transformation on quadrature points, i.e.,
      * the Jacobians of the transformation from the unit to the real cell.
      *
-     * Indexed by @p data_index_offsets.
+     * Indexed by @p compressed_data_index_offsets.
      */
     AlignedVector<DerivativeForm<1, dim, spacedim, Number>> jacobians;
 
@@ -614,7 +699,7 @@ namespace NonMatching
      * the inverse Jacobians of the transformation from the
      * unit to the real cell.
      *
-     * Indexed by @p data_index_offsets.
+     * Indexed by @p compressed_data_index_offsets.
      */
     AlignedVector<DerivativeForm<1, spacedim, dim, Number>> inverse_jacobians;
 
@@ -718,6 +803,8 @@ namespace NonMatching
     n_q_points_unvectorized.clear();
     unit_points_index.clear();
     data_index_offsets.clear();
+    compressed_data_index_offsets.clear();
+    cell_type.clear();
   }
 
 
@@ -783,12 +870,32 @@ namespace NonMatching
                                           internal_mapping_data,
                                           mapping_data);
 
+    // check for cartesian/affine cell
+    if (!quadrature.empty() &&
+        update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+      {
+        cell_type.push_back(
+          internal::compute_geometry_type(cell->diameter(),
+                                          mapping_data.inverse_jacobians));
+      }
+    else
+      cell_type.push_back(
+        dealii::internal::MatrixFreeFunctions::GeometryType::general);
+
     // store mapping data
-    store_mapping_data(0,
-                       n_q_points_data,
-                       n_q_points_unvectorized[0],
-                       mapping_data,
-                       quadrature.get_weights());
+    store_mapping_data(
+      0,
+      n_q_points_data,
+      n_q_points_unvectorized[0],
+      mapping_data,
+      quadrature.get_weights(),
+      0,
+      cell_type.back() <=
+        dealii::internal::MatrixFreeFunctions::GeometryType::affine);
+
+    unit_points_index.push_back(0);
+    data_index_offsets.push_back(0);
+    compressed_data_index_offsets.push_back(0);
 
     state = State::single_cell;
     is_reinitialized();
@@ -820,12 +927,16 @@ namespace NonMatching
 
 
   template <int dim, int spacedim, typename Number>
-  template <typename ContainerType>
+  template <typename ContainerType, typename QuadratureType>
   void
-  MappingInfo<dim, spacedim, Number>::reinit_cells(
-    const ContainerType                &cell_iterator_range,
-    const std::vector<Quadrature<dim>> &quadrature_vector,
-    const unsigned int                  n_unfiltered_cells)
+  MappingInfo<dim, spacedim, Number>::do_reinit_cells(
+    const ContainerType               &cell_iterator_range,
+    const std::vector<QuadratureType> &quadrature_vector,
+    const unsigned int                 n_unfiltered_cells,
+    const std::function<
+      void(const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+           const QuadratureType &quadrature,
+           MappingData          &mapping_data)> &compute_mapping_data)
   {
     clear();
 
@@ -839,6 +950,8 @@ namespace NonMatching
 
     n_q_points_unvectorized.reserve(n_cells);
 
+    cell_type.reserve(n_cells);
+
     if (additional_data.store_cells)
       cell_level_and_indices.resize(n_cells);
 
@@ -873,9 +986,8 @@ namespace NonMatching
       cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
                                                  numbers::invalid_unsigned_int);
 
-    MappingData                mapping_data;
-    CellSimilarity::Similarity cell_similarity =
-      CellSimilarity::Similarity::none;
+    MappingData  mapping_data;
+    MappingData  mapping_data_last_cell;
     unsigned int cell_index = 0;
     for (const auto &cell : cell_iterator_range)
       {
@@ -885,32 +997,97 @@ namespace NonMatching
             cell_level_and_indices[cell_index] = {cell->level(), cell->index()};
           }
 
+        const auto &quadrature = quadrature_vector[cell_index];
+        const bool  empty      = quadrature.empty();
+
         // store unit points
         const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
           n_q_points_unvectorized[cell_index]);
         store_unit_points(unit_points_index[cell_index],
                           n_q_points,
                           n_q_points_unvectorized[cell_index],
-                          quadrature_vector[cell_index].get_points());
+                          quadrature.get_points());
 
         // compute mapping data
-        internal::ComputeMappingDataHelper<dim, spacedim>::
-          compute_mapping_data_for_quadrature(mapping,
-                                              update_flags_mapping,
-                                              cell,
-                                              cell_similarity,
-                                              quadrature_vector[cell_index],
-                                              internal_mapping_data,
-                                              mapping_data);
+        compute_mapping_data(cell, quadrature, mapping_data);
 
         // store mapping data
         const unsigned int n_q_points_data =
           compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
+
+        // check for cartesian/affine cell
+        if (!empty &&
+            update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+          {
+            cell_type.push_back(
+              internal::compute_geometry_type(cell->diameter(),
+                                              mapping_data.inverse_jacobians));
+          }
+        else
+          cell_type.push_back(
+            dealii::internal::MatrixFreeFunctions::GeometryType::general);
+
+        if (cell_index > 0)
+          {
+            // check if current and last cell are affine
+            const bool affine_cells =
+              cell_type[cell_index] <=
+                dealii::internal::MatrixFreeFunctions::affine &&
+              cell_type[cell_index - 1] <=
+                dealii::internal::MatrixFreeFunctions::affine;
+
+            // create a comparator to compare inverse Jacobian of current
+            // and last cell
+            FloatingPointComparator<double> comparator(
+              1e4 / cell->diameter() * std::numeric_limits<double>::epsilon() *
+              1024.);
+
+            // we can only compare if current and last cell have at least
+            // one quadrature point and both cells are at least affine
+            const auto comparison_result =
+              (!affine_cells || mapping_data.inverse_jacobians.empty() ||
+               mapping_data_last_cell.inverse_jacobians.empty()) ?
+                FloatingPointComparator<double>::ComparisonResult::less :
+                comparator.compare(mapping_data.inverse_jacobians[0],
+                                   mapping_data_last_cell.inverse_jacobians[0]);
+
+            // we can compress the Jacobians and inverse Jacobians if
+            // inverse Jacobians are equal and cells are affine
+            if (affine_cells &&
+                comparison_result ==
+                  FloatingPointComparator<double>::ComparisonResult::equal)
+              {
+                compressed_data_index_offsets.push_back(
+                  compressed_data_index_offsets.back());
+              }
+            else
+              {
+                const unsigned int n_compressed_data_last_cell =
+                  cell_type[cell_index - 1] <=
+                      dealii::internal::MatrixFreeFunctions::affine ?
+                    1 :
+                    compute_n_q_points<Number>(
+                      n_q_points_unvectorized[cell_index - 1]);
+
+                compressed_data_index_offsets.push_back(
+                  compressed_data_index_offsets.back() +
+                  n_compressed_data_last_cell);
+              }
+          }
+        else
+          compressed_data_index_offsets.push_back(0);
+
+        // cache mapping_data from last cell
+        mapping_data_last_cell = mapping_data;
+
         store_mapping_data(data_index_offsets[cell_index],
                            n_q_points_data,
                            n_q_points_unvectorized[cell_index],
                            mapping_data,
-                           quadrature_vector[cell_index].get_weights());
+                           quadrature.get_weights(),
+                           compressed_data_index_offsets[cell_index],
+                           cell_type[cell_index] <=
+                             dealii::internal::MatrixFreeFunctions::affine);
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -919,12 +1096,47 @@ namespace NonMatching
         ++cell_index;
       }
 
+    // TODO: release allocated memory from compressed data vectors
+
     state = State::cell_vector;
     is_reinitialized();
   }
 
 
 
+  template <int dim, int spacedim, typename Number>
+  template <typename ContainerType>
+  void
+  MappingInfo<dim, spacedim, Number>::reinit_cells(
+    const ContainerType                &cell_iterator_range,
+    const std::vector<Quadrature<dim>> &quadrature_vector,
+    const unsigned int                  n_unfiltered_cells)
+  {
+    auto compute_mapping_data_for_cells =
+      [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+          const Quadrature<dim> &quadrature,
+          MappingData           &mapping_data) {
+        CellSimilarity::Similarity cell_similarity =
+          CellSimilarity::Similarity::none;
+        internal::ComputeMappingDataHelper<dim, spacedim>::
+          compute_mapping_data_for_quadrature(mapping,
+                                              update_flags_mapping,
+                                              cell,
+                                              cell_similarity,
+                                              quadrature,
+                                              internal_mapping_data,
+                                              mapping_data);
+      };
+
+    do_reinit_cells<ContainerType, Quadrature<dim>>(
+      cell_iterator_range,
+      quadrature_vector,
+      n_unfiltered_cells,
+      compute_mapping_data_for_cells);
+  }
+
+
+
   template <int dim, int spacedim, typename Number>
   template <typename Iterator>
   void
@@ -933,8 +1145,6 @@ namespace NonMatching
     const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
     const unsigned int                                 n_unfiltered_cells)
   {
-    clear();
-
     Assert(
       additional_data.use_global_weights == false,
       ExcMessage(
@@ -942,65 +1152,13 @@ namespace NonMatching
 
     Assert(additional_data.store_cells == false, ExcNotImplemented());
 
-    do_cell_index_compression =
-      n_unfiltered_cells != numbers::invalid_unsigned_int;
-
     if (update_flags_mapping & (update_JxW_values | update_normal_vectors))
       update_flags_mapping |= update_covariant_transformation;
 
-    const unsigned int n_cells = quadrature_vector.size();
-    AssertDimension(n_cells,
-                    std::distance(cell_iterator_range.begin(),
-                                  cell_iterator_range.end()));
-
-    n_q_points_unvectorized.reserve(n_cells);
-
-    // fill unit points index offset vector
-    unit_points_index.reserve(n_cells + 1);
-    unit_points_index.push_back(0);
-    data_index_offsets.reserve(n_cells + 1);
-    data_index_offsets.push_back(0);
-    for (const auto &quadrature : quadrature_vector)
-      {
-        const unsigned int n_points = quadrature.size();
-        n_q_points_unvectorized.push_back(n_points);
-
-        const unsigned int n_q_points =
-          compute_n_q_points<VectorizedArrayType>(n_points);
-        unit_points_index.push_back(unit_points_index.back() + n_q_points);
-
-        const unsigned int n_q_points_data =
-          compute_n_q_points<Number>(n_points);
-        data_index_offsets.push_back(data_index_offsets.back() +
-                                     n_q_points_data);
-      }
-
-    const unsigned int n_unit_points = unit_points_index.back();
-    const unsigned int n_data_points = data_index_offsets.back();
-
-    // resize data vectors
-    resize_unit_points(n_unit_points);
-    resize_data_fields(n_data_points);
-
-    if (do_cell_index_compression)
-      cell_index_to_compressed_cell_index.resize(n_unfiltered_cells,
-                                                 numbers::invalid_unsigned_int);
-
-    MappingData  mapping_data;
-    unsigned int cell_index = 0;
-    for (const auto &cell : cell_iterator_range)
-      {
-        const auto &quadrature = quadrature_vector[cell_index];
-
-        // store unit points
-        const unsigned int n_q_points = compute_n_q_points<VectorizedArrayType>(
-          n_q_points_unvectorized[cell_index]);
-        store_unit_points(unit_points_index[cell_index],
-                          n_q_points,
-                          n_q_points_unvectorized[cell_index],
-                          quadrature_vector[cell_index].get_points());
-
-        // compute mapping data
+    auto compute_mapping_data_for_surface =
+      [&](const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+          const ImmersedSurfaceQuadrature<dim> &quadrature,
+          MappingData                          &mapping_data) {
         internal::ComputeMappingDataHelper<dim, spacedim>::
           compute_mapping_data_for_immersed_surface_quadrature(
             mapping,
@@ -1009,25 +1167,13 @@ namespace NonMatching
             quadrature,
             internal_mapping_data,
             mapping_data);
+      };
 
-        // store mapping data
-        const unsigned int n_q_points_data =
-          compute_n_q_points<Number>(n_q_points_unvectorized[cell_index]);
-        store_mapping_data(data_index_offsets[cell_index],
-                           n_q_points_data,
-                           n_q_points_unvectorized[cell_index],
-                           mapping_data,
-                           quadrature_vector[cell_index].get_weights());
-
-        if (do_cell_index_compression)
-          cell_index_to_compressed_cell_index[cell->active_cell_index()] =
-            cell_index;
-
-        ++cell_index;
-      }
-
-    state = State::cell_vector;
-    is_reinitialized();
+    do_reinit_cells<IteratorRange<Iterator>, ImmersedSurfaceQuadrature<dim>>(
+      cell_iterator_range,
+      quadrature_vector,
+      n_unfiltered_cells,
+      compute_mapping_data_for_surface);
   }
 
 
@@ -1065,6 +1211,8 @@ namespace NonMatching
 
     n_q_points_unvectorized.reserve(n_faces);
 
+    cell_type.reserve(n_faces);
+
     // fill unit points index offset vector
     unit_points_index.resize(n_faces + 1);
     data_index_offsets.resize(n_faces + 1);
@@ -1156,13 +1304,21 @@ namespace NonMatching
                                                        internal_mapping_data,
                                                        mapping_data);
 
+            cell_type.push_back(
+              dealii::internal::MatrixFreeFunctions::GeometryType::general);
+
+            compressed_data_index_offsets.push_back(
+              data_index_offsets[current_face_index]);
+
             const unsigned int n_q_points_data = compute_n_q_points<Number>(
               n_q_points_unvectorized[current_face_index]);
             store_mapping_data(data_index_offsets[current_face_index],
                                n_q_points_data,
                                n_q_points_unvectorized[current_face_index],
                                mapping_data,
-                               quadrature_on_face.get_weights());
+                               quadrature_on_face.get_weights(),
+                               data_index_offsets[current_face_index],
+                               false);
           }
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -1189,22 +1345,18 @@ namespace NonMatching
   template <int dim, int spacedim, typename Number>
   unsigned int
   MappingInfo<dim, spacedim, Number>::get_n_q_points_unvectorized(
-    const unsigned int cell_index,
-    const unsigned int face_number) const
+    const unsigned int geometry_index) const
   {
-    if (cell_index == numbers::invalid_unsigned_int &&
-        face_number == numbers::invalid_unsigned_int)
-      {
-        Assert(state == State::single_cell,
-               ExcMessage(
-                 "This mapping info is not reinitialized for a single cell!"));
-        return n_q_points_unvectorized[0];
-      }
-    else
-      {
-        return n_q_points_unvectorized[compute_geometry_index_offset(
-          cell_index, face_number)];
-      }
+    return n_q_points_unvectorized[geometry_index];
+  }
+
+
+  template <int dim, int spacedim, typename Number>
+  dealii::internal::MatrixFreeFunctions::GeometryType
+  MappingInfo<dim, spacedim, Number>::get_cell_type(
+    const unsigned int geometry_index) const
+  {
+    return cell_type[geometry_index];
   }
 
 
@@ -1244,14 +1396,18 @@ namespace NonMatching
 
 
   template <int dim, int spacedim, typename Number>
+  template <bool is_face>
   unsigned int
   MappingInfo<dim, spacedim, Number>::compute_geometry_index_offset(
     const unsigned int cell_index,
     const unsigned int face_number) const
   {
+    if (cell_index == numbers::invalid_unsigned_int)
+      return 0;
+
     const unsigned int compressed_cell_index =
       compute_compressed_cell_index(cell_index);
-    if (face_number == numbers::invalid_unsigned_int)
+    if (!is_face)
       {
         Assert(state == State::cell_vector,
                ExcMessage(
@@ -1347,7 +1503,9 @@ namespace NonMatching
     const unsigned int              n_q_points,
     const unsigned int              n_q_points_unvectorized,
     const MappingInfo::MappingData &mapping_data,
-    const std::vector<double>      &weights)
+    const std::vector<double>      &weights,
+    const unsigned int              compressed_unit_point_index_offset,
+    const bool                      affine_cell)
   {
     const unsigned int n_lanes =
       dealii::internal::VectorizedArrayTrait<Number>::width();
@@ -1355,22 +1513,28 @@ namespace NonMatching
     for (unsigned int q = 0; q < n_q_points; ++q)
       {
         const unsigned int offset = unit_points_index_offset + q;
+        const unsigned int compressed_offset =
+          compressed_unit_point_index_offset + q;
         for (unsigned int v = 0;
              v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
              ++v)
           {
-            if (update_flags_mapping & UpdateFlags::update_jacobians)
-              for (unsigned int d = 0; d < dim; ++d)
-                for (unsigned int s = 0; s < spacedim; ++s)
-                  dealii::internal::VectorizedArrayTrait<Number>::get(
-                    jacobians[offset][d][s], v) =
-                    mapping_data.jacobians[q * n_lanes + v][d][s];
-            if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
-              for (unsigned int d = 0; d < dim; ++d)
-                for (unsigned int s = 0; s < spacedim; ++s)
-                  dealii::internal::VectorizedArrayTrait<Number>::get(
-                    inverse_jacobians[offset][s][d], v) =
-                    mapping_data.inverse_jacobians[q * n_lanes + v][s][d];
+            if (q == 0 || !affine_cell)
+              {
+                if (update_flags_mapping & UpdateFlags::update_jacobians)
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int s = 0; s < spacedim; ++s)
+                      dealii::internal::VectorizedArrayTrait<Number>::get(
+                        jacobians[compressed_offset][d][s], v) =
+                        mapping_data.jacobians[q * n_lanes + v][d][s];
+                if (update_flags_mapping &
+                    UpdateFlags::update_inverse_jacobians)
+                  for (unsigned int d = 0; d < dim; ++d)
+                    for (unsigned int s = 0; s < spacedim; ++s)
+                      dealii::internal::VectorizedArrayTrait<Number>::get(
+                        inverse_jacobians[compressed_offset][s][d], v) =
+                        mapping_data.inverse_jacobians[q * n_lanes + v][s][d];
+              }
             if (update_flags_mapping & UpdateFlags::update_JxW_values)
               {
                 if (additional_data.use_global_weights)
@@ -1447,7 +1611,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_unit_point(
     const unsigned int offset) const
   {
-    return &unit_points[offset];
+    return unit_points.data() + offset;
   }
 
 
@@ -1459,7 +1623,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_unit_point_faces(
     const unsigned int offset) const
   {
-    return &unit_points_faces[offset];
+    return unit_points_faces.data() + offset;
   }
 
 
@@ -1469,7 +1633,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_real_point(
     const unsigned int offset) const
   {
-    return &real_points[offset];
+    return real_points.data() + offset;
   }
 
 
@@ -1477,23 +1641,9 @@ namespace NonMatching
   template <int dim, int spacedim, typename Number>
   unsigned int
   MappingInfo<dim, spacedim, Number>::compute_unit_point_index_offset(
-    const unsigned int cell_index,
-    const unsigned int face_number) const
+    const unsigned int geometry_index) const
   {
-    if (cell_index == numbers::invalid_unsigned_int &&
-        face_number == numbers::invalid_unsigned_int)
-      {
-        Assert(state == State::single_cell,
-               ExcMessage(
-                 "This mapping info is not reinitialized for a single cell!"));
-        return 0;
-      }
-    else
-      {
-        const unsigned int offset =
-          compute_geometry_index_offset(cell_index, face_number);
-        return unit_points_index[offset];
-      }
+    return unit_points_index[geometry_index];
   }
 
 
@@ -1501,23 +1651,18 @@ namespace NonMatching
   template <int dim, int spacedim, typename Number>
   unsigned int
   MappingInfo<dim, spacedim, Number>::compute_data_index_offset(
-    const unsigned int cell_index,
-    const unsigned int face_number) const
+    const unsigned int geometry_index) const
   {
-    if (cell_index == numbers::invalid_unsigned_int &&
-        face_number == numbers::invalid_unsigned_int)
-      {
-        Assert(state == State::single_cell,
-               ExcMessage(
-                 "This mapping info is not reinitialized for a single cell!"));
-        return 0;
-      }
-    else
-      {
-        const unsigned int offset =
-          compute_geometry_index_offset(cell_index, face_number);
-        return data_index_offsets[offset];
-      }
+    return data_index_offsets[geometry_index];
+  }
+
+
+  template <int dim, int spacedim, typename Number>
+  unsigned int
+  MappingInfo<dim, spacedim, Number>::compute_compressed_data_index_offset(
+    const unsigned int geometry_index) const
+  {
+    return compressed_data_index_offsets[geometry_index];
   }
 
 
@@ -1527,7 +1672,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_jacobian(
     const unsigned int offset) const
   {
-    return &jacobians[offset];
+    return jacobians.data() + offset;
   }
 
 
@@ -1537,7 +1682,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_inverse_jacobian(
     const unsigned int offset) const
   {
-    return &inverse_jacobians[offset];
+    return inverse_jacobians.data() + offset;
   }
 
 
@@ -1546,7 +1691,7 @@ namespace NonMatching
   MappingInfo<dim, spacedim, Number>::get_normal_vector(
     const unsigned int offset) const
   {
-    return &normal_vectors[offset];
+    return normal_vectors.data() + offset;
   }
 
 
@@ -1555,7 +1700,7 @@ namespace NonMatching
   inline const Number *
   MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
   {
-    return &JxW_values[offset];
+    return JxW_values.data() + offset;
   }
 
 

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.