]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Vectorize MappingInfo 15137/head
authorMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Fri, 31 Mar 2023 14:58:24 +0000 (16:58 +0200)
committerMaximilian Bergbauer <maximilian.bergbauer@tum.de>
Wed, 3 May 2023 07:55:53 +0000 (09:55 +0200)
include/deal.II/base/derivative_form.h
include/deal.II/base/vectorization.h
include/deal.II/fe/mapping.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
source/fe/mapping_q.cc
tests/non_matching/mapping_info.cc
tests/non_matching/mapping_info.output

index 3cb74f5e9e203fbec55dcf9c588d0afe6fef9c60..9c36923c665f348db3d3fbe37b4b53f93f059e0c 100644 (file)
@@ -106,6 +106,13 @@ public:
   DerivativeForm &
   operator=(const Tensor<1, dim, Number> &);
 
+  /**
+   * Number conversion operator.
+   */
+  template <typename OtherNumber>
+  DerivativeForm &
+  operator=(const DerivativeForm<order, dim, spacedim, OtherNumber> &df);
+
   /**
    * Converts a DerivativeForm <order, dim, dim, Number> to Tensor<order+1, dim,
    * Number>. In particular, if order == 1 and the derivative is the Jacobian of
@@ -252,6 +259,19 @@ DerivativeForm<order, dim, spacedim, Number>::operator=(
 
 
 
+template <int order, int dim, int spacedim, typename Number>
+template <typename OtherNumber>
+inline DerivativeForm<order, dim, spacedim, Number> &
+DerivativeForm<order, dim, spacedim, Number>::operator=(
+  const DerivativeForm<order, dim, spacedim, OtherNumber> &df)
+{
+  for (unsigned int j = 0; j < spacedim; ++j)
+    (*this)[j] = df[j];
+  return *this;
+}
+
+
+
 template <int order, int dim, int spacedim, typename Number>
 inline Tensor<order, dim, Number> &
 DerivativeForm<order, dim, spacedim, Number>::operator[](const unsigned int i)
index 9aea1cb6a4367baa6808c03083d3a83497ef6285..7bac5b49cfb03a22123981f62a89bc51d44185bb 100644 (file)
@@ -5478,49 +5478,150 @@ namespace internal
   template <typename T>
   struct VectorizedArrayTrait
   {
-    using value_type                   = T;
+    /**
+     * Define scalar value type.
+     */
+    using value_type = T;
+
+    /**
+     * Define width of template type.
+     */
     static constexpr std::size_t width = 1;
 
-    static T &
-    get(T &value, unsigned int c)
+    /**
+     * Define vectorized value type for internal vectorization.
+     */
+    using vectorized_value_type = VectorizedArray<T>;
+
+    /**
+     * Define a stride which defines how often the template type T fits into the
+     * vectorized_value_type. This is useful to write vectorized templated code
+     * where the internal computation is vectorized and the user interface is
+     * optionally scalar or also vectorized.
+     */
+    static constexpr std::size_t stride = vectorized_value_type::size();
+
+    /**
+     * Get a reference to scalar value (on lane 0).
+     */
+    static value_type &
+    get(value_type &value, unsigned int c)
     {
-      AssertDimension(c, 0);
+      AssertIndexRange(c, 1);
       (void)c;
 
       return value;
     }
 
-    static const T &
-    get(const T &value, unsigned int c)
+    /**
+     * Get a read-only reference to scalar value (on lane 0).
+     */
+    static const value_type &
+    get(const value_type &value, unsigned int c)
     {
-      AssertDimension(c, 0);
+      AssertIndexRange(c, 1);
       (void)c;
 
       return value;
     }
+
+    /**
+     * Get a reference to scalar value on lane c from a vectorized values field.
+     */
+    static value_type &
+    get_from_vectorized(vectorized_value_type &values, unsigned int c)
+    {
+      AssertIndexRange(c, stride);
+
+      return values[c];
+    }
+
+    /**
+     * Get a read-only reference to scalar value on lane c from a vectorized
+     * values field.
+     */
+    static const value_type &
+    get_from_vectorized(const vectorized_value_type &values, unsigned int c)
+    {
+      AssertIndexRange(c, stride);
+
+      return values[c];
+    }
   };
 
   template <typename T, std::size_t width_>
   struct VectorizedArrayTrait<VectorizedArray<T, width_>>
   {
-    using value_type                   = T;
+    /**
+     * Define scalar value type.
+     */
+    using value_type = T;
+
+    /**
+     * Define width of template type.
+     */
     static constexpr std::size_t width = width_;
 
-    static T &
-    get(VectorizedArray<T, width_> &values, unsigned int c)
+    /**
+     * Define vectorized value type for internal vectorization.
+     */
+    using vectorized_value_type = VectorizedArray<T, width_>;
+
+    /**
+     * Define a stride which defines how often the template type
+     * VectorizedArray<T, width_> fits into the vectorized value type. This is
+     * useful to write vectorized templated code where the internal computation
+     * is vectorized and the user interface is optionally scalar or also
+     * vectorized.
+     */
+    static constexpr std::size_t stride = 1;
+
+    /**
+     * Get a reference to scalar value on lane c.
+     */
+    static value_type &
+    get(vectorized_value_type &values, unsigned int c)
     {
       AssertIndexRange(c, width_);
 
       return values[c];
     }
 
-    static const T &
-    get(const VectorizedArray<T, width_> &values, unsigned int c)
+    /**
+     * Get a read-only reference to scalar value on lane c.
+     */
+    static const value_type &
+    get(const vectorized_value_type &values, unsigned int c)
     {
       AssertIndexRange(c, width_);
 
       return values[c];
     }
+
+    /**
+     * Get a reference to vectorized values from a vectorized values field.
+     */
+    static vectorized_value_type &
+    get_from_vectorized(vectorized_value_type &values, unsigned int c)
+    {
+      (void)c;
+      AssertIndexRange(c, stride);
+
+      return values;
+    }
+
+    /**
+     * Get a read-only reference to vectorized values from a vectorized values
+     * field.
+     */
+    static const vectorized_value_type &
+    get_from_vectorized(const vectorized_value_type &values, unsigned int c)
+    {
+      (void)c;
+      AssertIndexRange(c, stride);
+
+      return values;
+    }
   };
 } // namespace internal
 
index b42684821563393b6940bde353a2503c35b7e91f..f2af0920679e827f9b86a84f55e795a4a89181ac 100644 (file)
@@ -55,8 +55,11 @@ namespace NonMatching
 {
   template <int dim>
   class FEImmersedSurfaceValues;
-  template <int dim, int spacedim>
-  class MappingInfo;
+  namespace internal
+  {
+    template <int dim, int spacedim>
+    class ComputeMappingDataHelper;
+  }
 } // namespace NonMatching
 
 
@@ -1320,7 +1323,7 @@ public:
   friend class FEFaceValues<dim, spacedim>;
   friend class FESubfaceValues<dim, spacedim>;
   friend class NonMatching::FEImmersedSurfaceValues<dim>;
-  friend class NonMatching::MappingInfo<dim, spacedim>;
+  friend class NonMatching::internal::ComputeMappingDataHelper<dim, spacedim>;
 };
 
 
index 4e07a9a447414ce8a6b1c21dce6f6c6df4b765d1..00efffd8563e49b54d8a07f8e01bf38023a3db24 100644 (file)
@@ -40,20 +40,16 @@ namespace internal
 {
   namespace FEPointEvaluation
   {
-    template <typename Number>
-    struct VectorizedPointsTypeTraits
-    {
-      using value_type                    = VectorizedArray<Number>;
-      static constexpr std::size_t stride = value_type::size();
-    };
-
-    template <typename Number, std::size_t width_>
-    struct VectorizedPointsTypeTraits<VectorizedArray<Number, width_>>
-    {
-      using value_type                    = VectorizedArray<Number, width_>;
-      static constexpr std::size_t stride = 1;
-    };
-
+    DeclException1(
+      ExcFEPointEvaluationAccessToUninitializedMappingField,
+      std::string,
+      << "You are requesting information from an FEPointEvaluation "
+      << "object for which this kind of information has not been computed. "
+      << "What information these objects compute is determined by the update_* "
+      << "flags you pass to MappingInfo() in the Constructor. Here, "
+      << "the operation you are attempting requires the <" << arg1
+      << "> flag to be set, but it was apparently not specified "
+      << "upon initialization.");
 
     /**
      * Struct to distinguish between the value and gradient types of different
@@ -65,8 +61,8 @@ namespace internal
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
       using VectorizedArrayType =
-        typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
-          Number>::value_type;
+        typename dealii::internal::VectorizedArrayTrait<
+          Number>::vectorized_value_type;
       using value_type        = Tensor<1, n_components, Number>;
       using scalar_value_type = Tensor<1, n_components, ScalarNumber>;
       using vectorized_value_type =
@@ -80,18 +76,18 @@ namespace internal
         Tensor<1, dim, Tensor<1, n_components, VectorizedArrayType>>;
 
       static void
-      read_value(const Number       vector_entry,
+      read_value(const ScalarNumber vector_entry,
                  const unsigned int component,
-                 value_type &       result)
+                 scalar_value_type &result)
       {
         AssertIndexRange(component, n_components);
         result[component] = vector_entry;
       }
 
       static void
-      write_value(Number &           vector_entry,
-                  const unsigned int component,
-                  const value_type & result)
+      write_value(VectorizedArrayType &        vector_entry,
+                  const unsigned int           component,
+                  const vectorized_value_type &result)
       {
         AssertIndexRange(component, n_components);
         vector_entry = result[component];
@@ -100,41 +96,24 @@ namespace internal
       static void
       set_gradient(const interface_vectorized_gradient_type &value,
                    const unsigned int                        vector_lane,
-                   scalar_gradient_type &                    result)
+                   gradient_type &                           result)
       {
         for (unsigned int i = 0; i < n_components; ++i)
           for (unsigned int d = 0; d < dim; ++d)
-            result[i][d] = value[d][i][vector_lane];
-      }
-
-      static void
-      set_gradient(const interface_vectorized_gradient_type &value,
-                   const unsigned int,
-                   vectorized_gradient_type &result)
-      {
-        for (unsigned int i = 0; i < n_components; ++i)
-          for (unsigned int d = 0; d < dim; ++d)
-            result[i][d] = value[d][i];
+            result[i][d] =
+              internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+                value[d][i], vector_lane);
       }
 
       static void
       get_gradient(interface_vectorized_gradient_type &value,
                    const unsigned int                  vector_lane,
-                   const scalar_gradient_type &        result)
+                   const gradient_type &               result)
       {
         for (unsigned int i = 0; i < n_components; ++i)
           for (unsigned int d = 0; d < dim; ++d)
-            value[d][i][vector_lane] = result[i][d];
-      }
-
-      static void
-      get_gradient(interface_vectorized_gradient_type &value,
-                   const unsigned int,
-                   const vectorized_gradient_type &result)
-      {
-        for (unsigned int i = 0; i < n_components; ++i)
-          for (unsigned int d = 0; d < dim; ++d)
-            value[d][i] = result[i][d];
+            internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+              value[d][i], vector_lane) = result[i][d];
       }
 
       static void
@@ -239,8 +218,8 @@ namespace internal
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
       using VectorizedArrayType =
-        typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
-          Number>::value_type;
+        typename dealii::internal::VectorizedArrayTrait<
+          Number>::vectorized_value_type;
       using value_type               = Number;
       using scalar_value_type        = ScalarNumber;
       using vectorized_value_type    = VectorizedArrayType;
@@ -249,17 +228,17 @@ namespace internal
       using vectorized_gradient_type = Tensor<1, dim, VectorizedArrayType>;
 
       static void
-      read_value(const Number vector_entry,
+      read_value(const ScalarNumber vector_entry,
                  const unsigned int,
-                 value_type &result)
+                 scalar_value_type &result)
       {
         result = vector_entry;
       }
 
       static void
-      write_value(Number &vector_entry,
+      write_value(VectorizedArrayType &vector_entry,
                   const unsigned int,
-                  const value_type &result)
+                  const vectorized_value_type &result)
       {
         vector_entry = result;
       }
@@ -392,8 +371,8 @@ namespace internal
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
       using VectorizedArrayType =
-        typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
-          Number>::value_type;
+        typename dealii::internal::VectorizedArrayTrait<
+          Number>::vectorized_value_type;
       using value_type               = Tensor<1, dim, Number>;
       using scalar_value_type        = Tensor<1, dim, ScalarNumber>;
       using vectorized_value_type    = Tensor<1, dim, VectorizedArrayType>;
@@ -404,17 +383,17 @@ namespace internal
         Tensor<1, dim, Tensor<1, dim, VectorizedArrayType>>;
 
       static void
-      read_value(const Number       vector_entry,
+      read_value(const ScalarNumber vector_entry,
                  const unsigned int component,
-                 value_type &       result)
+                 scalar_value_type &result)
       {
         result[component] = vector_entry;
       }
 
       static void
-      write_value(Number &           vector_entry,
-                  const unsigned int component,
-                  const value_type & result)
+      write_value(VectorizedArrayType &        vector_entry,
+                  const unsigned int           component,
+                  const vectorized_value_type &result)
       {
         vector_entry = result[component];
       }
@@ -422,41 +401,24 @@ namespace internal
       static void
       set_gradient(const interface_vectorized_gradient_type &value,
                    const unsigned int                        vector_lane,
-                   scalar_gradient_type &                    result)
+                   gradient_type &                           result)
       {
         for (unsigned int i = 0; i < dim; ++i)
           for (unsigned int d = 0; d < dim; ++d)
-            result[i][d] = value[d][i][vector_lane];
-      }
-
-      static void
-      set_gradient(const interface_vectorized_gradient_type &value,
-                   const unsigned int,
-                   vectorized_gradient_type &result)
-      {
-        for (unsigned int i = 0; i < dim; ++i)
-          for (unsigned int d = 0; d < dim; ++d)
-            result[i][d] = value[d][i];
+            result[i][d] =
+              internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+                value[d][i], vector_lane);
       }
 
       static void
       get_gradient(interface_vectorized_gradient_type &value,
                    const unsigned int                  vector_lane,
-                   const scalar_gradient_type &        result)
-      {
-        for (unsigned int i = 0; i < dim; ++i)
-          for (unsigned int d = 0; d < dim; ++d)
-            value[d][i][vector_lane] = result[i][d];
-      }
-
-      static void
-      get_gradient(interface_vectorized_gradient_type &value,
-                   const unsigned int,
-                   const vectorized_gradient_type &result)
+                   const gradient_type &               result)
       {
         for (unsigned int i = 0; i < dim; ++i)
           for (unsigned int d = 0; d < dim; ++d)
-            value[d][i] = result[i][d];
+            internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+              value[d][i], vector_lane) = result[i][d];
       }
 
       static void
@@ -561,8 +523,8 @@ namespace internal
       using ScalarNumber =
         typename internal::VectorizedArrayTrait<Number>::value_type;
       using VectorizedArrayType =
-        typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
-          Number>::value_type;
+        typename dealii::internal::VectorizedArrayTrait<
+          Number>::vectorized_value_type;
       using value_type               = Number;
       using scalar_value_type        = ScalarNumber;
       using vectorized_value_type    = VectorizedArrayType;
@@ -571,17 +533,17 @@ namespace internal
       using vectorized_gradient_type = Tensor<1, 1, VectorizedArrayType>;
 
       static void
-      read_value(const Number vector_entry,
+      read_value(const ScalarNumber vector_entry,
                  const unsigned int,
-                 value_type &result)
+                 scalar_value_type &result)
       {
         result = vector_entry;
       }
 
       static void
-      write_value(Number &vector_entry,
+      write_value(VectorizedArrayType &vector_entry,
                   const unsigned int,
-                  const value_type &result)
+                  const vectorized_value_type &result)
       {
         vector_entry = result;
       }
@@ -763,9 +725,8 @@ public:
 
   using ScalarNumber =
     typename internal::VectorizedArrayTrait<Number>::value_type;
-  using VectorizedArrayType =
-    typename internal::FEPointEvaluation::VectorizedPointsTypeTraits<
-      Number>::value_type;
+  using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+    Number>::vectorized_value_type;
   using ETT = typename internal::FEPointEvaluation::
     EvaluatorTypeTraits<dim, n_components, Number>;
   using value_type            = typename ETT::value_type;
@@ -812,9 +773,10 @@ public:
    * objects, this parameter allows to select a range of `n_components`
    * components starting from this parameter.
    */
-  FEPointEvaluation(NonMatching::MappingInfo<dim, spacedim> &mapping_info,
-                    const FiniteElement<dim> &               fe,
-                    const unsigned int first_selected_component = 0);
+  FEPointEvaluation(
+    NonMatching::MappingInfo<dim, spacedim, Number> &mapping_info,
+    const FiniteElement<dim> &                       fe,
+    const unsigned int first_selected_component = 0);
 
   /**
    * Copy constructor.
@@ -1021,6 +983,13 @@ public:
   quadrature_point_indices() const;
 
 private:
+  static constexpr std::size_t n_lanes_user_interface =
+    internal::VectorizedArrayTrait<Number>::width;
+  static constexpr std::size_t n_lanes_internal =
+    internal::VectorizedArrayTrait<VectorizedArrayType>::width;
+  static constexpr std::size_t stride =
+    internal::VectorizedArrayTrait<Number>::stride;
+
   /**
    * Common setup function for both constructors. Does the setup for both fast
    * and slow path.
@@ -1068,15 +1037,16 @@ private:
   integrate_slow(const ArrayView<ScalarNumber> &         solution_values,
                  const EvaluationFlags::EvaluationFlags &integration_flags);
 
+
   /**
-   * Number of quadrature points of the current cell/face.
+   * Number of quadrature batches of the current cell/face.
    */
   const unsigned int n_q_points;
 
   /**
-   * Number of active quadrature points of the last quadrature point batch.
+   * Number of quadrature points of the current cell/face.
    */
-  const unsigned int n_filled_lanes_last_batch;
+  const unsigned int n_q_points_scalar;
 
   /**
    * Pointer to the Mapping object passed to the constructor.
@@ -1137,6 +1107,43 @@ private:
    */
   std::vector<gradient_type> gradients;
 
+  /**
+   * Pointer to first unit point batch of current cell/face from MappingInfo,
+   * set internally during do_reinit().
+   */
+  const Point<dim, VectorizedArrayType> *unit_point_ptr;
+
+  /**
+   * Pointer to real point of first quadrature point of current cell/face from
+   * MappingInfo, set internally during do_reinit().
+   */
+  const Point<spacedim, Number> *real_point_ptr;
+
+  /**
+   * Pointer to Jacobian of first quadrature point of current cell/face from
+   * MappingInfo, set internally during do_reinit().
+   */
+  const DerivativeForm<1, dim, spacedim, Number> *jacobian_ptr;
+
+  /**
+   * Pointer to inverse Jacobian of first quadrature point of current cell/face
+   * from MappingInfo, set internally during do_reinit().
+   */
+  const DerivativeForm<1, spacedim, dim, Number> *inverse_jacobian_ptr;
+
+  /**
+   * Pointer to normal vector of first quadrature point of current cell/face
+   * from MappingInfo, set internally during do_reinit().
+   */
+  const Tensor<1, spacedim, Number> *normal_ptr;
+
+  /**
+   * Pointer to Jacobian determinant times quadrature weight of first quadrature
+   * point of current cell/face from MappingInfo, set internally during
+   * do_reinit().
+   */
+  const Number *JxW_ptr;
+
   /**
    * Number of unknowns per component, i.e., number of unique basis functions,
    * for the chosen FiniteElement (or base element).
@@ -1168,14 +1175,14 @@ private:
   /**
    * Pointer to mapping info on the fly computed during reinit.
    */
-  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim>>
+  std::unique_ptr<NonMatching::MappingInfo<dim, spacedim, Number>>
     mapping_info_on_the_fly;
 
   /**
    * Pointer to currently used mapping info (either on the fly or external
    * precomputed).
    */
-  SmartPointer<NonMatching::MappingInfo<dim, spacedim>> mapping_info;
+  SmartPointer<NonMatching::MappingInfo<dim, spacedim, Number>> mapping_info;
 
   /**
    * The current cell index to access mapping data from mapping info.
@@ -1220,13 +1227,14 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
   const UpdateFlags         update_flags,
   const unsigned int        first_selected_component)
   : n_q_points(numbers::invalid_unsigned_int)
-  , n_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(numbers::invalid_unsigned_int)
   , mapping(&mapping)
   , fe(&fe)
   , update_flags(update_flags)
   , mapping_info_on_the_fly(
-      std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(mapping,
-                                                                update_flags))
+      std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
+        mapping,
+        update_flags))
   , mapping_info(mapping_info_on_the_fly.get())
   , current_cell_index(numbers::invalid_unsigned_int)
   , current_face_number(numbers::invalid_unsigned_int)
@@ -1239,11 +1247,11 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
 
 template <int n_components, int dim, int spacedim, typename Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::FEPointEvaluation(
-  NonMatching::MappingInfo<dim, spacedim> &mapping_info,
-  const FiniteElement<dim> &               fe,
-  const unsigned int                       first_selected_component)
+  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_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(numbers::invalid_unsigned_int)
   , mapping(&mapping_info.get_mapping())
   , fe(&fe)
   , update_flags(mapping_info.get_update_flags())
@@ -1263,7 +1271,7 @@ 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_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(other.n_q_points_scalar)
   , mapping(other.mapping)
   , fe(other.fe)
   , poly(other.poly)
@@ -1281,7 +1289,7 @@ FEPointEvaluation<n_components_, dim, spacedim, Number>::FEPointEvaluation(
   , fe_values(other.fe_values)
   , mapping_info_on_the_fly(
       other.mapping_info_on_the_fly ?
-        std::make_unique<NonMatching::MappingInfo<dim, spacedim>>(
+        std::make_unique<NonMatching::MappingInfo<dim, spacedim, Number>>(
           *mapping,
           update_flags) :
         nullptr)
@@ -1302,7 +1310,7 @@ 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_filled_lanes_last_batch(numbers::invalid_unsigned_int)
+  , n_q_points_scalar(other.n_q_points_scalar)
   , mapping(other.mapping)
   , fe(other.fe)
   , poly(other.poly)
@@ -1482,54 +1490,55 @@ template <int n_components, int dim, int spacedim, typename Number>
 void
 FEPointEvaluation<n_components, dim, spacedim, Number>::do_reinit()
 {
-  const auto unit_points =
-    mapping_info->get_unit_points(current_cell_index, current_face_number);
-
-  const unsigned int n_q_points_unvectorized = unit_points.size();
-
-  if (std::is_same<ScalarNumber, Number>::value)
-    {
-      const_cast<unsigned int &>(n_q_points) = unit_points.size();
-    }
-  else
-    {
-      const unsigned int n_lanes =
-        internal::VectorizedArrayTrait<VectorizedArrayType>::width;
-      const_cast<unsigned int &>(n_filled_lanes_last_batch) =
-        n_q_points_unvectorized % n_lanes;
-      const_cast<unsigned int &>(n_q_points) =
-        n_q_points_unvectorized / n_lanes;
-      if (n_filled_lanes_last_batch > 0)
-        ++const_cast<unsigned int &>(n_q_points);
-    }
+  const_cast<unsigned int &>(n_q_points_scalar) =
+    mapping_info->get_n_q_points_unvectorized(current_cell_index,
+                                              current_face_number);
+
+  const_cast<unsigned int &>(n_q_points) =
+    n_q_points_scalar / n_lanes_user_interface +
+    (n_q_points_scalar % n_lanes_user_interface > 0 ? 1 : 0);
+
+  // set unit point pointer
+  const unsigned int unit_point_offset =
+    mapping_info->compute_unit_point_index_offset(current_cell_index,
+                                                  current_face_number);
+  unit_point_ptr = mapping_info->get_unit_point(unit_point_offset);
+
+  // set data pointers
+  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);
+  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)
+    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);
 
   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 (!polynomials_are_hat_functions)
+  if (fast_path && !polynomials_are_hat_functions)
     {
-      const std::size_t n_points = unit_points.size();
-      const std::size_t n_lanes  = VectorizedArrayType::size();
       const std::size_t n_batches =
-        n_points / n_lanes + (n_points % n_lanes > 0 ? 1 : 0);
+        n_q_points_scalar / n_lanes_internal +
+        (n_q_points_scalar % n_lanes_internal > 0 ? 1 : 0);
       const std::size_t n_shapes = poly.size();
       shapes.resize_fast(n_batches * n_shapes);
-      for (unsigned int i = 0, qb = 0; i < n_points; i += n_lanes, ++qb)
-        {
-          // convert to vectorized format
-          Point<dim, VectorizedArrayType> vectorized_points;
-          for (unsigned int j = 0; j < n_lanes && i + j < n_points; ++j)
-            for (unsigned int d = 0; d < dim; ++d)
-              vectorized_points[d][j] = unit_points[i + j][d];
-
-          auto view =
-            make_array_view(shapes.begin() + qb * n_shapes,
-                            shapes.begin() + (qb * n_shapes + n_shapes));
-
-          internal::compute_values_of_array(view, poly, vectorized_points);
-        }
+      for (unsigned int qb = 0; qb < n_batches; ++qb)
+        internal::compute_values_of_array(make_array_view(shapes,
+                                                          qb * n_shapes,
+                                                          n_shapes),
+                                          poly,
+                                          unit_point_ptr[qb]);
     }
 
   is_reinitialized = true;
@@ -1548,103 +1557,54 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate_fast(
     solution_renumbered.resize(dofs_per_component);
   for (unsigned int comp = 0; comp < n_components; ++comp)
     for (unsigned int i = 0; i < dofs_per_component; ++i)
-      internal::FEPointEvaluation::
-        EvaluatorTypeTraits<dim, n_components, ScalarNumber>::read_value(
-          solution_values[renumber[(component_in_base_element + comp) *
-                                     dofs_per_component +
-                                   i]],
-          comp,
-          solution_renumbered[i]);
+      ETT::read_value(
+        solution_values[renumber[(component_in_base_element + comp) *
+                                   dofs_per_component +
+                                 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>());
 
-  const auto unit_points =
-    mapping_info->get_unit_points(current_cell_index, current_face_number);
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const std::size_t n_points = unit_points.size();
-  const std::size_t n_lanes =
-    internal::VectorizedArrayTrait<VectorizedArrayType>::width;
-  const std::size_t stride =
-    internal::FEPointEvaluation::VectorizedPointsTypeTraits<Number>::stride;
-
   // loop over quadrature batches qb / points q
-  for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+  const unsigned int n_shapes = poly.size();
+  for (unsigned int qb = 0, q = 0; q < n_q_points_scalar;
+       ++qb, q += n_lanes_internal)
     {
       // compute
-      const unsigned int n_shapes     = poly.size();
-      const auto         val_and_grad = [&]() {
-        if (polynomials_are_hat_functions)
-          {
-            // convert quadrature points to vectorized format
-            Point<dim, VectorizedArrayType> vectorized_points;
-            for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
-              for (unsigned int d = 0; d < dim; ++d)
-                vectorized_points[d][v] = unit_points[q + v][d];
-
-            return internal::evaluate_tensor_product_value_and_gradient_linear(
-              poly, solution_renumbered, vectorized_points);
-          }
-        else
-          return internal::evaluate_tensor_product_value_and_gradient_shapes<
+      const auto val_and_grad =
+        polynomials_are_hat_functions ?
+          internal::evaluate_tensor_product_value_and_gradient_linear(
+            poly, solution_renumbered, unit_point_ptr[qb]) :
+          internal::evaluate_tensor_product_value_and_gradient_shapes<
             dim,
             scalar_value_type,
-            VectorizedArrayType>(make_array_view(shapes.begin() + qb * n_shapes,
-                                                 shapes.begin() +
-                                                   (qb * n_shapes + n_shapes)),
-                                 poly.size(),
+            VectorizedArrayType>(make_array_view(shapes,
+                                                 qb * n_shapes,
+                                                 n_shapes),
+                                 n_shapes,
                                  solution_renumbered);
-      }();
 
       if (evaluation_flags & EvaluationFlags::values)
         {
-          for (unsigned int v = 0; v < stride && q + v < n_points; ++v)
-            internal::FEPointEvaluation::
-              EvaluatorTypeTraits<dim, n_components, Number>::set_value(
-                val_and_grad.first, v, values[qb * stride + v]);
+          for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
+            ETT::set_value(val_and_grad.first, v, values[qb * stride + v]);
         }
       if (evaluation_flags & EvaluationFlags::gradients)
         {
           Assert(update_flags & update_gradients ||
                    update_flags & update_inverse_jacobians,
                  ExcNotInitialized());
-          if (std::is_same<Number, ScalarNumber>::value)
-            {
-              // convert back to standard format
-              for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
-                {
-                  internal::FEPointEvaluation::EvaluatorTypeTraits<
-                    dim,
-                    n_components,
-                    Number>::set_gradient(val_and_grad.second,
-                                          v,
-                                          unit_gradients[q + v]);
-                  gradients[q + v] = apply_transformation(
-                    mapping_data.inverse_jacobians[q + v].transpose(),
-                    unit_gradients[q + v]);
-                }
-            }
-          else
+
+          for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
             {
-              // convert to vectorized format
-              DerivativeForm<1, spacedim, dim, Number>
-                vectorized_inverse_transposed_jacobians;
-              for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
-                for (unsigned int d = 0; d < dim; ++d)
-                  for (unsigned int s = 0; s < spacedim; ++s)
-                    internal::VectorizedArrayTrait<Number>::get(
-                      vectorized_inverse_transposed_jacobians[s][d], v) =
-                      mapping_data.inverse_jacobians[q + v].transpose()[s][d];
-
-              internal::FEPointEvaluation::
-                EvaluatorTypeTraits<dim, n_components, Number>::set_gradient(
-                  val_and_grad.second, 0, unit_gradients[qb]);
-              gradients[qb] =
-                apply_transformation(vectorized_inverse_transposed_jacobians,
-                                     unit_gradients[qb]);
+              const unsigned int offset = qb * stride + v;
+              ETT::set_gradient(val_and_grad.second, v, unit_gradients[offset]);
+              gradients[offset] =
+                apply_transformation(inverse_jacobian_ptr[offset].transpose(),
+                                     unit_gradients[offset]);
             }
         }
     }
@@ -1664,7 +1624,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate_slow(
            "Not initialized. Please call FEPointEvaluation::reinit()!"));
 
   const std::size_t n_points = fe_values->get_quadrature().size();
-  const std::size_t n_lanes  = internal::VectorizedArrayTrait<Number>::width;
 
   if (evaluation_flags & EvaluationFlags::values)
     {
@@ -1676,29 +1635,26 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate_slow(
           for (unsigned int d = 0; d < n_components; ++d)
             if (nonzero_shape_function_component[i][d] &&
                 (fe->is_primitive(i) || fe->is_primitive()))
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  internal::FEPointEvaluation::
-                    EvaluatorTypeTraits<dim, n_components, Number>::access(
-                      values[qb],
-                      v,
-                      d,
-                      fe_values->shape_value(i, q + v) * value);
+                  ETT::access(values[qb],
+                              v,
+                              d,
+                              fe_values->shape_value(i, q + v) * value);
             else if (nonzero_shape_function_component[i][d])
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  internal::FEPointEvaluation::
-                    EvaluatorTypeTraits<dim, n_components, Number>::access(
-                      values[qb],
-                      v,
-                      d,
-                      fe_values->shape_value_component(i, q + v, d) * value);
+                  ETT::access(values[qb],
+                              v,
+                              d,
+                              fe_values->shape_value_component(i, q + v, d) *
+                                value);
         }
     }
 
@@ -1712,29 +1668,26 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::evaluate_slow(
           for (unsigned int d = 0; d < n_components; ++d)
             if (nonzero_shape_function_component[i][d] &&
                 (fe->is_primitive(i) || fe->is_primitive()))
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  internal::FEPointEvaluation::
-                    EvaluatorTypeTraits<dim, n_components, Number>::access(
-                      gradients[qb],
-                      v,
-                      d,
-                      fe_values->shape_grad(i, q + v) * value);
+                  ETT::access(gradients[qb],
+                              v,
+                              d,
+                              fe_values->shape_grad(i, q + v) * value);
             else if (nonzero_shape_function_component[i][d])
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  internal::FEPointEvaluation::
-                    EvaluatorTypeTraits<dim, n_components, Number>::access(
-                      gradients[qb],
-                      v,
-                      d,
-                      fe_values->shape_grad_component(i, q + v, d) * value);
+                  ETT::access(gradients[qb],
+                              v,
+                              d,
+                              fe_values->shape_grad_component(i, q + v, d) *
+                                value);
         }
     }
 }
@@ -1773,135 +1726,83 @@ inline void
 FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_fast(
   const ArrayView<ScalarNumber> &         solution_values,
   const EvaluationFlags::EvaluationFlags &integration_flags)
-
 {
   // fast path with tensor product integration
   if (solution_renumbered_vectorized.size() != dofs_per_component)
     solution_renumbered_vectorized.resize(dofs_per_component);
   // zero content
-  solution_renumbered_vectorized.fill(
-    typename internal::FEPointEvaluation::EvaluatorTypeTraits<
-      dim,
-      n_components,
-      VectorizedArrayType>::value_type());
-
-  const auto unit_points =
-    mapping_info->get_unit_points(current_cell_index, current_face_number);
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const std::size_t n_points = unit_points.size();
-  const std::size_t n_lanes =
-    internal::VectorizedArrayTrait<VectorizedArrayType>::width;
-  const std::size_t stride =
-    internal::FEPointEvaluation::VectorizedPointsTypeTraits<Number>::stride;
+  solution_renumbered_vectorized.fill(vectorized_value_type());
 
   // loop over quadrature batches qb / points q
-  for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+  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 =
-        (qb == (n_q_points - 1)) && (n_filled_lanes_last_batch > 0);
+        q + n_lanes_user_interface > n_q_points_scalar;
 
-      typename internal::ProductTypeNoPoint<value_type,
-                                            VectorizedArrayType>::type value =
-        {};
-      Tensor<1,
-             dim,
-             typename internal::ProductTypeNoPoint<value_type,
-                                                   VectorizedArrayType>::type>
-        gradient;
+      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)
             {
-              for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes; ++v)
-                internal::FEPointEvaluation::EvaluatorTypeTraits<
-                  dim,
-                  n_components,
-                  Number>::set_zero_value(values[qb], v);
+              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; v < stride && q + v < n_points; ++v)
-            internal::FEPointEvaluation::
-              EvaluatorTypeTraits<dim, n_components, Number>::get_value(
-                value, v, values[qb * stride + v]);
+
+          for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
+            ETT::get_value(value, v, values[qb * stride + v]);
         }
       if (integration_flags & EvaluationFlags::gradients)
         {
-          if (std::is_same<Number, ScalarNumber>::value)
+          // zero out lanes of incomplete last quadrature point batch
+          if (incomplete_last_batch)
             {
-              for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
-                {
-                  gradients[q + v] =
-                    apply_transformation(mapping_data.inverse_jacobians[q + v],
-                                         gradients[q + v]);
-                  internal::FEPointEvaluation::EvaluatorTypeTraits<
-                    dim,
-                    n_components,
-                    Number>::get_gradient(gradient, v, gradients[q + v]);
-                }
+              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);
             }
-          else
+
+          for (unsigned int v = 0; v < stride && q + v < n_q_points_scalar; ++v)
             {
-              // convert to vectorized format
-              DerivativeForm<1, spacedim, dim, Number>
-                vectorized_inverse_jacobians;
-              for (unsigned int v = 0; v < n_lanes && q + v < n_points; ++v)
-                for (unsigned int d = 0; d < dim; ++d)
-                  for (unsigned int s = 0; s < spacedim; ++s)
-                    internal::VectorizedArrayTrait<Number>::get(
-                      vectorized_inverse_jacobians[s][d], v) =
-                      mapping_data.inverse_jacobians[q + v][s][d];
-
-              // zero out lanes of incomplete last quadrature point batch
-              if (incomplete_last_batch)
-                {
-                  for (unsigned int v = n_filled_lanes_last_batch; v < n_lanes;
-                       ++v)
-                    internal::FEPointEvaluation::EvaluatorTypeTraits<
-                      dim,
-                      n_components,
-                      Number>::set_zero_gradient(gradients[qb], v);
-                }
-
-              gradients[qb] = apply_transformation(vectorized_inverse_jacobians,
-                                                   gradients[qb]);
-              internal::FEPointEvaluation::
-                EvaluatorTypeTraits<dim, n_components, Number>::get_gradient(
-                  gradient, 0, gradients[qb]);
+              const unsigned int offset = qb * stride + v;
+              ETT::get_gradient(
+                gradient,
+                v,
+                apply_transformation(inverse_jacobian_ptr[offset],
+                                     gradients[offset]));
             }
         }
 
       // compute
-      const unsigned int n_shapes = poly.size();
       if (polynomials_are_hat_functions)
-        {
-          // convert to vectorized format
-          Point<dim, VectorizedArrayType> vectorized_points;
-          for (unsigned int j = 0; j < n_lanes && q + j < n_points; ++j)
-            for (unsigned int d = 0; d < dim; ++d)
-              vectorized_points[d][j] = unit_points[q + j][d];
-
-          internal::integrate_add_tensor_product_value_and_gradient_linear(
-            poly,
-            value,
-            gradient,
-            solution_renumbered_vectorized,
-            vectorized_points);
-        }
+        internal::integrate_add_tensor_product_value_and_gradient_linear(
+          poly,
+          value,
+          gradient,
+          solution_renumbered_vectorized,
+          unit_point_ptr[qb]);
       else
         internal::integrate_add_tensor_product_value_and_gradient_shapes<
           dim,
           VectorizedArrayType,
-          typename internal::ProductTypeNoPoint<value_type,
-                                                VectorizedArrayType>::type>(
-          make_array_view(shapes.begin() + qb * n_shapes,
-                          shapes.begin() + (qb * n_shapes + n_shapes)),
-          n_shapes,
-          value,
-          gradient,
-          solution_renumbered_vectorized);
+          vectorized_value_type>(make_array_view(shapes,
+                                                 qb * n_shapes,
+                                                 n_shapes),
+                                 n_shapes,
+                                 value,
+                                 gradient,
+                                 solution_renumbered_vectorized);
     }
 
   // add between the lanes and write into the result
@@ -1910,13 +1811,8 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_fast(
     for (unsigned int i = 0; i < dofs_per_component; ++i)
       {
         VectorizedArrayType result;
-        internal::FEPointEvaluation::EvaluatorTypeTraits<
-          dim,
-          n_components,
-          VectorizedArrayType>::write_value(result,
-                                            comp,
-                                            solution_renumbered_vectorized[i]);
-        for (unsigned int lane = n_lanes / 2; lane > 0; lane /= 2)
+        ETT::write_value(result, comp, solution_renumbered_vectorized[i]);
+        for (unsigned int lane = n_lanes_internal / 2; lane > 0; lane /= 2)
           for (unsigned int j = 0; j < lane; ++j)
             result[j] += result[lane + j];
         solution_values[renumber[comp * dofs_per_component + i]] = result[0];
@@ -1938,7 +1834,6 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_slow(
   std::fill(solution_values.begin(), solution_values.end(), 0.0);
 
   const std::size_t n_points = fe_values->get_quadrature().size();
-  const std::size_t n_lanes  = internal::VectorizedArrayTrait<Number>::width;
 
   if (integration_flags & EvaluationFlags::values)
     {
@@ -1948,27 +1843,22 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_slow(
           for (unsigned int d = 0; d < n_components; ++d)
             if (nonzero_shape_function_component[i][d] &&
                 (fe->is_primitive(i) || fe->is_primitive()))
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  solution_values[i] +=
-                    fe_values->shape_value(i, q + v) *
-                    internal::FEPointEvaluation::
-                      EvaluatorTypeTraits<dim, n_components, Number>::access(
-                        values[qb], v, d);
+                  solution_values[i] += fe_values->shape_value(i, q + v) *
+                                        ETT::access(values[qb], v, d);
             else if (nonzero_shape_function_component[i][d])
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
                   solution_values[i] +=
                     fe_values->shape_value_component(i, q + v, d) *
-                    internal::FEPointEvaluation::
-                      EvaluatorTypeTraits<dim, n_components, Number>::access(
-                        values[qb], v, d);
+                    ETT::access(values[qb], v, d);
         }
     }
 
@@ -1980,27 +1870,22 @@ FEPointEvaluation<n_components, dim, spacedim, Number>::integrate_slow(
           for (unsigned int d = 0; d < n_components; ++d)
             if (nonzero_shape_function_component[i][d] &&
                 (fe->is_primitive(i) || fe->is_primitive()))
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
-                  solution_values[i] +=
-                    fe_values->shape_grad(i, q + v) *
-                    internal::FEPointEvaluation::
-                      EvaluatorTypeTraits<dim, n_components, Number>::access(
-                        gradients[qb], v, d);
+                  solution_values[i] += fe_values->shape_grad(i, q + v) *
+                                        ETT::access(gradients[qb], v, d);
             else if (nonzero_shape_function_component[i][d])
-              for (unsigned int qb = 0, q = 0; q < n_points; ++qb, q += n_lanes)
+              for (unsigned int qb = 0, q = 0; q < n_points;
+                   ++qb, q += n_lanes_user_interface)
                 for (unsigned int v = 0;
-                     v < (q + n_lanes > n_points ? n_filled_lanes_last_batch :
-                                                   n_lanes);
+                     v < n_lanes_user_interface && q + v < n_points;
                      ++v)
                   solution_values[i] +=
                     fe_values->shape_grad_component(i, q + v, d) *
-                    internal::FEPointEvaluation::
-                      EvaluatorTypeTraits<dim, n_components, Number>::access(
-                        gradients[qb], v, d);
+                    ETT::access(gradients[qb], v, d);
         }
     }
 }
@@ -2110,20 +1995,12 @@ inline DerivativeForm<1, dim, spacedim, Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::jacobian(
   const unsigned int point_index) const
 {
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  DerivativeForm<1, dim, spacedim, Number> vectorized_jacobian;
-  for (unsigned int v = 0;
-       v < n_lanes && point_index * n_lanes + v < mapping_data.jacobians.size();
-       ++v)
-    for (unsigned int d = 0; d < dim; ++d)
-      for (unsigned int s = 0; s < spacedim; ++s)
-        internal::VectorizedArrayTrait<Number>::get(vectorized_jacobian[d][s],
-                                                    v) =
-          mapping_data.jacobians[point_index * n_lanes + v][d][s];
-  return vectorized_jacobian;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(jacobian_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_jacobians"));
+  return jacobian_ptr[point_index];
 }
 
 
@@ -2133,21 +2010,12 @@ inline DerivativeForm<1, spacedim, dim, Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::inverse_jacobian(
   const unsigned int point_index) const
 {
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  DerivativeForm<1, dim, spacedim, Number> vectorized_inverse_jacobian;
-  for (unsigned int v = 0;
-       v < n_lanes &&
-       point_index * n_lanes + v < mapping_data.inverse_jacobians.size();
-       ++v)
-    for (unsigned int d = 0; d < dim; ++d)
-      for (unsigned int s = 0; s < spacedim; ++s)
-        internal::VectorizedArrayTrait<Number>::get(
-          vectorized_inverse_jacobian[d][s], v) =
-          mapping_data.inverse_jacobians[point_index * n_lanes + v][d][s];
-  return vectorized_inverse_jacobian;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(inverse_jacobian_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_inverse_jacobians"));
+  return inverse_jacobian_ptr[point_index];
 }
 
 
@@ -2157,38 +2025,27 @@ inline Number
 FEPointEvaluation<n_components, dim, spacedim, Number>::JxW(
   const unsigned int point_index) const
 {
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  Number             vectorized_JxW = 0;
-  for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v <
-                                            mapping_data.JxW_values.size();
-       ++v)
-    internal::VectorizedArrayTrait<Number>::get(vectorized_JxW, v) =
-      mapping_data.JxW_values[point_index * n_lanes + v];
-  return vectorized_JxW;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(JxW_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_JxW_values"));
+  return JxW_ptr[point_index];
 }
 
 
+
 template <int n_components, int dim, int spacedim, typename Number>
 inline Tensor<1, spacedim, Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::normal_vector(
   const unsigned int point_index) const
 {
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  Tensor<1, spacedim, Number> vectorized_normal_vectors;
-  for (unsigned int v = 0; v < n_lanes && point_index * n_lanes + v <
-                                            mapping_data.normal_vectors.size();
-       ++v)
-    for (unsigned int s = 0; s < spacedim; ++s)
-      internal::VectorizedArrayTrait<Number>::get(vectorized_normal_vectors[s],
-                                                  v) =
-        mapping_data.normal_vectors[point_index * n_lanes + v][s];
-  return vectorized_normal_vectors;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(normal_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_normal_vectors"));
+  return normal_ptr[point_index];
 }
 
 
@@ -2198,19 +2055,12 @@ inline Point<spacedim, Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::real_point(
   const unsigned int point_index) const
 {
-  const auto &mapping_data =
-    mapping_info->get_mapping_data(current_cell_index, current_face_number);
-
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  Point<spacedim, Number> vectorized_real_point;
-  for (unsigned int v = 0;
-       v < n_lanes &&
-       point_index * n_lanes + v < mapping_data.quadrature_points.size();
-       ++v)
-    for (unsigned int s = 0; s < spacedim; ++s)
-      internal::VectorizedArrayTrait<Number>::get(vectorized_real_point[s], v) =
-        mapping_data.quadrature_points[point_index * n_lanes + v][s];
-  return vectorized_real_point;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(real_point_ptr != nullptr,
+         internal::FEPointEvaluation::
+           ExcFEPointEvaluationAccessToUninitializedMappingField(
+             "update_quadrature_points"));
+  return real_point_ptr[point_index];
 }
 
 
@@ -2220,17 +2070,13 @@ inline Point<dim, Number>
 FEPointEvaluation<n_components, dim, spacedim, Number>::unit_point(
   const unsigned int point_index) const
 {
-  const auto unit_points =
-    mapping_info->get_unit_points(current_cell_index, current_face_number);
-  const unsigned int n_lanes = internal::VectorizedArrayTrait<Number>::width;
-  Point<spacedim, Number> vectorized_unit_point;
-  for (unsigned int v = 0;
-       v < n_lanes && point_index * n_lanes + v < unit_points.size();
-       ++v)
-    for (unsigned int d = 0; d < dim; ++d)
-      internal::VectorizedArrayTrait<Number>::get(vectorized_unit_point[d], v) =
-        unit_points[point_index * n_lanes + v][d];
-  return vectorized_unit_point;
+  AssertIndexRange(point_index, n_q_points);
+  Assert(unit_point_ptr != nullptr, ExcMessage("unit_point_ptr is not set!"));
+  Point<dim, Number> unit_point;
+  for (unsigned int d = 0; d < dim; ++d)
+    unit_point[d] = internal::VectorizedArrayTrait<Number>::get_from_vectorized(
+      unit_point_ptr[point_index / stride][d], point_index % stride);
+  return unit_point;
 }
 
 
index 29af641639a368f04fc8a0b09e6920e246dd47e4..99b69eb97f0688fb2c67987137f2e655484c4006 100644 (file)
@@ -2989,7 +2989,7 @@ namespace internal
   template <int dim, typename Number>
   inline void
   compute_values_of_array(
-    ArrayView<dealii::ndarray<Number, 2, dim>> &        shapes,
+    ArrayView<dealii::ndarray<Number, 2, dim>>          shapes,
     const std::vector<Polynomials::Polynomial<double>> &poly,
     const Point<dim, Number> &                          p)
   {
index 7a106c18ace0bdb2129f9c2d32d3f2e8a2b7ea18..418cd25f7808555609dc92bef86229dd757b53b7 100644 (file)
@@ -38,17 +38,174 @@ DEAL_II_NAMESPACE_OPEN
 
 namespace NonMatching
 {
+  namespace internal
+  {
+    template <int dim, int spacedim = dim>
+    class ComputeMappingDataHelper
+    {
+      using MappingData =
+        dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+                                                                     spacedim>;
+
+    public:
+      static UpdateFlags
+      required_update_flags(
+        const SmartPointer<const Mapping<dim, spacedim>> mapping,
+        const UpdateFlags &                              update_flags)
+      {
+        return mapping->requires_update_flags(update_flags);
+      }
+
+      static void
+      compute_mapping_data_for_quadrature(
+        const SmartPointer<const Mapping<dim, spacedim>> mapping,
+        const UpdateFlags &                              update_flags_mapping,
+        const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+        CellSimilarity::Similarity &cell_similarity,
+        const Quadrature<dim> &     quadrature,
+        std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+                     internal_mapping_data,
+        MappingData &mapping_data)
+      {
+        mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+        // reuse internal_mapping_data for MappingQ to avoid memory allocations
+        if (const MappingQ<dim, spacedim> *mapping_q =
+              dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+          {
+            (void)mapping_q;
+            auto &data =
+              dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+                *internal_mapping_data);
+            data.initialize(update_flags_mapping,
+                            quadrature,
+                            quadrature.size());
+          }
+        else
+          {
+            internal_mapping_data =
+              mapping->get_data(update_flags_mapping, quadrature);
+          }
+
+        cell_similarity = mapping->fill_fe_values(cell,
+                                                  cell_similarity,
+                                                  quadrature,
+                                                  *internal_mapping_data,
+                                                  mapping_data);
+      }
+
+
+
+      static void
+      compute_mapping_data_for_immersed_surface_quadrature(
+        const SmartPointer<const Mapping<dim, spacedim>> mapping,
+        const UpdateFlags &                              update_flags_mapping,
+        const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+        const ImmersedSurfaceQuadrature<dim> &                      quadrature,
+        std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+                     internal_mapping_data,
+        MappingData &mapping_data)
+      {
+        mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+        // reuse internal_mapping_data for MappingQ to avoid memory allocations
+        if (dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+          {
+            auto &data =
+              dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+                *internal_mapping_data);
+            data.initialize(update_flags_mapping,
+                            quadrature,
+                            quadrature.size());
+          }
+        else
+          {
+            internal_mapping_data =
+              mapping->get_data(update_flags_mapping, quadrature);
+          }
+
+        mapping->fill_fe_immersed_surface_values(cell,
+                                                 quadrature,
+                                                 *internal_mapping_data,
+                                                 mapping_data);
+      }
+
+
+
+      static void
+      compute_mapping_data_for_face_quadrature(
+        const SmartPointer<const Mapping<dim, spacedim>> mapping,
+        const UpdateFlags &                              update_flags_mapping,
+        const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+        const unsigned int                                          face_no,
+        const Quadrature<dim - 1> &                                 quadrature,
+        std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+                     internal_mapping_data,
+        MappingData &mapping_data)
+      {
+        mapping_data.initialize(quadrature.size(), update_flags_mapping);
+
+        // reuse internal_mapping_data for MappingQ to avoid memory allocations
+        if (const MappingQ<dim, spacedim> *mapping_q =
+              dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+          {
+            auto &data =
+              dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
+                *internal_mapping_data);
+            data.initialize_face(update_flags_mapping,
+                                 QProjector<dim>::project_to_oriented_face(
+                                   ReferenceCells::get_hypercube<dim>(),
+                                   quadrature,
+                                   face_no,
+                                   cell->face_orientation(face_no),
+                                   cell->face_flip(face_no),
+                                   cell->face_rotation(face_no)),
+                                 quadrature.size());
+
+            mapping_q->fill_mapping_data_for_face_quadrature(
+              cell, face_no, quadrature, *internal_mapping_data, mapping_data);
+          }
+        else
+          {
+            auto internal_mapping_data =
+              mapping->get_face_data(update_flags_mapping,
+                                     hp::QCollection<dim - 1>(quadrature));
+
+            mapping->fill_fe_face_values(cell,
+                                         face_no,
+                                         hp::QCollection<dim - 1>(quadrature),
+                                         *internal_mapping_data,
+                                         mapping_data);
+          }
+      }
+    };
+  } // namespace internal
+
   /**
    * This class provides the mapping information computation and mapping data
    * storage to be used together with FEPointEvaluation.
+   *
+   * MappingInfo is vectorized across quadrature points, which means data can be
+   * provided in vectorized format.
+   *
+   * Two different modes are available: partially vectorized and fully
+   * vectorized across quadrature points. Partially vectorized (scalar Number
+   * template argument) means the computed mapping data is provided in scalar
+   * format and only unit points are provided vectorized (for seamless
+   * interaction with FEPointEvaluation). Fully vectorized (vectorized Number
+   * template argument, e.g. VectorizedArray<double>) provides both mapping data
+   * and unit points in vectorized format. The Number template parameter of
+   * MappingInfo and FEPointEvaluation has to be identical.
    */
-  template <int dim, int spacedim = dim>
+  template <int dim, int spacedim = dim, typename Number = double>
   class MappingInfo : public Subscriptor
   {
   public:
-    using MappingData =
-      dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
-                                                                   spacedim>;
+    /**
+     * The VectorizedArray type the unit points are stored in inside this class.
+     */
+    using VectorizedArrayType = typename dealii::internal::VectorizedArrayTrait<
+      Number>::vectorized_value_type;
 
     /**
      * Constructor.
@@ -139,37 +296,46 @@ namespace NonMatching
       const unsigned int n_unfiltered_cells = numbers::invalid_unsigned_int);
 
     /**
-     * Getter function for current unit points.
-     *
-     * @p cell_index and @p face_number are the indices
-     * into the compressed, CRS like data storage of unit points.
-     *
-     * If you have initialized this object with reinit_cells() you can access
-     * the stored unit points of the cell with the respective @p cell_index
-     * (and the default argument for @p face_number).
-     *
-     * If you have initialized this object with reinit_faces() you can access
-     * the stored unit points of the faces on the cell with the respective @p cell_index
-     * and the respective local @p face_number.
-     *
-     * If you have initialized this object with reinit() you can access the
-     * stored unit points of a single cell with the default arguments.
-     *
-     * The correct state of this object is checked in this call (in debug mode).
+     * Getter function for unit points. The offset can be obtained with
+     * compute_unit_point_index_offset().
+     */
+    const Point<dim, VectorizedArrayType> *
+    get_unit_point(const unsigned int offset) const;
+
+    /**
+     * Getter function for Jacobians. The offset can be obtained with
+     * compute_data_index_offset().
+     */
+    const DerivativeForm<1, dim, spacedim, Number> *
+    get_jacobian(const unsigned int offset) const;
+
+    /**
+     * Getter function for inverse Jacobians. The offset can be obtained with
+     * compute_data_index_offset().
      */
-    const ArrayView<const Point<dim>>
-    get_unit_points(
-      const unsigned int cell_index  = numbers::invalid_unsigned_int,
-      const unsigned int face_number = numbers::invalid_unsigned_int) const;
+    const DerivativeForm<1, spacedim, dim, Number> *
+    get_inverse_jacobian(const unsigned int offset) const;
 
     /**
-     * Getter function for computed mapping data. This function accesses
-     * internal data and is therefore not a stable interface.
+     * Getter function for normal vectors. The offset can be obtained with
+     * compute_data_index_offset().
      */
-    const MappingData &
-    get_mapping_data(
-      const unsigned int cell_index  = numbers::invalid_unsigned_int,
-      const unsigned int face_number = numbers::invalid_unsigned_int) const;
+    const Tensor<1, spacedim, Number> *
+    get_normal_vector(const unsigned int offset) const;
+
+    /**
+     * Getter function for Jacobian times quadrature weight (JxW). The offset
+     * can be obtained with compute_data_index_offset().
+     */
+    const Number *
+    get_JxW(const unsigned int offset) const;
+
+    /**
+     * Getter function for real points. The offset can be obtained with
+     * compute_data_index_offset().
+     */
+    const Point<dim, Number> *
+    get_real_point(const unsigned int offset) const;
 
     /**
      * Getter function for underlying mapping.
@@ -183,79 +349,127 @@ namespace NonMatching
     UpdateFlags
     get_update_flags() const;
 
+    /**
+     * Getter function for the mapping update flags.
+     */
+    UpdateFlags
+    get_update_flags_mapping() const;
+
     /**
      * Connects to is_reinitialized().
      */
     boost::signals2::connection
     connect_is_reinitialized(const std::function<void()> &set_is_reinitialized);
 
+    /**
+     * 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 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;
+
+    /**
+     * Get number of unvectorized quadrature points.
+     */
+    unsigned int
+    get_n_q_points_unvectorized(const unsigned int cell_index,
+                                const unsigned int face_number) const;
+
   private:
+    using MappingData =
+      dealii::internal::FEValuesImplementation::MappingRelatedData<dim,
+                                                                   spacedim>;
+
     /**
-     * Enum class for reinitialized states.
+     * Compute number of quadrature point batches depending on NumberType.
      */
-    enum class State
-    {
-      invalid,
-      single_cell,
-      cell_vector,
-      faces_on_cells_in_vector
-    };
+    template <typename NumberType>
+    unsigned int
+    compute_n_q_points(const unsigned int n_q_points_unvectorized);
 
     /**
-     * Enum class that stores the currently initialized state
-     * upon the last call of reinit().
+     * Resize the unit_point data field.
      */
-    State state;
+    void
+    resize_unit_points(const unsigned int n_unit_point_batches);
 
     /**
-     * Compute the mapping related data for the given @p mapping,
-     * @p cell and @p unit_points that is required by the FEPointEvaluation
-     * class.
+     * Resize the mapping data fields.
      */
     void
-    compute_mapping_data_for_quadrature(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      CellSimilarity::Similarity &cell_similarity,
-      const Quadrature<dim> &     quadrature,
-      MappingData &               mapping_data);
+    resize_data_fields(const unsigned int n_data_point_batches);
 
     /**
-     * Compute the mapping related data for the given @p mapping,
-     * @p cell and @p quadrature that is required by the FEPointEvaluation
-     * class.
+     * Store the unit points.
      */
     void
-    compute_mapping_data_for_immersed_surface_quadrature(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const ImmersedSurfaceQuadrature<dim> &                      quadrature,
-      MappingData &                                               mapping_data);
+    store_unit_points(const unsigned int             unit_points_index_offset,
+                      const unsigned int             n_q_points,
+                      const unsigned int             n_q_points_unvectorized,
+                      const std::vector<Point<dim>> &points);
 
     /**
-     * Compute the mapping related data for the given @p mapping, @p cell,
-     * @p face_no and @p quadrature that is required by the FEPointEvaluation
-     * class.
+     * Store the requested mapping data.
      */
     void
-    compute_mapping_data_for_face_quadrature(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const unsigned int                                          face_no,
-      const Quadrature<dim - 1> &                                 quadrature,
-      MappingData &                                               mapping_data);
+    store_mapping_data(const unsigned int unit_points_index_offset,
+                       const unsigned int n_q_points,
+                       const unsigned int n_q_points_unvectorized,
+                       const MappingData &mapping_data);
+
+    /**
+     * Compute the compressed cell index.
+     */
+    unsigned int
+    compute_compressed_cell_index(const unsigned int cell_index) const;
+
+    /**
+     * Compute the geometry index offset of the current cell/face.
+     */
+    unsigned int
+    compute_geometry_index_offset(const unsigned int cell_index,
+                                  const unsigned int face_number) const;
+
+    /**
+     * Enum class for reinitialized states.
+     */
+    enum class State
+    {
+      invalid,
+      single_cell,
+      cell_vector,
+      faces_on_cells_in_vector
+    };
+
+    /**
+     * Enum class that stores the currently initialized state
+     * upon the last call of reinit().
+     */
+    State state;
 
     /**
      * The reference points specified at reinit().
+     *
+     * Indexed by @p unit_points_index.
      */
-    std::vector<Point<dim>> unit_points;
+    AlignedVector<Point<dim, VectorizedArrayType>> unit_points;
 
     /**
-     * Offset to point to the first unit point of a cell/face
+     * Offset to point to the first unit point of a cell/face.
      */
-    std::vector<unsigned int> unit_points_index;
+    AlignedVector<unsigned int> unit_points_index;
 
     /**
      * A pointer to the internal data of the underlying mapping.
      */
-    std::unique_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
+    std::shared_ptr<typename Mapping<dim, spacedim>::InternalDataBase>
       internal_mapping_data;
 
     /**
@@ -274,10 +488,54 @@ namespace NonMatching
     UpdateFlags update_flags_mapping;
 
     /**
-     * The internal data container for mapping information. The implementation
-     * is subject to future changes.
+     * Stores the index offset into the arrays @p JxW_values, @p jacobians,
+     * @p inverse_jacobians and @p normal_vectors.
      */
-    std::vector<MappingData> mapping_data;
+    AlignedVector<unsigned int> data_index_offsets;
+
+    /**
+     * The storage of the Jacobian determinant times the quadrature weight on
+     * quadrature points.
+     *
+     * Indexed by @p data_index_offsets.
+     */
+    AlignedVector<Number> JxW_values;
+
+    /**
+     * Stores the normal vectors.
+     *
+     * Indexed by @p data_index_offsets.
+     */
+    AlignedVector<Tensor<1, spacedim, Number>> normal_vectors;
+
+    /**
+     * 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.
+     */
+    AlignedVector<DerivativeForm<1, dim, spacedim, Number>> jacobians;
+
+    /**
+     * The storage of covariant transformation on quadrature points, i.e.,
+     * the inverse Jacobians of the transformation from the
+     * unit to the real cell.
+     *
+     * Indexed by @p data_index_offsets.
+     */
+    AlignedVector<DerivativeForm<1, spacedim, dim, Number>> inverse_jacobians;
+
+    /**
+     * The mapped real points.
+     *
+     * Indexed by @p data_index_offsets.
+     */
+    AlignedVector<Point<spacedim, Number>> real_points;
+
+    /**
+     * Number of unvectorized unit points per geometric entity (cell/face).
+     */
+    std::vector<unsigned int> n_q_points_unvectorized;
 
     /**
      * Offset to point to the first element of a cell in internal data
@@ -306,15 +564,16 @@ namespace NonMatching
   // ----------------------- template functions ----------------------
 
 
-  template <int dim, int spacedim>
-  MappingInfo<dim, spacedim>::MappingInfo(const Mapping<dim> &mapping,
-                                          const UpdateFlags   update_flags)
+  template <int dim, int spacedim, typename Number>
+  MappingInfo<dim, spacedim, Number>::MappingInfo(
+    const Mapping<dim> &mapping,
+    const UpdateFlags   update_flags)
     : mapping(&mapping)
     , update_flags(update_flags)
+    , update_flags_mapping(update_default)
   {
-    update_flags_mapping = update_default;
     // translate update flags
-    if (update_flags & update_jacobians || update_flags & update_JxW_values)
+    if (update_flags & update_jacobians)
       update_flags_mapping |= update_jacobians;
     if (update_flags & update_JxW_values)
       update_flags_mapping |= update_JxW_values;
@@ -327,6 +586,10 @@ namespace NonMatching
     // always save quadrature points for now
     update_flags_mapping |= update_quadrature_points;
 
+    update_flags_mapping =
+      internal::ComputeMappingDataHelper<dim, spacedim>::required_update_flags(
+        this->mapping, update_flags_mapping);
+
     // construct internal_mapping_data for MappingQ to be able to reuse it in
     // reinit() calls to avoid memory allocations
     if (const MappingQ<dim, spacedim> *mapping_q =
@@ -340,9 +603,9 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   void
-  MappingInfo<dim, spacedim>::reinit(
+  MappingInfo<dim, spacedim, Number>::reinit(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const std::vector<Point<dim>> &                             unit_points_in)
   {
@@ -351,9 +614,9 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   void
-  MappingInfo<dim, spacedim>::reinit(
+  MappingInfo<dim, spacedim, Number>::reinit(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const ArrayView<const Point<dim>> &                         unit_points_in)
   {
@@ -364,21 +627,48 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   void
-  MappingInfo<dim, spacedim>::reinit(
+  MappingInfo<dim, spacedim, Number>::reinit(
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const Quadrature<dim> &                                     quadrature)
   {
-    unit_points = quadrature.get_points();
-
-    mapping_data.resize(1);
-    CellSimilarity::Similarity cell_similarity =
-      CellSimilarity::Similarity::none;
-    compute_mapping_data_for_quadrature(cell,
-                                        cell_similarity,
-                                        quadrature,
-                                        mapping_data[0]);
+    n_q_points_unvectorized.resize(1);
+    n_q_points_unvectorized[0] = quadrature.size();
+
+    const unsigned int n_q_points =
+      compute_n_q_points<VectorizedArrayType>(n_q_points_unvectorized[0]);
+
+    const unsigned int n_q_points_data =
+      compute_n_q_points<Number>(n_q_points_unvectorized[0]);
+
+    // resize data vectors
+    resize_unit_points(n_q_points);
+    resize_data_fields(n_q_points_data);
+
+    // store unit points
+    store_unit_points(0,
+                      n_q_points,
+                      n_q_points_unvectorized[0],
+                      quadrature.get_points());
+
+    // compute mapping data
+    MappingData                mapping_data;
+    CellSimilarity::Similarity cell_similarity = CellSimilarity::none;
+    internal::ComputeMappingDataHelper<dim, spacedim>::
+      compute_mapping_data_for_quadrature(mapping,
+                                          update_flags_mapping,
+                                          cell,
+                                          cell_similarity,
+                                          quadrature,
+                                          internal_mapping_data,
+                                          mapping_data);
+
+    // store mapping data
+    store_mapping_data(0,
+                       n_q_points_data,
+                       n_q_points_unvectorized[0],
+                       mapping_data);
 
     state = State::single_cell;
     is_reinitialized();
@@ -386,10 +676,10 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   template <typename ContainerType>
   void
-  MappingInfo<dim, spacedim>::reinit_cells(
+  MappingInfo<dim, spacedim, Number>::reinit_cells(
     const ContainerType &                       cell_iterator_range,
     const std::vector<std::vector<Point<dim>>> &unit_points_vector,
     const unsigned int                          n_unfiltered_cells)
@@ -409,10 +699,10 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   template <typename ContainerType>
   void
-  MappingInfo<dim, spacedim>::reinit_cells(
+  MappingInfo<dim, spacedim, Number>::reinit_cells(
     const ContainerType &               cell_iterator_range,
     const std::vector<Quadrature<dim>> &quadrature_vector,
     const unsigned int                  n_unfiltered_cells)
@@ -425,37 +715,70 @@ namespace NonMatching
                     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)
-      unit_points_index.push_back(unit_points_index.back() + quadrature.size());
+      {
+        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();
 
-    unit_points.resize(n_unit_points);
-    mapping_data.resize(n_cells);
+    // 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;
     CellSimilarity::Similarity cell_similarity =
       CellSimilarity::Similarity::none;
     unsigned int cell_index = 0;
     for (const auto &cell : cell_iterator_range)
       {
-        auto it = unit_points.begin() + unit_points_index[cell_index];
-        for (const auto &unit_point :
-             quadrature_vector[cell_index].get_points())
-          {
-            *it = unit_point;
-            ++it;
-          }
-
-        compute_mapping_data_for_quadrature(cell,
-                                            cell_similarity,
-                                            quadrature_vector[cell_index],
-                                            mapping_data[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
+        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);
+
+        // 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);
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -470,10 +793,10 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   template <typename Iterator>
   void
-  MappingInfo<dim, spacedim>::reinit_surface(
+  MappingInfo<dim, spacedim, Number>::reinit_surface(
     const IteratorRange<Iterator> &                    cell_iterator_range,
     const std::vector<ImmersedSurfaceQuadrature<dim>> &quadrature_vector,
     const unsigned int                                 n_unfiltered_cells)
@@ -489,35 +812,70 @@ namespace NonMatching
                     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)
-      unit_points_index.push_back(unit_points_index.back() +
-                                  quadrature.get_points().size());
+      {
+        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();
 
-    unit_points.resize(n_unit_points);
-    mapping_data.resize(n_cells);
+    // 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];
 
-        auto it = unit_points.begin() + unit_points_index[cell_index];
-        for (const auto &unit_point : quadrature.get_points())
-          {
-            *it = unit_point;
-            ++it;
-          }
-
-        compute_mapping_data_for_immersed_surface_quadrature(
-          cell, quadrature, mapping_data[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
+        internal::ComputeMappingDataHelper<dim, spacedim>::
+          compute_mapping_data_for_immersed_surface_quadrature(
+            mapping,
+            update_flags_mapping,
+            cell,
+            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);
 
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -532,10 +890,10 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
+  template <int dim, int spacedim, typename Number>
   template <typename Iterator>
   void
-  MappingInfo<dim, spacedim>::reinit_faces(
+  MappingInfo<dim, spacedim, Number>::reinit_faces(
     const IteratorRange<Iterator> &                      cell_iterator_range,
     const std::vector<std::vector<Quadrature<dim - 1>>> &quadrature_vector,
     const unsigned int                                   n_unfiltered_cells)
@@ -559,10 +917,14 @@ namespace NonMatching
         ++cell_index;
       }
 
+    n_q_points_unvectorized.reserve(n_faces);
+
     // fill unit points index offset vector
     unit_points_index.resize(n_faces + 1);
+    data_index_offsets.resize(n_faces + 1);
     cell_index                 = 0;
     unsigned int n_unit_points = 0;
+    unsigned int n_data_points = 0;
     for (const auto &cell : cell_iterator_range)
       {
         for (const auto &f : cell->face_indices())
@@ -570,14 +932,26 @@ namespace NonMatching
             const unsigned int current_face_index =
               cell_index_offset[cell_index] + f;
 
-            unit_points_index[current_face_index] = n_unit_points;
-            n_unit_points +=
-              quadrature_vector[cell_index][f].get_points().size();
+            unit_points_index[current_face_index]  = n_unit_points;
+            data_index_offsets[current_face_index] = n_data_points;
+
+            const unsigned int n_points =
+              quadrature_vector[cell_index][f].size();
+            n_q_points_unvectorized.push_back(n_points);
+
+            const unsigned int n_q_points =
+              compute_n_q_points<VectorizedArrayType>(n_points);
+            n_unit_points += n_q_points;
+
+            const unsigned int n_q_points_data =
+              compute_n_q_points<Number>(n_points);
+            n_data_points += n_q_points_data;
           }
 
         ++cell_index;
       }
-    unit_points_index[n_faces] = n_unit_points;
+    unit_points_index[n_faces]  = n_unit_points;
+    data_index_offsets[n_faces] = n_data_points;
 
     // compress indices
     if (do_cell_index_compression)
@@ -585,8 +959,11 @@ namespace NonMatching
                                                  numbers::invalid_unsigned_int);
 
     // fill unit points and mapping data for every face of all cells
-    unit_points.resize(n_unit_points);
-    mapping_data.resize(n_faces);
+    // resize data vectors
+    resize_unit_points(n_unit_points);
+    resize_data_fields(n_data_points);
+
+    MappingData mapping_data;
     cell_index = 0;
     QProjector<dim> q_projector;
     for (const auto &cell : cell_iterator_range)
@@ -611,16 +988,30 @@ namespace NonMatching
             const unsigned int current_face_index =
               cell_index_offset[cell_index] + f;
 
-            auto it =
-              unit_points.begin() + unit_points_index[current_face_index];
-            for (const auto &unit_point : unit_points_on_cell)
-              {
-                *it = unit_point;
-                ++it;
-              }
-
-            compute_mapping_data_for_face_quadrature(
-              cell, f, quadrature_on_face, mapping_data[current_face_index]);
+            // store unit points
+            const unsigned int n_q_points =
+              compute_n_q_points<VectorizedArrayType>(
+                n_q_points_unvectorized[current_face_index]);
+            store_unit_points(unit_points_index[current_face_index],
+                              n_q_points,
+                              n_q_points_unvectorized[current_face_index],
+                              quadrature_on_cell.get_points());
+
+            internal::ComputeMappingDataHelper<dim, spacedim>::
+              compute_mapping_data_for_face_quadrature(mapping,
+                                                       update_flags_mapping,
+                                                       cell,
+                                                       f,
+                                                       quadrature_on_face,
+                                                       internal_mapping_data,
+                                                       mapping_data);
+
+            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);
           }
         if (do_cell_index_compression)
           cell_index_to_compressed_cell_index[cell->active_cell_index()] =
@@ -635,9 +1026,9 @@ namespace NonMatching
 
 
 
-  template <int dim, int spacedim>
-  inline const ArrayView<const Point<dim>>
-  MappingInfo<dim, spacedim>::get_unit_points(
+  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
   {
@@ -647,287 +1038,327 @@ namespace NonMatching
         Assert(state == State::single_cell,
                ExcMessage(
                  "This mapping info is not reinitialized for a single cell!"));
-        return unit_points;
+        return n_q_points_unvectorized[0];
+      }
+    else
+      {
+        return n_q_points_unvectorized[compute_geometry_index_offset(
+          cell_index, face_number)];
       }
-    else if (face_number == numbers::invalid_unsigned_int)
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  template <typename NumberType>
+  unsigned int
+  MappingInfo<dim, spacedim, Number>::compute_n_q_points(
+    const unsigned int n_q_points_unvectorized)
+  {
+    const unsigned int n_lanes =
+      dealii::internal::VectorizedArrayTrait<NumberType>::width;
+    const unsigned int n_filled_lanes_last_batch =
+      n_q_points_unvectorized % n_lanes;
+    unsigned int n_q_points = n_q_points_unvectorized / n_lanes;
+    if (n_filled_lanes_last_batch > 0)
+      ++n_q_points;
+    return n_q_points;
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  unsigned int
+  MappingInfo<dim, spacedim, Number>::compute_geometry_index_offset(
+    const unsigned int cell_index,
+    const unsigned int face_number) const
+  {
+    const unsigned int compressed_cell_index =
+      compute_compressed_cell_index(cell_index);
+    if (face_number == numbers::invalid_unsigned_int)
       {
         Assert(state == State::cell_vector,
                ExcMessage(
                  "This mapping info is not reinitialized for a cell vector!"));
-        if (do_cell_index_compression)
-          {
-            Assert(cell_index_to_compressed_cell_index[cell_index] !=
-                     numbers::invalid_unsigned_int,
-                   ExcMessage("Mapping info object was not initialized for this"
-                              " active cell index!"));
-            const auto it_begin =
-              unit_points.begin() +
-              unit_points_index
-                [cell_index_to_compressed_cell_index[cell_index]];
-            const auto it_end =
-              unit_points.begin() +
-              unit_points_index
-                [cell_index_to_compressed_cell_index[cell_index] + 1];
-            return make_array_view(it_begin, it_end);
-          }
-        else
-          {
-            const auto it_begin =
-              unit_points.begin() + unit_points_index[cell_index];
-            const auto it_end =
-              unit_points.begin() + unit_points_index[cell_index + 1];
-            return make_array_view(it_begin, it_end);
-          }
+        return compressed_cell_index;
       }
-    else if (cell_index != numbers::invalid_unsigned_int)
+    else
       {
+        Assert(cell_index != numbers::invalid_unsigned_int,
+               ExcMessage(
+                 "cell_index has to be set if face_number is specified!"));
         Assert(state == State::faces_on_cells_in_vector,
                ExcMessage("This mapping info is not reinitialized for faces"
                           " on cells in a vector!"));
-        if (do_cell_index_compression)
-          {
-            Assert(
-              cell_index_to_compressed_cell_index[cell_index] !=
-                numbers::invalid_unsigned_int,
-              ExcMessage(
-                "Mapping info object was not initialized for this active cell index"
-                " and corresponding face numbers!"));
-            const unsigned int current_face_index =
-              cell_index_offset
-                [cell_index_to_compressed_cell_index[cell_index]] +
-              face_number;
-            const auto it_begin =
-              unit_points.begin() + unit_points_index[current_face_index];
-            const auto it_end =
-              unit_points.begin() + unit_points_index[current_face_index + 1];
-            return make_array_view(it_begin, it_end);
-          }
-        else
-          {
-            const unsigned int current_face_index =
-              cell_index_offset[cell_index] + face_number;
-            const auto it_begin =
-              unit_points.begin() + unit_points_index[current_face_index];
-            const auto it_end =
-              unit_points.begin() + unit_points_index[current_face_index + 1];
-            return make_array_view(it_begin, it_end);
-          }
+        return cell_index_offset[compressed_cell_index] + face_number;
       }
-    else
-      AssertThrow(
-        false,
-        ExcMessage(
-          "cell_index has to be specified if face_number is specified!"));
   }
 
 
 
-  template <int dim, int spacedim>
-  inline const typename MappingInfo<dim, spacedim>::MappingData &
-  MappingInfo<dim, spacedim>::get_mapping_data(
-    const unsigned int cell_index,
-    const unsigned int face_number) const
+  template <int dim, int spacedim, typename Number>
+  unsigned int
+  MappingInfo<dim, spacedim, Number>::compute_compressed_cell_index(
+    const unsigned int cell_index) const
   {
-    if (cell_index == numbers::invalid_unsigned_int &&
-        face_number == numbers::invalid_unsigned_int)
+    if (do_cell_index_compression)
       {
-        Assert(state == State::single_cell,
-               ExcMessage(
-                 "This mapping info is not reinitialized for a single cell!"));
-        return mapping_data[0];
+        Assert(cell_index_to_compressed_cell_index[cell_index] !=
+                 numbers::invalid_unsigned_int,
+               ExcMessage("Mapping info object was not initialized for this"
+                          " active cell index!"));
+        return cell_index_to_compressed_cell_index[cell_index];
       }
-    else if (face_number == numbers::invalid_unsigned_int)
+    else
+      return cell_index;
+  }
+
+
+  template <int dim, int spacedim, typename Number>
+  void
+  MappingInfo<dim, spacedim, Number>::store_unit_points(
+    const unsigned int             unit_points_index_offset,
+    const unsigned int             n_q_points,
+    const unsigned int             n_q_points_unvectorized,
+    const std::vector<Point<dim>> &points)
+  {
+    const unsigned int n_lanes =
+      dealii::internal::VectorizedArrayTrait<VectorizedArrayType>::width;
+
+    for (unsigned int q = 0; q < n_q_points; ++q)
       {
-        Assert(state == State::cell_vector,
-               ExcMessage(
-                 "This mapping info is not reinitialized for a cell vector!"));
-        if (do_cell_index_compression)
-          {
-            Assert(cell_index_to_compressed_cell_index[cell_index] !=
-                     numbers::invalid_unsigned_int,
-                   ExcMessage("Mapping info object was not initialized for this"
-                              " active cell index!"));
-            return mapping_data
-              [cell_index_to_compressed_cell_index[cell_index]];
-          }
-        else
-          return mapping_data[cell_index];
+        const unsigned int offset = unit_points_index_offset + q;
+        for (unsigned int v = 0;
+             v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
+             ++v)
+          for (unsigned int d = 0; d < dim; ++d)
+            dealii::internal::VectorizedArrayTrait<VectorizedArrayType>::get(
+              unit_points[offset][d], v) = points[q * n_lanes + v][d];
       }
-    else if (cell_index != numbers::invalid_unsigned_int)
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  void
+  MappingInfo<dim, spacedim, Number>::store_mapping_data(
+    const unsigned int              unit_points_index_offset,
+    const unsigned int              n_q_points,
+    const unsigned int              n_q_points_unvectorized,
+    const MappingInfo::MappingData &mapping_data)
+  {
+    const unsigned int n_lanes =
+      dealii::internal::VectorizedArrayTrait<Number>::width;
+
+    for (unsigned int q = 0; q < n_q_points; ++q)
       {
-        Assert(state == State::faces_on_cells_in_vector,
-               ExcMessage("This mapping info is not reinitialized for faces"
-                          " on cells in a vector!"));
-        if (do_cell_index_compression)
+        const unsigned int offset = unit_points_index_offset + q;
+        for (unsigned int v = 0;
+             v < n_lanes && q * n_lanes + v < n_q_points_unvectorized;
+             ++v)
           {
-            Assert(
-              cell_index_to_compressed_cell_index[cell_index] !=
-                numbers::invalid_unsigned_int,
-              ExcMessage(
-                "Mapping info object was not initialized for this active cell index"
-                " and corresponding face numbers!"));
-            return mapping_data
-              [cell_index_offset
-                 [cell_index_to_compressed_cell_index[cell_index]] +
-               face_number];
+            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 (update_flags_mapping & UpdateFlags::update_JxW_values)
+              dealii::internal::VectorizedArrayTrait<Number>::get(
+                JxW_values[offset], v) =
+                mapping_data.JxW_values[q * n_lanes + v];
+            if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+              for (unsigned int s = 0; s < spacedim; ++s)
+                dealii::internal::VectorizedArrayTrait<Number>::get(
+                  normal_vectors[offset][s], v) =
+                  mapping_data.normal_vectors[q * n_lanes + v][s];
+            if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+              for (unsigned int s = 0; s < spacedim; ++s)
+                dealii::internal::VectorizedArrayTrait<Number>::get(
+                  real_points[offset][s], v) =
+                  mapping_data.quadrature_points[q * n_lanes + v][s];
           }
-        else
-          return mapping_data[cell_index_offset[cell_index] + face_number];
       }
-    else
-      AssertThrow(
-        false,
-        ExcMessage(
-          "cell_index has to be specified if face_number is specified!"));
   }
 
 
 
-  template <int dim, int spacedim>
-  const Mapping<dim, spacedim> &
-  MappingInfo<dim, spacedim>::get_mapping() const
+  template <int dim, int spacedim, typename Number>
+  void
+  MappingInfo<dim, spacedim, Number>::resize_unit_points(
+    const unsigned int n_unit_point_batches)
   {
-    return *mapping;
+    unit_points.resize(n_unit_point_batches);
   }
 
 
 
-  template <int dim, int spacedim>
-  UpdateFlags
-  MappingInfo<dim, spacedim>::get_update_flags() const
+  template <int dim, int spacedim, typename Number>
+  void
+  MappingInfo<dim, spacedim, Number>::resize_data_fields(
+    const unsigned int n_data_point_batches)
   {
-    return update_flags;
+    if (update_flags_mapping & UpdateFlags::update_jacobians)
+      jacobians.resize(n_data_point_batches);
+    if (update_flags_mapping & UpdateFlags::update_inverse_jacobians)
+      inverse_jacobians.resize(n_data_point_batches);
+    if (update_flags_mapping & UpdateFlags::update_JxW_values)
+      JxW_values.resize(n_data_point_batches);
+    if (update_flags_mapping & UpdateFlags::update_normal_vectors)
+      normal_vectors.resize(n_data_point_batches);
+    if (update_flags_mapping & UpdateFlags::update_quadrature_points)
+      real_points.resize(n_data_point_batches);
   }
 
 
 
-  template <int dim, int spacedim>
-  boost::signals2::connection
-  MappingInfo<dim, spacedim>::connect_is_reinitialized(
-    const std::function<void()> &set_is_reinitialized)
+  template <int dim, int spacedim, typename Number>
+  inline const Point<
+    dim,
+    typename MappingInfo<dim, spacedim, Number>::VectorizedArrayType> *
+  MappingInfo<dim, spacedim, Number>::get_unit_point(
+    const unsigned int offset) const
   {
-    return is_reinitialized.connect(set_is_reinitialized);
+    return &unit_points[offset];
   }
 
 
 
-  template <int dim, int spacedim>
-  void
-  MappingInfo<dim, spacedim>::compute_mapping_data_for_quadrature(
-    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    CellSimilarity::Similarity &                                cell_similarity,
-    const Quadrature<dim> &                                     quadrature,
-    MappingData &                                               mapping_data)
+  template <int dim, int spacedim, typename Number>
+  inline const Point<dim, Number> *
+  MappingInfo<dim, spacedim, Number>::get_real_point(
+    const unsigned int offset) const
   {
-    update_flags_mapping |=
-      mapping->requires_update_flags(update_flags_mapping);
+    return &real_points[offset];
+  }
 
-    mapping_data.initialize(quadrature.size(), update_flags_mapping);
 
-    // reuse internal_mapping_data for MappingQ to avoid memory allocations
-    if (const MappingQ<dim, spacedim> *mapping_q =
-          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+
+  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
+  {
+    if (cell_index == numbers::invalid_unsigned_int &&
+        face_number == numbers::invalid_unsigned_int)
       {
-        (void)mapping_q;
-        auto &data =
-          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
-            *internal_mapping_data);
-        data.initialize(update_flags_mapping, quadrature, quadrature.size());
+        Assert(state == State::single_cell,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a single cell!"));
+        return 0;
       }
     else
       {
-        internal_mapping_data =
-          mapping->get_data(update_flags_mapping, quadrature);
+        const unsigned int offset =
+          compute_geometry_index_offset(cell_index, face_number);
+        return unit_points_index[offset];
       }
-
-    cell_similarity = mapping->fill_fe_values(
-      cell, cell_similarity, quadrature, *internal_mapping_data, mapping_data);
   }
 
 
 
-  template <int dim, int spacedim>
-  void
-  MappingInfo<dim, spacedim>::
-    compute_mapping_data_for_immersed_surface_quadrature(
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      const ImmersedSurfaceQuadrature<dim> &                      quadrature,
-      MappingData &                                               mapping_data)
+  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
   {
-    update_flags_mapping |=
-      mapping->requires_update_flags(update_flags_mapping);
-
-    mapping_data.initialize(quadrature.size(), update_flags_mapping);
-
-    // reuse internal_mapping_data for MappingQ to avoid memory allocations
-    if (const MappingQ<dim, spacedim> *mapping_q =
-          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
+    if (cell_index == numbers::invalid_unsigned_int &&
+        face_number == numbers::invalid_unsigned_int)
       {
-        (void)mapping_q;
-        auto &data =
-          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
-            *internal_mapping_data);
-        data.initialize(update_flags_mapping, quadrature, quadrature.size());
+        Assert(state == State::single_cell,
+               ExcMessage(
+                 "This mapping info is not reinitialized for a single cell!"));
+        return 0;
       }
     else
       {
-        internal_mapping_data =
-          mapping->get_data(update_flags_mapping, quadrature);
+        const unsigned int offset =
+          compute_geometry_index_offset(cell_index, face_number);
+        return data_index_offsets[offset];
       }
+  }
+
 
-    mapping->fill_fe_immersed_surface_values(cell,
-                                             quadrature,
-                                             *internal_mapping_data,
-                                             mapping_data);
+
+  template <int dim, int spacedim, typename Number>
+  inline const DerivativeForm<1, dim, spacedim, Number> *
+  MappingInfo<dim, spacedim, Number>::get_jacobian(
+    const unsigned int offset) const
+  {
+    return &jacobians[offset];
   }
 
 
 
-  template <int dim, int spacedim>
-  void
-  MappingInfo<dim, spacedim>::compute_mapping_data_for_face_quadrature(
-    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    const unsigned int                                          face_no,
-    const Quadrature<dim - 1> &                                 quadrature,
-    MappingData &                                               mapping_data)
+  template <int dim, int spacedim, typename Number>
+  inline const DerivativeForm<1, spacedim, dim, Number> *
+  MappingInfo<dim, spacedim, Number>::get_inverse_jacobian(
+    const unsigned int offset) const
   {
-    update_flags_mapping |=
-      mapping->requires_update_flags(update_flags_mapping);
+    return &inverse_jacobians[offset];
+  }
 
-    mapping_data.initialize(quadrature.size(), update_flags_mapping);
 
-    // reuse internal_mapping_data for MappingQ to avoid memory allocations
-    if (const MappingQ<dim, spacedim> *mapping_q =
-          dynamic_cast<const MappingQ<dim, spacedim> *>(&(*mapping)))
-      {
-        auto &data =
-          dynamic_cast<typename MappingQ<dim, spacedim>::InternalData &>(
-            *internal_mapping_data);
-        data.initialize_face(update_flags_mapping,
-                             QProjector<dim>::project_to_oriented_face(
-                               ReferenceCells::get_hypercube<dim>(),
-                               quadrature,
-                               face_no,
-                               cell->face_orientation(face_no),
-                               cell->face_flip(face_no),
-                               cell->face_rotation(face_no)),
-                             quadrature.size());
-
-        mapping_q->fill_mapping_data_for_face_quadrature(
-          cell, face_no, quadrature, *internal_mapping_data, mapping_data);
-      }
-    else
-      {
-        auto internal_mapping_data =
-          mapping->get_face_data(update_flags_mapping,
-                                 hp::QCollection<dim - 1>(quadrature));
-
-        mapping->fill_fe_face_values(cell,
-                                     face_no,
-                                     hp::QCollection<dim - 1>(quadrature),
-                                     *internal_mapping_data,
-                                     mapping_data);
-      }
+  template <int dim, int spacedim, typename Number>
+  inline const Tensor<1, spacedim, Number> *
+  MappingInfo<dim, spacedim, Number>::get_normal_vector(
+    const unsigned int offset) const
+  {
+    return &normal_vectors[offset];
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  inline const Number *
+  MappingInfo<dim, spacedim, Number>::get_JxW(const unsigned int offset) const
+  {
+    return &JxW_values[offset];
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  const Mapping<dim, spacedim> &
+  MappingInfo<dim, spacedim, Number>::get_mapping() const
+  {
+    return *mapping;
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  UpdateFlags
+  MappingInfo<dim, spacedim, Number>::get_update_flags() const
+  {
+    return update_flags;
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  UpdateFlags
+  MappingInfo<dim, spacedim, Number>::get_update_flags_mapping() const
+  {
+    return update_flags_mapping;
+  }
+
+
+
+  template <int dim, int spacedim, typename Number>
+  boost::signals2::connection
+  MappingInfo<dim, spacedim, Number>::connect_is_reinitialized(
+    const std::function<void()> &set_is_reinitialized)
+  {
+    return is_reinitialized.connect(set_is_reinitialized);
   }
 } // namespace NonMatching
 
index b1e8eb1676db3b30361c6ac2af648f62e137c52c..bd3d1dcf507a9e1b74432c4f2edc37f0f888d7b8 100644 (file)
@@ -186,8 +186,10 @@ MappingQ<dim, spacedim>::InternalData::initialize_face(
       if (this->update_each &
           (update_boundary_forms | update_normal_vectors | update_JxW_values))
         {
-          aux.resize(dim - 1,
-                     AlignedVector<Tensor<1, spacedim>>(n_original_q_points));
+          aux.resize(dim - 1);
+          aux[0].resize(n_original_q_points);
+          if (dim > 2)
+            aux[1].resize(n_original_q_points);
 
           // Compute tangentials to the unit cell.
           for (const unsigned int i : GeometryInfo<dim>::face_indices())
index 612ed9185e9872db121d658d32bbb6d9beef6361..a2b5b98dce142881e4a12c7835799a1e11c942ee 100644 (file)
 
 using namespace dealii;
 
+template <int dim>
 void
 test(const bool filtered_compression)
 {
-  constexpr unsigned int dim    = 2;
   constexpr unsigned int degree = 1;
 
   FE_Q<dim> fe_q(degree);
@@ -471,7 +471,11 @@ main(int argc, char **argv)
 
   initlog();
 
-  test(true);
+  test<2>(true);
   deallog << std::endl;
-  test(false);
+  test<2>(false);
+  deallog << std::endl;
+  test<3>(true);
+  deallog << std::endl;
+  test<3>(false);
 }
index b4737523f919846ccd4588c4df66e3db7a3c05e8..5da28926fd8315a99801af98fa42c439e61169bf 100644 (file)
@@ -6,3 +6,11 @@ DEAL::
 DEAL::check difference l2 norm cell: 0.00000
 DEAL::check difference l2 norm surface: 0.00000
 DEAL::check difference l2 norm faces: 0.00000
+DEAL::
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000
+DEAL::
+DEAL::check difference l2 norm cell: 0.00000
+DEAL::check difference l2 norm surface: 0.00000
+DEAL::check difference l2 norm faces: 0.00000

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.