]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Put implementation details of FEFaceEvaluation into separate file 8713/head
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 9 Sep 2019 08:35:55 +0000 (10:35 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 9 Sep 2019 08:35:55 +0000 (10:35 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/operators.h
include/deal.II/matrix_free/vector_access_internal.h [new file with mode: 0644]

index de9e359c8292fadc6f9072ff1ab4f5e66297f807..e77ae68e113b0a045cbeb2ebcd1474a07a8b1f3a 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/utilities.h>
 #include <deal.II/base/vectorization.h>
 
+#include <deal.II/matrix_free/dof_info.h>
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 
@@ -1730,6 +1731,1420 @@ namespace internal
         }
     }
   };
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_read(const Number2 *src_ptr, VectorizedArrayType &dst)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+      dst[v] = src_ptr[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated load function
+  template <typename Number, unsigned int width>
+  void
+  do_vectorized_read(const Number *src_ptr, VectorizedArray<Number, width> &dst)
+  {
+    dst.load(src_ptr);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_gather(const Number2 *      src_ptr,
+                       const unsigned int * indices,
+                       VectorizedArrayType &dst)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+      dst[v] = src_ptr[indices[v]];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated gather function
+  template <typename Number, unsigned int width>
+  void
+  do_vectorized_gather(const Number *                  src_ptr,
+                       const unsigned int *            indices,
+                       VectorizedArray<Number, width> &dst)
+  {
+    dst.gather(src_ptr, indices);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_add(const VectorizedArrayType src, Number2 *dst_ptr)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+      dst_ptr[v] += src[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated load function
+  template <typename Number, unsigned int width>
+  void
+  do_vectorized_add(const VectorizedArray<Number, width> src, Number *dst_ptr)
+  {
+    VectorizedArray<Number, width> tmp;
+    tmp.load(dst_ptr);
+    (tmp + src).store(dst_ptr);
+  }
+
+
+
+  // internal helper function for reading data; base version of different types
+  template <typename VectorizedArrayType, typename Number2>
+  void
+  do_vectorized_scatter_add(const VectorizedArrayType src,
+                            const unsigned int *      indices,
+                            Number2 *                 dst_ptr)
+  {
+    for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+      dst_ptr[indices[v]] += src[v];
+  }
+
+
+
+  // internal helper function for reading data; specialized version where we
+  // can use a dedicated gather function
+  template <typename Number, unsigned int width>
+  void
+  do_vectorized_scatter_add(const VectorizedArray<Number, width> src,
+                            const unsigned int *                 indices,
+                            Number *                             dst_ptr)
+  {
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
+    for (unsigned int v = 0; v < width; ++v)
+      dst_ptr[indices[v]] += src[v];
+#else
+    VectorizedArray<Number, width> tmp;
+    tmp.gather(dst_ptr, indices);
+    (tmp + src).scatter(indices, dst_ptr);
+#endif
+  }
+
+
+
+  template <int dim,
+            int fe_degree,
+            int n_q_points_1d,
+            int n_components,
+            typename Number,
+            typename VectorizedArrayType,
+            typename Number2 = Number>
+  struct FEFaceEvaluationSelector
+  {
+    static void
+    evaluate(const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+             const VectorizedArrayType *   values_array,
+             VectorizedArrayType *         values_quad,
+             VectorizedArrayType *         gradients_quad,
+             VectorizedArrayType *         scratch_data,
+             const bool                    evaluate_values,
+             const bool                    evaluate_gradients,
+             const unsigned int            face_no,
+             const unsigned int            subface_index,
+             const unsigned int            face_orientation,
+             const Table<2, unsigned int> &orientation_map)
+    {
+      constexpr unsigned int static_dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         numbers::invalid_unsigned_int;
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? static_dofs_per_face :
+                         Utilities::pow(data.fe_degree + 1, dim - 1);
+
+      // we allocate small amounts of data on the stack to signal the compiler
+      // that this temporary data is only needed for the calculations but the
+      // final results can be discarded and need not be written back to
+      // memory. For large sizes or when the dofs per face is not a compile-time
+      // constant, however, we want to go to the heap in the `scratch_data`
+      // variable to not risk a stack overflow.
+      constexpr unsigned int stack_array_size_threshold = 100;
+
+      VectorizedArrayType
+                           temp_data[static_dofs_per_face < stack_array_size_threshold ?
+                    n_components * 2 * static_dofs_per_face :
+                    1];
+      VectorizedArrayType *temp1;
+      if (static_dofs_per_face < stack_array_size_threshold)
+        temp1 = &temp_data[0];
+      else
+        temp1 = scratch_data;
+
+      FEFaceNormalEvaluationImpl<dim,
+                                 fe_degree,
+                                 n_components,
+                                 VectorizedArrayType>::
+        template interpolate<true, false>(
+          data, values_array, temp1, evaluate_gradients, face_no);
+
+      const unsigned int n_q_points_1d_actual =
+        fe_degree > -1 ? n_q_points_1d : 0;
+      if (fe_degree > -1 &&
+          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+          data.element_type <= MatrixFreeFunctions::tensor_symmetric)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          n_components,
+          VectorizedArrayType>::evaluate_in_face(data,
+                                                 temp1,
+                                                 values_quad,
+                                                 gradients_quad,
+                                                 scratch_data + 2 *
+                                                                  n_components *
+                                                                  dofs_per_face,
+                                                 evaluate_values,
+                                                 evaluate_gradients,
+                                                 subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          n_components,
+          VectorizedArrayType>::evaluate_in_face(data,
+                                                 temp1,
+                                                 values_quad,
+                                                 gradients_quad,
+                                                 scratch_data + 2 *
+                                                                  n_components *
+                                                                  dofs_per_face,
+                                                 evaluate_values,
+                                                 evaluate_gradients,
+                                                 subface_index);
+
+      if (face_orientation)
+        adjust_for_face_orientation(face_orientation,
+                                    orientation_map,
+                                    false,
+                                    evaluate_values,
+                                    evaluate_gradients,
+                                    data.n_q_points_face,
+                                    scratch_data,
+                                    values_quad,
+                                    gradients_quad);
+    }
+
+    static void
+    integrate(const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+              VectorizedArrayType *         values_array,
+              VectorizedArrayType *         values_quad,
+              VectorizedArrayType *         gradients_quad,
+              VectorizedArrayType *         scratch_data,
+              const bool                    integrate_values,
+              const bool                    integrate_gradients,
+              const unsigned int            face_no,
+              const unsigned int            subface_index,
+              const unsigned int            face_orientation,
+              const Table<2, unsigned int> &orientation_map)
+    {
+      if (face_orientation)
+        adjust_for_face_orientation(face_orientation,
+                                    orientation_map,
+                                    true,
+                                    integrate_values,
+                                    integrate_gradients,
+                                    data.n_q_points_face,
+                                    scratch_data,
+                                    values_quad,
+                                    gradients_quad);
+
+      constexpr unsigned int static_dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         numbers::invalid_unsigned_int;
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? static_dofs_per_face :
+                         Utilities::pow(data.fe_degree + 1, dim - 1);
+
+      constexpr unsigned int stack_array_size_threshold = 100;
+
+      VectorizedArrayType
+                           temp_data[static_dofs_per_face < stack_array_size_threshold ?
+                    n_components * 2 * static_dofs_per_face :
+                    1];
+      VectorizedArrayType *temp1;
+      if (static_dofs_per_face < stack_array_size_threshold)
+        temp1 = &temp_data[0];
+      else
+        temp1 = scratch_data;
+
+      const unsigned int n_q_points_1d_actual =
+        fe_degree > -1 ? n_q_points_1d : 0;
+      if (fe_degree > -1 &&
+          subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
+          data.element_type <= MatrixFreeFunctions::tensor_symmetric)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          n_components,
+          VectorizedArrayType>::integrate_in_face(data,
+                                                  temp1,
+                                                  values_quad,
+                                                  gradients_quad,
+                                                  scratch_data +
+                                                    2 * n_components *
+                                                      dofs_per_face,
+                                                  integrate_values,
+                                                  integrate_gradients,
+                                                  subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d_actual,
+          n_components,
+          VectorizedArrayType>::integrate_in_face(data,
+                                                  temp1,
+                                                  values_quad,
+                                                  gradients_quad,
+                                                  scratch_data +
+                                                    2 * n_components *
+                                                      dofs_per_face,
+                                                  integrate_values,
+                                                  integrate_gradients,
+                                                  subface_index);
+
+      FEFaceNormalEvaluationImpl<dim,
+                                 fe_degree,
+                                 n_components,
+                                 VectorizedArrayType>::
+        template interpolate<false, false>(
+          data, temp1, values_array, integrate_gradients, face_no);
+    }
+
+    static bool
+    gather_evaluate(
+      const Number2 *                                            src_ptr,
+      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+      const MatrixFreeFunctions::DoFInfo &                       dof_info,
+      VectorizedArrayType *                                      values_quad,
+      VectorizedArrayType *                                      gradients_quad,
+      VectorizedArrayType *                                      scratch_data,
+      const bool         evaluate_values,
+      const bool         evaluate_gradients,
+      const unsigned int active_fe_index,
+      const unsigned int first_selected_component,
+      const unsigned int cell,
+      const unsigned int face_no,
+      const unsigned int subface_index,
+      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
+      const unsigned int                                 face_orientation,
+      const Table<2, unsigned int> &                     orientation_map)
+    {
+      const unsigned int side = face_no % 2;
+
+      constexpr unsigned int static_dofs_per_component =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
+                         numbers::invalid_unsigned_int;
+      constexpr unsigned int static_dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         numbers::invalid_unsigned_int;
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? static_dofs_per_face :
+                         Utilities::pow(data.fe_degree + 1, dim - 1);
+
+      constexpr unsigned int stack_array_size_threshold = 100;
+
+      VectorizedArrayType
+        temp_data[static_dofs_per_face < stack_array_size_threshold ?
+                    n_components * 2 * dofs_per_face :
+                    1];
+      VectorizedArrayType *__restrict temp1;
+      if (static_dofs_per_face < stack_array_size_threshold)
+        temp1 = &temp_data[0];
+      else
+        temp1 = scratch_data;
+
+      // case 1: contiguous and interleaved indices
+      if (((evaluate_gradients == false &&
+            data.nodal_at_cell_boundaries == true) ||
+           (data.element_type ==
+              MatrixFreeFunctions::tensor_symmetric_hermite &&
+            fe_degree > 1)) &&
+          dof_info.index_storage_variants[dof_access_index][cell] ==
+            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+              interleaved_contiguous)
+        {
+          AssertDimension(
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+            VectorizedArrayType::n_array_elements);
+          const unsigned int dof_index =
+            dof_info.dof_indices_contiguous
+              [dof_access_index][cell * VectorizedArrayType::n_array_elements] +
+            dof_info.component_dof_indices_offset[active_fe_index]
+                                                 [first_selected_component] *
+              VectorizedArrayType::n_array_elements;
+
+          if (fe_degree > 1 && evaluate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 1 + side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 = index_array[2 * i];
+                  const unsigned int ind2 = index_array[2 * i + 1];
+                  AssertIndexRange(ind1, data.dofs_per_component_on_cell);
+                  AssertIndexRange(ind2, data.dofs_per_component_on_cell);
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      do_vectorized_read(
+                        src_ptr + dof_index +
+                          (ind1 + comp * static_dofs_per_component) *
+                            VectorizedArrayType::n_array_elements,
+                        temp1[i + 2 * comp * dofs_per_face]);
+                      do_vectorized_read(
+                        src_ptr + dof_index +
+                          (ind2 + comp * static_dofs_per_component) *
+                            VectorizedArrayType::n_array_elements,
+                        temp1[dofs_per_face + i + 2 * comp * dofs_per_face]);
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+                        grad_weight *
+                        (temp1[i + 2 * comp * dofs_per_face] -
+                         temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind = index_array[i];
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    do_vectorized_read(
+                      src_ptr + dof_index +
+                        (ind + comp * static_dofs_per_component) *
+                          VectorizedArrayType::n_array_elements,
+                      temp1[i + 2 * comp * dofs_per_face]);
+                }
+            }
+        }
+
+      // case 2: contiguous and interleaved indices with fixed stride
+      else if (((evaluate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   interleaved_contiguous_strided)
+        {
+          AssertDimension(
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+            VectorizedArrayType::n_array_elements);
+          const unsigned int *indices =
+            &dof_info.dof_indices_contiguous
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+          if (fe_degree > 1 && evaluate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 1 + side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 =
+                    index_array[2 * i] * VectorizedArrayType::n_array_elements;
+                  const unsigned int ind2 =
+                    index_array[2 * i + 1] *
+                    VectorizedArrayType::n_array_elements;
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      do_vectorized_gather(
+                        src_ptr + ind1 +
+                          comp * static_dofs_per_component *
+                            VectorizedArrayType::n_array_elements +
+                          dof_info.component_dof_indices_offset
+                              [active_fe_index][first_selected_component] *
+                            VectorizedArrayType::n_array_elements,
+                        indices,
+                        temp1[i + 2 * comp * dofs_per_face]);
+                      VectorizedArrayType grad;
+                      do_vectorized_gather(
+                        src_ptr + ind2 +
+                          comp * static_dofs_per_component *
+                            VectorizedArrayType::n_array_elements +
+                          dof_info.component_dof_indices_offset
+                              [active_fe_index][first_selected_component] *
+                            VectorizedArrayType::n_array_elements,
+                        indices,
+                        grad);
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+                        grad_weight *
+                        (temp1[i + 2 * comp * dofs_per_face] - grad);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind =
+                    index_array[i] * VectorizedArrayType::n_array_elements;
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    do_vectorized_gather(
+                      src_ptr + ind +
+                        comp * static_dofs_per_component *
+                          VectorizedArrayType::n_array_elements +
+                        dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component] *
+                          VectorizedArrayType::n_array_elements,
+                      indices,
+                      temp1[i + 2 * comp * dofs_per_face]);
+                }
+            }
+        }
+
+      // case 3: contiguous and interleaved indices with mixed stride
+      else if (((evaluate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   interleaved_contiguous_mixed_strides)
+        {
+          const unsigned int *strides =
+            &dof_info.dof_indices_interleave_strides
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+          unsigned int indices[VectorizedArrayType::n_array_elements];
+          for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements;
+               ++v)
+            indices[v] =
+              dof_info.dof_indices_contiguous
+                [dof_access_index]
+                [cell * VectorizedArrayType::n_array_elements + v] +
+              dof_info.component_dof_indices_offset[active_fe_index]
+                                                   [first_selected_component] *
+                strides[v];
+          const unsigned int nvec =
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+          if (fe_degree > 1 && evaluate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 1 + side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              if (nvec == VectorizedArrayType::n_array_elements)
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      unsigned int ind1[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind1[v] =
+                          indices[v] + (comp * static_dofs_per_component +
+                                        index_array[2 * i]) *
+                                         strides[v];
+                      unsigned int ind2[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind2[v] =
+                          indices[v] + (comp * static_dofs_per_component +
+                                        index_array[2 * i + 1]) *
+                                         strides[v];
+                      do_vectorized_gather(src_ptr,
+                                           ind1,
+                                           temp1[i + 2 * comp * dofs_per_face]);
+                      VectorizedArrayType grad;
+                      do_vectorized_gather(src_ptr, ind2, grad);
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+                        grad_weight *
+                        (temp1[i + 2 * comp * dofs_per_face] - grad);
+                    }
+              else
+                {
+                  for (unsigned int i = 0; i < n_components * 2 * dofs_per_face;
+                       ++i)
+                    temp1[i] = VectorizedArrayType();
+                  for (unsigned int v = 0; v < nvec; ++v)
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        {
+                          const unsigned int ind1 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[2 * i]) *
+                                           strides[v];
+                          const unsigned int ind2 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[2 * i + 1]) *
+                                           strides[v];
+                          temp1[i + 2 * comp * dofs_per_face][v] =
+                            src_ptr[ind1];
+                          const Number grad = src_ptr[ind2];
+                          temp1[i + dofs_per_face +
+                                2 * comp * dofs_per_face][v] =
+                            grad_weight[0] *
+                            (temp1[i + 2 * comp * dofs_per_face][v] - grad);
+                        }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              if (nvec == VectorizedArrayType::n_array_elements)
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      unsigned int ind[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind[v] =
+                          indices[v] +
+                          (comp * static_dofs_per_component + index_array[i]) *
+                            strides[v];
+                      do_vectorized_gather(src_ptr,
+                                           ind,
+                                           temp1[i + 2 * comp * dofs_per_face]);
+                    }
+              else
+                {
+                  for (unsigned int i = 0; i < n_components * dofs_per_face;
+                       ++i)
+                    temp1[i] = VectorizedArrayType();
+                  for (unsigned int v = 0; v < nvec; ++v)
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        {
+                          const unsigned int ind1 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[i]) *
+                                           strides[v];
+                          temp1[i + 2 * comp * dofs_per_face][v] =
+                            src_ptr[ind1];
+                        }
+                }
+            }
+        }
+
+      // case 4: contiguous indices without interleaving
+      else if (((evaluate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   contiguous &&
+               dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
+                 VectorizedArrayType::n_array_elements)
+        {
+          const unsigned int *indices =
+            &dof_info.dof_indices_contiguous
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+          if (evaluate_gradients == true &&
+              data.element_type ==
+                MatrixFreeFunctions::tensor_symmetric_hermite)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 1 + side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 = index_array[2 * i];
+                  const unsigned int ind2 = index_array[2 * i + 1];
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      do_vectorized_gather(
+                        src_ptr + ind1 + comp * static_dofs_per_component +
+                          dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component],
+                        indices,
+                        temp1[i + 2 * comp * dofs_per_face]);
+                      VectorizedArrayType grad;
+                      do_vectorized_gather(
+                        src_ptr + ind2 + comp * static_dofs_per_component +
+                          dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component],
+                        indices,
+                        grad);
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
+                        grad_weight *
+                        (temp1[i + 2 * comp * dofs_per_face] - grad);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  {
+                    const unsigned int ind = index_array[i];
+                    do_vectorized_gather(
+                      src_ptr + ind + comp * static_dofs_per_component +
+                        dof_info.component_dof_indices_offset
+                          [active_fe_index][first_selected_component],
+                      indices,
+                      temp1[i + comp * 2 * dofs_per_face]);
+                  }
+            }
+        }
+
+      // case 5: default vector access
+      else
+        {
+          return false;
+        }
+
+      if (fe_degree > -1 &&
+          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+          data.element_type <= MatrixFreeFunctions::tensor_symmetric)
+        FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          n_components,
+          VectorizedArrayType>::evaluate_in_face(data,
+                                                 temp1,
+                                                 values_quad,
+                                                 gradients_quad,
+                                                 scratch_data + 2 *
+                                                                  n_components *
+                                                                  dofs_per_face,
+                                                 evaluate_values,
+                                                 evaluate_gradients,
+                                                 subface_index);
+      else
+        FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          n_components,
+          VectorizedArrayType>::evaluate_in_face(data,
+                                                 temp1,
+                                                 values_quad,
+                                                 gradients_quad,
+                                                 scratch_data + 2 *
+                                                                  n_components *
+                                                                  dofs_per_face,
+                                                 evaluate_values,
+                                                 evaluate_gradients,
+                                                 subface_index);
+
+      if (face_orientation)
+        adjust_for_face_orientation(face_orientation,
+                                    orientation_map,
+                                    false,
+                                    evaluate_values,
+                                    evaluate_gradients,
+                                    data.n_q_points_face,
+                                    scratch_data,
+                                    values_quad,
+                                    gradients_quad);
+
+      return true;
+    }
+
+    static bool
+    integrate_scatter(
+      Number2 *                                                  dst_ptr,
+      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+      const MatrixFreeFunctions::DoFInfo &                       dof_info,
+      VectorizedArrayType *                                      values_array,
+      VectorizedArrayType *                                      values_quad,
+      VectorizedArrayType *                                      gradients_quad,
+      VectorizedArrayType *                                      scratch_data,
+      const bool         integrate_values,
+      const bool         integrate_gradients,
+      const unsigned int active_fe_index,
+      const unsigned int first_selected_component,
+      const unsigned int cell,
+      const unsigned int face_no,
+      const unsigned int subface_index,
+      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
+      const unsigned int                                 face_orientation,
+      const Table<2, unsigned int> &                     orientation_map)
+    {
+      if (face_orientation)
+        adjust_for_face_orientation(face_orientation,
+                                    orientation_map,
+                                    true,
+                                    integrate_values,
+                                    integrate_gradients,
+                                    data.n_q_points_face,
+                                    scratch_data,
+                                    values_quad,
+                                    gradients_quad);
+
+      const unsigned int     side = face_no % 2;
+      constexpr unsigned int static_dofs_per_component =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
+                         numbers::invalid_unsigned_int;
+      const unsigned int dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         Utilities::pow(data.fe_degree + 1, dim - 1);
+
+      constexpr unsigned int stack_array_size_threshold = 100;
+
+      VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ?
+                                      n_components * 2 * dofs_per_face :
+                                      1];
+      VectorizedArrayType *__restrict temp1;
+      if (dofs_per_face < stack_array_size_threshold)
+        temp1 = &temp_data[0];
+      else
+        temp1 = scratch_data;
+
+      if (fe_degree > -1 &&
+          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+          data.element_type <= internal::MatrixFreeFunctions::tensor_symmetric)
+        internal::FEFaceEvaluationImpl<
+          true,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          n_components,
+          VectorizedArrayType>::integrate_in_face(data,
+                                                  temp1,
+                                                  values_quad,
+                                                  gradients_quad,
+                                                  scratch_data +
+                                                    2 * n_components *
+                                                      dofs_per_face,
+                                                  integrate_values,
+                                                  integrate_gradients,
+                                                  subface_index);
+      else
+        internal::FEFaceEvaluationImpl<
+          false,
+          dim,
+          fe_degree,
+          n_q_points_1d,
+          n_components,
+          VectorizedArrayType>::integrate_in_face(data,
+                                                  temp1,
+                                                  values_quad,
+                                                  gradients_quad,
+                                                  scratch_data +
+                                                    2 * n_components *
+                                                      dofs_per_face,
+                                                  integrate_values,
+                                                  integrate_gradients,
+                                                  subface_index);
+
+      // case 1: contiguous and interleaved indices
+      if (((integrate_gradients == false &&
+            data.nodal_at_cell_boundaries == true) ||
+           (data.element_type ==
+              internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+            fe_degree > 1)) &&
+          dof_info.index_storage_variants[dof_access_index][cell] ==
+            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+              interleaved_contiguous)
+        {
+          AssertDimension(
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+            VectorizedArrayType::n_array_elements);
+          const unsigned int dof_index =
+            dof_info.dof_indices_contiguous
+              [dof_access_index][cell * VectorizedArrayType::n_array_elements] +
+            dof_info.component_dof_indices_offset[active_fe_index]
+                                                 [first_selected_component] *
+              VectorizedArrayType::n_array_elements;
+
+          if (fe_degree > 1 && integrate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 2 - side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 = index_array[2 * i];
+                  const unsigned int ind2 = index_array[2 * i + 1];
+                  AssertIndexRange(ind1, data.dofs_per_component_on_cell);
+                  AssertIndexRange(ind2, data.dofs_per_component_on_cell);
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      VectorizedArrayType val =
+                        temp1[i + 2 * comp * dofs_per_face] -
+                        grad_weight *
+                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      VectorizedArrayType grad =
+                        grad_weight *
+                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      do_vectorized_add(
+                        val,
+                        dst_ptr + dof_index +
+                          (ind1 + comp * static_dofs_per_component) *
+                            VectorizedArrayType::n_array_elements);
+                      do_vectorized_add(
+                        grad,
+                        dst_ptr + dof_index +
+                          (ind2 + comp * static_dofs_per_component) *
+                            VectorizedArrayType::n_array_elements);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind = index_array[i];
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    do_vectorized_add(
+                      temp1[i + 2 * comp * dofs_per_face],
+                      dst_ptr + dof_index +
+                        (ind + comp * static_dofs_per_component) *
+                          VectorizedArrayType::n_array_elements);
+                }
+            }
+        }
+
+      // case 2: contiguous and interleaved indices with fixed stride
+      else if (((integrate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   interleaved_contiguous_strided)
+        {
+          AssertDimension(
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+            VectorizedArrayType::n_array_elements);
+          const unsigned int *indices =
+            &dof_info.dof_indices_contiguous
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+          if (fe_degree > 1 && integrate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 2 - side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 =
+                    index_array[2 * i] * VectorizedArrayType::n_array_elements;
+                  const unsigned int ind2 =
+                    index_array[2 * i + 1] *
+                    VectorizedArrayType::n_array_elements;
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      VectorizedArrayType val =
+                        temp1[i + 2 * comp * dofs_per_face] -
+                        grad_weight *
+                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      VectorizedArrayType grad =
+                        grad_weight *
+                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      do_vectorized_scatter_add(
+                        val,
+                        indices,
+                        dst_ptr + ind1 +
+                          comp * static_dofs_per_component *
+                            VectorizedArrayType::n_array_elements +
+                          dof_info.component_dof_indices_offset
+                              [active_fe_index][first_selected_component] *
+                            VectorizedArrayType::n_array_elements);
+                      do_vectorized_scatter_add(
+                        grad,
+                        indices,
+                        dst_ptr + ind2 +
+                          comp * static_dofs_per_component *
+                            VectorizedArrayType::n_array_elements +
+                          dof_info.component_dof_indices_offset
+                              [active_fe_index][first_selected_component] *
+                            VectorizedArrayType::n_array_elements);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind =
+                    index_array[i] * VectorizedArrayType::n_array_elements;
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    do_vectorized_scatter_add(
+                      temp1[i + 2 * comp * dofs_per_face],
+                      indices,
+                      dst_ptr + ind +
+                        comp * static_dofs_per_component *
+                          VectorizedArrayType::n_array_elements +
+                        dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component] *
+                          VectorizedArrayType::n_array_elements);
+                }
+            }
+        }
+
+      // case 3: contiguous and interleaved indices with mixed stride
+      else if (((integrate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   interleaved_contiguous_mixed_strides)
+        {
+          const unsigned int *strides =
+            &dof_info.dof_indices_interleave_strides
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+          unsigned int indices[VectorizedArrayType::n_array_elements];
+          for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements;
+               ++v)
+            indices[v] =
+              dof_info.dof_indices_contiguous
+                [dof_access_index]
+                [cell * VectorizedArrayType::n_array_elements + v] +
+              dof_info.component_dof_indices_offset[active_fe_index]
+                                                   [first_selected_component] *
+                strides[v];
+          const unsigned int nvec =
+            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+
+          if (fe_degree > 1 && integrate_gradients == true)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 2 - side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              if (nvec == VectorizedArrayType::n_array_elements)
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      unsigned int ind1[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind1[v] =
+                          indices[v] + (comp * static_dofs_per_component +
+                                        index_array[2 * i]) *
+                                         strides[v];
+                      unsigned int ind2[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind2[v] =
+                          indices[v] + (comp * static_dofs_per_component +
+                                        index_array[2 * i + 1]) *
+                                         strides[v];
+                      VectorizedArrayType val =
+                        temp1[i + 2 * comp * dofs_per_face] -
+                        grad_weight *
+                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      VectorizedArrayType grad =
+                        grad_weight *
+                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      do_vectorized_scatter_add(val, ind1, dst_ptr);
+                      do_vectorized_scatter_add(grad, ind2, dst_ptr);
+                    }
+              else
+                {
+                  for (unsigned int v = 0; v < nvec; ++v)
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        {
+                          const unsigned int ind1 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[2 * i]) *
+                                           strides[v];
+                          const unsigned int ind2 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[2 * i + 1]) *
+                                           strides[v];
+                          Number val =
+                            temp1[i + 2 * comp * dofs_per_face][v] -
+                            grad_weight[0] * temp1[i + dofs_per_face +
+                                                   2 * comp * dofs_per_face][v];
+                          Number grad =
+                            grad_weight[0] * temp1[i + dofs_per_face +
+                                                   2 * comp * dofs_per_face][v];
+                          dst_ptr[ind1] += val;
+                          dst_ptr[ind2] += grad;
+                        }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              if (nvec == VectorizedArrayType::n_array_elements)
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
+                    {
+                      unsigned int ind[VectorizedArrayType::n_array_elements];
+                      DEAL_II_OPENMP_SIMD_PRAGMA
+                      for (unsigned int v = 0;
+                           v < VectorizedArrayType::n_array_elements;
+                           ++v)
+                        ind[v] =
+                          indices[v] +
+                          (comp * static_dofs_per_component + index_array[i]) *
+                            strides[v];
+                      do_vectorized_scatter_add(
+                        temp1[i + 2 * comp * dofs_per_face], ind, dst_ptr);
+                    }
+              else
+                {
+                  for (unsigned int v = 0; v < nvec; ++v)
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        {
+                          const unsigned int ind1 =
+                            indices[v] + (comp * static_dofs_per_component +
+                                          index_array[i]) *
+                                           strides[v];
+                          dst_ptr[ind1] +=
+                            temp1[i + 2 * comp * dofs_per_face][v];
+                        }
+                }
+            }
+        }
+
+      // case 4: contiguous indices without interleaving
+      else if (((integrate_gradients == false &&
+                 data.nodal_at_cell_boundaries == true) ||
+                (data.element_type ==
+                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                 fe_degree > 1)) &&
+               dof_info.index_storage_variants[dof_access_index][cell] ==
+                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                   contiguous &&
+               dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
+                 VectorizedArrayType::n_array_elements)
+        {
+          const unsigned int *indices =
+            &dof_info.dof_indices_contiguous
+               [dof_access_index][cell * VectorizedArrayType::n_array_elements];
+
+          if (integrate_gradients == true &&
+              data.element_type ==
+                internal::MatrixFreeFunctions::tensor_symmetric_hermite)
+            {
+              // we know that the gradient weights for the Hermite case on the
+              // right (side==1) are the negative from the value at the left
+              // (side==0), so we only read out one of them.
+              const VectorizedArrayType grad_weight =
+                data.shape_data_on_face[0][fe_degree + 2 - side];
+              AssertDimension(data.face_to_cell_index_hermite.size(1),
+                              2 * dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_hermite(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind1 = index_array[2 * i];
+                  const unsigned int ind2 = index_array[2 * i + 1];
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    {
+                      VectorizedArrayType val =
+                        temp1[i + 2 * comp * dofs_per_face] -
+                        grad_weight *
+                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      VectorizedArrayType grad =
+                        grad_weight *
+                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                      do_vectorized_scatter_add(
+                        val,
+                        indices,
+                        dst_ptr + comp * static_dofs_per_component + ind1 +
+                          dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component]);
+                      do_vectorized_scatter_add(
+                        grad,
+                        indices,
+                        dst_ptr + comp * static_dofs_per_component + ind2 +
+                          dof_info.component_dof_indices_offset
+                            [active_fe_index][first_selected_component]);
+                    }
+                }
+            }
+          else
+            {
+              AssertDimension(data.face_to_cell_index_nodal.size(1),
+                              dofs_per_face);
+              const unsigned int *index_array =
+                &data.face_to_cell_index_nodal(face_no, 0);
+              for (unsigned int i = 0; i < dofs_per_face; ++i)
+                {
+                  const unsigned int ind = index_array[i];
+                  for (unsigned int comp = 0; comp < n_components; ++comp)
+                    do_vectorized_scatter_add(
+                      temp1[i + 2 * comp * dofs_per_face],
+                      indices,
+                      dst_ptr + comp * static_dofs_per_component + ind +
+                        dof_info.component_dof_indices_offset
+                          [active_fe_index][first_selected_component]);
+                }
+            }
+        }
+
+      // case 5: default vector access, must be handled separately, just do
+      // the face-normal interpolation
+      else
+        {
+          FEFaceNormalEvaluationImpl<dim,
+                                     fe_degree,
+                                     n_components,
+                                     VectorizedArrayType>::
+            template interpolate<false, false>(
+              data, temp1, values_array, integrate_gradients, face_no);
+          return false;
+        }
+
+      return true;
+    }
+
+    static void
+    adjust_for_face_orientation(const unsigned int            face_orientation,
+                                const Table<2, unsigned int> &orientation_map,
+                                const bool                    integrate,
+                                const bool                    values,
+                                const bool                    gradients,
+                                const unsigned int            n_q_points,
+                                VectorizedArrayType *         tmp_values,
+                                VectorizedArrayType *         values_quad,
+                                VectorizedArrayType *         gradients_quad)
+    {
+      Assert(face_orientation, ExcInternalError());
+      const unsigned int *orientation = &orientation_map[face_orientation][0];
+      for (unsigned int c = 0; c < n_components; ++c)
+        {
+          if (values == true)
+            {
+              if (integrate)
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[q] = values_quad[c * n_q_points + orientation[q]];
+              else
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  tmp_values[orientation[q]] = values_quad[c * n_q_points + q];
+              for (unsigned int q = 0; q < n_q_points; ++q)
+                values_quad[c * n_q_points + q] = tmp_values[q];
+            }
+          if (gradients == true)
+            for (unsigned int d = 0; d < dim; ++d)
+              {
+                if (integrate)
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[q] = gradients_quad[(c * dim + d) * n_q_points +
+                                                   orientation[q]];
+                else
+                  for (unsigned int q = 0; q < n_q_points; ++q)
+                    tmp_values[orientation[q]] =
+                      gradients_quad[(c * dim + d) * n_q_points + q];
+                for (unsigned int q = 0; q < n_q_points; ++q)
+                  gradients_quad[(c * dim + d) * n_q_points + q] =
+                    tmp_values[q];
+              }
+        }
+    }
+  };
+
+
+
+  template <int dim, int fe_degree, int n_components, typename Number>
+  struct CellwiseInverseMassMatrixImpl
+  {
+    static void
+    apply(const AlignedVector<Number> &inverse_shape,
+          const AlignedVector<Number> &inverse_coefficients,
+          const unsigned int           n_desired_components,
+          const Number *               in_array,
+          Number *                     out_array)
+    {
+      constexpr unsigned int dofs_per_component =
+        Utilities::pow(fe_degree + 1, dim);
+      Assert(inverse_coefficients.size() > 0 &&
+               inverse_coefficients.size() % dofs_per_component == 0,
+             ExcMessage(
+               "Expected diagonal to be a multiple of scalar dof per cells"));
+      if (inverse_coefficients.size() != dofs_per_component)
+        AssertDimension(n_desired_components * dofs_per_component,
+                        inverse_coefficients.size());
+
+      Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
+
+      internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
+                                       dim,
+                                       fe_degree + 1,
+                                       fe_degree + 1,
+                                       Number>
+        evaluator(AlignedVector<Number>(),
+                  AlignedVector<Number>(),
+                  inverse_shape);
+
+      const unsigned int shift_coefficient =
+        inverse_coefficients.size() > dofs_per_component ? dofs_per_component :
+                                                           0;
+      const Number *inv_coefficient = inverse_coefficients.data();
+      for (unsigned int d = 0; d < n_desired_components; ++d)
+        {
+          const Number *in  = in_array + d * dofs_per_component;
+          Number *      out = out_array + d * dofs_per_component;
+          // Need to select 'apply' method with hessian slot because values
+          // assume symmetries that do not exist in the inverse shapes
+          evaluator.template hessians<0, false, false>(in, out);
+          if (dim > 1)
+            {
+              evaluator.template hessians<1, false, false>(out, out);
+
+              if (dim == 3)
+                {
+                  evaluator.template hessians<2, false, false>(out, out);
+                  for (unsigned int q = 0; q < dofs_per_component; ++q)
+                    out[q] *= inv_coefficient[q];
+                  evaluator.template hessians<2, true, false>(out, out);
+                }
+              else if (dim == 2)
+                for (unsigned int q = 0; q < dofs_per_component; ++q)
+                  out[q] *= inv_coefficient[q];
+
+              evaluator.template hessians<1, true, false>(out, out);
+            }
+          else
+            {
+              for (unsigned int q = 0; q < dofs_per_component; ++q)
+                out[q] *= inv_coefficient[q];
+            }
+          evaluator.template hessians<0, true, false>(out, out);
+
+          inv_coefficient += shift_coefficient;
+        }
+    }
+
+    static void
+    transform_from_q_points_to_basis(const AlignedVector<Number> &inverse_shape,
+                                     const unsigned int n_desired_components,
+                                     const Number *     in_array,
+                                     Number *           out_array)
+    {
+      constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
+      internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
+                                       dim,
+                                       fe_degree + 1,
+                                       fe_degree + 1,
+                                       Number>
+        evaluator(AlignedVector<Number>(),
+                  AlignedVector<Number>(),
+                  inverse_shape);
+
+      for (unsigned int d = 0; d < n_desired_components; ++d)
+        {
+          const Number *in  = in_array + d * dofs_per_cell;
+          Number *      out = out_array + d * dofs_per_cell;
+
+          if (dim == 3)
+            {
+              evaluator.template hessians<2, true, false>(in, out);
+              evaluator.template hessians<1, true, false>(out, out);
+              evaluator.template hessians<0, true, false>(out, out);
+            }
+          if (dim == 2)
+            {
+              evaluator.template hessians<1, true, false>(in, out);
+              evaluator.template hessians<0, true, false>(out, out);
+            }
+          if (dim == 1)
+            evaluator.template hessians<0, true, false>(in, out);
+        }
+    }
+  };
+
 } // end of namespace internal
 
 
index 753fb8fc5797484e75979a457ea880e1fd1952e9..9716922d9d11c8a5e0268633fca66072f6a045cf 100644 (file)
@@ -36,6 +36,7 @@
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/tensor_product_kernels.h>
 #include <deal.II/matrix_free/type_traits.h>
+#include <deal.II/matrix_free/vector_access_internal.h>
 
 
 DEAL_II_NAMESPACE_OPEN
@@ -2918,17 +2919,6 @@ public:
    * points is inaccurate and this value must be used instead.
    */
   const unsigned int n_q_points;
-
-protected:
-  /**
-   * For faces not oriented in the standard way, this method applies
-   * re-indexing on quadrature points. Called at the end of evaluate() and at
-   * the beginning of integrate().
-   */
-  void
-  adjust_for_face_orientation(const bool integrate,
-                              const bool values,
-                              const bool gradients);
 };
 
 
@@ -3590,654 +3580,6 @@ FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
 
 namespace internal
 {
-  // below we use type-traits from matrix-free/type_traits.h
-
-  // access to generic const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
-  // where we should first do Partitioner::local_to_global()
-  template <typename VectorType,
-            typename std::enable_if<!has_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline typename VectorType::value_type
-  vector_access(const VectorType &vec, const unsigned int entry)
-  {
-    return vec(entry);
-  }
-
-
-
-  // access to generic non-const vectors that have operator ().
-  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
-  // where we should first do Partitioner::local_to_global()
-  template <typename VectorType,
-            typename std::enable_if<!has_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline typename VectorType::value_type &
-  vector_access(VectorType &vec, const unsigned int entry)
-  {
-    return vec(entry);
-  }
-
-
-
-  // access to distributed MPI vectors that have a local_element(uint)
-  // method to access data in local index space, which is what we use in
-  // DoFInfo and hence in read_dof_values etc.
-  template <typename VectorType,
-            typename std::enable_if<has_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline typename VectorType::value_type &
-  vector_access(VectorType &vec, const unsigned int entry)
-  {
-    return vec.local_element(entry);
-  }
-
-
-
-  // same for const access
-  template <typename VectorType,
-            typename std::enable_if<has_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline typename VectorType::value_type
-  vector_access(const VectorType &vec, const unsigned int entry)
-  {
-    return vec.local_element(entry);
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<has_add_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_add(VectorType &                           vec,
-                    const unsigned int                     entry,
-                    const typename VectorType::value_type &val)
-  {
-    vec.add_local_element(entry, val);
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<!has_add_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_add(VectorType &                           vec,
-                    const unsigned int                     entry,
-                    const typename VectorType::value_type &val)
-  {
-    vector_access(vec, entry) += val;
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<has_add_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_add_global(VectorType &                           vec,
-                           const types::global_dof_index          entry,
-                           const typename VectorType::value_type &val)
-  {
-    vec.add(entry, val);
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<!has_add_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_add_global(VectorType &                           vec,
-                           const types::global_dof_index          entry,
-                           const typename VectorType::value_type &val)
-  {
-    vec(entry) += val;
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<has_set_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_set(VectorType &                           vec,
-                    const unsigned int                     entry,
-                    const typename VectorType::value_type &val)
-  {
-    vec.set_local_element(entry, val);
-  }
-
-
-
-  template <typename VectorType,
-            typename std::enable_if<!has_set_local_element<VectorType>::value,
-                                    VectorType>::type * = nullptr>
-  inline void
-  vector_access_set(VectorType &                           vec,
-                    const unsigned int                     entry,
-                    const typename VectorType::value_type &val)
-  {
-    vector_access(vec, entry) = val;
-  }
-
-
-
-  // this is to make sure that the parallel partitioning in VectorType
-  // is really the same as stored in MatrixFree.
-  // version below is when has_partitioners_are_compatible == false
-  // FIXME: this is incorrect for PETSc/Trilinos MPI vectors
-  template <
-    typename VectorType,
-    typename std::enable_if<!has_partitioners_are_compatible<VectorType>::value,
-                            VectorType>::type * = nullptr>
-  inline void
-  check_vector_compatibility(
-    const VectorType &                            vec,
-    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
-  {
-    (void)vec;
-    (void)dof_info;
-
-    AssertDimension(vec.size(), dof_info.vector_partitioner->size());
-  }
-
-
-
-  // same as above for has_partitioners_are_compatible == true
-  template <
-    typename VectorType,
-    typename std::enable_if<has_partitioners_are_compatible<VectorType>::value,
-                            VectorType>::type * = nullptr>
-  inline void
-  check_vector_compatibility(
-    const VectorType &                            vec,
-    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
-  {
-    (void)vec;
-    (void)dof_info;
-    Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
-           ExcMessage(
-             "The parallel layout of the given vector is not "
-             "compatible with the parallel partitioning in MatrixFree. "
-             "Use MatrixFree::initialize_dof_vector to get a "
-             "compatible vector."));
-  }
-
-
-
-  // Below, three classes (VectorReader, VectorSetter,
-  // VectorDistributorLocalToGlobal) implement the same interface and can be
-  // used to to read from vector, set elements of a vector and add to elements
-  // of the vector.
-
-  // 1. A class to read data from vector
-  template <typename Number, typename VectorizedArrayType>
-  struct VectorReader
-  {
-    template <typename VectorType>
-    void
-    process_dof(const unsigned int index,
-                const VectorType & vec,
-                Number &           res) const
-    {
-      res = vector_access(vec, index);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            VectorType &         vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, true>) const
-    {
-      const Number *vec_ptr = vec.begin() + dof_index;
-      for (unsigned int i = 0; i < dofs_per_cell;
-           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
-        dof_values[i].load(vec_ptr);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            const VectorType &   vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, false>) const
-    {
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          dof_values[i][v] = vector_access(
-            vec, dof_index + v + i * VectorizedArrayType::n_array_elements);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      VectorType &         vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, true>) const
-    {
-      dealii::vectorized_load_and_transpose(dofs_per_cell,
-                                            vec.begin(),
-                                            dof_indices,
-                                            dof_values);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      const VectorType &   vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, false>) const
-    {
-      for (unsigned int d = 0; d < dofs_per_cell; ++d)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          dof_values[d][v] = vector_access(vec, dof_indices[v] + d);
-    }
-
-
-
-    // variant where VectorType::value_type is the same as Number -> can call
-    // gather
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       VectorType &         vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, true>) const
-    {
-      res.gather(vec.begin() + constant_offset, indices);
-    }
-
-
-
-    // variant where VectorType::value_type is not the same as Number -> must
-    // manually load the data
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       const VectorType &   vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, false>) const
-    {
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        res[v] = vector_access(vec, indices[v] + constant_offset);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dof_global(const types::global_dof_index index,
-                       const VectorType &            vec,
-                       Number &                      res) const
-    {
-      res = vec(index);
-    }
-
-
-
-    void
-    pre_constraints(const Number &, Number &res) const
-    {
-      res = Number();
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_constraint(const unsigned int index,
-                       const Number       weight,
-                       const VectorType & vec,
-                       Number &           res) const
-    {
-      res += weight * vector_access(vec, index);
-    }
-
-
-
-    void
-    post_constraints(const Number &sum, Number &write_pos) const
-    {
-      write_pos = sum;
-    }
-
-
-
-    void
-    process_empty(VectorizedArrayType &res) const
-    {
-      res = VectorizedArrayType();
-    }
-  };
-
-
-
-  // 2. A class to add values to the vector during
-  // FEEvaluation::distribute_local_to_global() call
-  template <typename Number, typename VectorizedArrayType>
-  struct VectorDistributorLocalToGlobal
-  {
-    template <typename VectorType>
-    void
-    process_dof(const unsigned int index, VectorType &vec, Number &res) const
-    {
-      vector_access_add(vec, index, res);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            VectorType &         vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, true>) const
-    {
-      Number *vec_ptr = vec.begin() + dof_index;
-      for (unsigned int i = 0; i < dofs_per_cell;
-           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
-        {
-          VectorizedArrayType tmp;
-          tmp.load(vec_ptr);
-          tmp += dof_values[i];
-          tmp.store(vec_ptr);
-        }
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            VectorType &         vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, false>) const
-    {
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          vector_access_add(vec,
-                            dof_index + v +
-                              i * VectorizedArrayType::n_array_elements,
-                            dof_values[i][v]);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      VectorType &         vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, true>) const
-    {
-      vectorized_transpose_and_store(
-        true, dofs_per_cell, dof_values, dof_indices, vec.begin());
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      VectorType &         vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, false>) const
-    {
-      for (unsigned int d = 0; d < dofs_per_cell; ++d)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]);
-    }
-
-
-
-    // variant where VectorType::value_type is the same as Number -> can call
-    // scatter
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       VectorType &         vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, true>) const
-    {
-#  if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        vector_access(vec, indices[v] + constant_offset) += res[v];
-#  else
-      // only use gather in case there is also scatter.
-      VectorizedArrayType tmp;
-      tmp.gather(vec.begin() + constant_offset, indices);
-      tmp += res;
-      tmp.scatter(indices, vec.begin() + constant_offset);
-#  endif
-    }
-
-
-
-    // variant where VectorType::value_type is not the same as Number -> must
-    // manually append all data
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       VectorType &         vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, false>) const
-    {
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        vector_access_add(vec, indices[v] + constant_offset, res[v]);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dof_global(const types::global_dof_index index,
-                       VectorType &                  vec,
-                       Number &                      res) const
-    {
-      vector_access_add_global(vec, index, res);
-    }
-
-
-
-    void
-    pre_constraints(const Number &input, Number &res) const
-    {
-      res = input;
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_constraint(const unsigned int index,
-                       const Number       weight,
-                       VectorType &       vec,
-                       Number &           res) const
-    {
-      vector_access_add(vec, index, weight * res);
-    }
-
-
-
-    void
-    post_constraints(const Number &, Number &) const
-    {}
-
-
-
-    void
-    process_empty(VectorizedArrayType &) const
-    {}
-  };
-
-
-
-  // 3. A class to set elements of the vector
-  template <typename Number, typename VectorizedArrayType>
-  struct VectorSetter
-  {
-    template <typename VectorType>
-    void
-    process_dof(const unsigned int index, VectorType &vec, Number &res) const
-    {
-      vector_access(vec, index) = res;
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            VectorType &         vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, true>) const
-    {
-      Number *vec_ptr = vec.begin() + dof_index;
-      for (unsigned int i = 0; i < dofs_per_cell;
-           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
-        dof_values[i].store(vec_ptr);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized(const unsigned int   dofs_per_cell,
-                            const unsigned int   dof_index,
-                            VectorType &         vec,
-                            VectorizedArrayType *dof_values,
-                            std::integral_constant<bool, false>) const
-    {
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          vector_access(vec,
-                        dof_index + v +
-                          i * VectorizedArrayType::n_array_elements) =
-            dof_values[i][v];
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      VectorType &         vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, true>) const
-    {
-      vectorized_transpose_and_store(
-        false, dofs_per_cell, dof_values, dof_indices, vec.begin());
-    }
-
-
-
-    template <typename VectorType, bool booltype>
-    void
-    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
-                                      const unsigned int * dof_indices,
-                                      VectorType &         vec,
-                                      VectorizedArrayType *dof_values,
-                                      std::integral_constant<bool, false>) const
-    {
-      for (unsigned int i = 0; i < dofs_per_cell; ++i)
-        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-          vector_access(vec, dof_indices[v] + i) = dof_values[i][v];
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       VectorType &         vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, true>) const
-    {
-      res.scatter(indices, vec.begin() + constant_offset);
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dof_gather(const unsigned int * indices,
-                       VectorType &         vec,
-                       const unsigned int   constant_offset,
-                       VectorizedArrayType &res,
-                       std::integral_constant<bool, false>) const
-    {
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        vector_access(vec, indices[v] + constant_offset) = res[v];
-    }
-
-
-
-    template <typename VectorType>
-    void
-    process_dof_global(const types::global_dof_index index,
-                       VectorType &                  vec,
-                       Number &                      res) const
-    {
-      vec(index) = res;
-    }
-
-
-
-    void
-    pre_constraints(const Number &, Number &) const
-    {}
-
-
-
-    template <typename VectorType>
-    void
-    process_constraint(const unsigned int,
-                       const Number,
-                       VectorType &,
-                       Number &) const
-    {}
-
-
-
-    void
-    post_constraints(const Number &, Number &) const
-    {}
-
-
-
-    void
-    process_empty(VectorizedArrayType &) const
-    {}
-  };
-
-
-
   // allows to select between block vectors and non-block vectors, which
   // allows to use a unified interface for extracting blocks on block vectors
   // and doing nothing on usual vectors
@@ -8017,80 +7359,25 @@ FEFaceEvaluation<dim,
   if (!(evaluate_values + evaluate_gradients))
     return;
 
-  constexpr unsigned int static_dofs_per_face =
-    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                     numbers::invalid_unsigned_int;
-  const unsigned int dofs_per_face =
-    fe_degree > -1 ? static_dofs_per_face :
-                     Utilities::pow(this->data->fe_degree + 1, dim - 1);
-
-  // we allocate small amounts of data on the stack to signal the compiler
-  // that this temporary data is only needed for the calculations but the
-  // final results can be discarded and need not be written back to
-  // memory. For large sizes or when the dofs per face is not a compile-time
-  // constant, however, we want to go to the heap in the `scratch_data`
-  // variable to not risk a stack overflow.
-  constexpr unsigned int stack_array_size_threshold = 100;
-
-  VectorizedArrayType
-                       temp_data[static_dofs_per_face < stack_array_size_threshold ?
-                n_components * 2 * static_dofs_per_face :
-                1];
-  VectorizedArrayType *temp1;
-  if (static_dofs_per_face < stack_array_size_threshold)
-    temp1 = &temp_data[0];
-  else
-    temp1 = this->scratch_data;
-
-  internal::FEFaceNormalEvaluationImpl<
+  internal::FEFaceEvaluationSelector<
     dim,
     fe_degree,
+    n_q_points_1d,
     n_components,
-    VectorizedArrayType>::template interpolate<true, false>(*this->data,
-                                                            values_array,
-                                                            temp1,
-                                                            evaluate_gradients,
-                                                            this->face_no);
-
-  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
-  if (fe_degree > -1 &&
-      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-      this->data->element_type <=
-        internal::MatrixFreeFunctions::tensor_symmetric)
-    internal::FEFaceEvaluationImpl<
-      true,
-      dim,
-      fe_degree,
-      n_q_points_1d_actual,
-      n_components,
-      VectorizedArrayType>::evaluate_in_face(*this->data,
-                                             temp1,
-                                             this->begin_values(),
-                                             this->begin_gradients(),
-                                             this->scratch_data +
-                                               2 * n_components * dofs_per_face,
-                                             evaluate_values,
-                                             evaluate_gradients,
-                                             this->subface_index);
-  else
-    internal::FEFaceEvaluationImpl<
-      false,
-      dim,
-      fe_degree,
-      n_q_points_1d_actual,
-      n_components,
-      VectorizedArrayType>::evaluate_in_face(*this->data,
-                                             temp1,
-                                             this->begin_values(),
-                                             this->begin_gradients(),
-                                             this->scratch_data +
-                                               2 * n_components * dofs_per_face,
-                                             evaluate_values,
-                                             evaluate_gradients,
-                                             this->subface_index);
-
-  if (this->face_orientation)
-    adjust_for_face_orientation(false, evaluate_values, evaluate_gradients);
+    Number,
+    VectorizedArrayType>::evaluate(*this->data,
+                                   values_array,
+                                   this->begin_values(),
+                                   this->begin_gradients(),
+                                   this->scratch_data,
+                                   evaluate_values,
+                                   evaluate_gradients,
+                                   this->face_no,
+                                   this->subface_index,
+                                   this->face_orientation,
+                                   this->mapping_data
+                                     ->descriptor[this->active_fe_index]
+                                     .face_orientations);
 
 #  ifdef DEBUG
   if (evaluate_values == true)
@@ -8146,73 +7433,25 @@ FEFaceEvaluation<dim,
   if (!(integrate_values + integrate_gradients))
     return;
 
-  if (this->face_orientation)
-    adjust_for_face_orientation(true, integrate_values, integrate_gradients);
-
-  constexpr unsigned int static_dofs_per_face =
-    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                     numbers::invalid_unsigned_int;
-  const unsigned int dofs_per_face =
-    fe_degree > -1 ? static_dofs_per_face :
-                     Utilities::pow(this->data->fe_degree + 1, dim - 1);
-
-  constexpr unsigned int stack_array_size_threshold = 100;
-
-  VectorizedArrayType
-                       temp_data[static_dofs_per_face < stack_array_size_threshold ?
-                n_components * 2 * static_dofs_per_face :
-                1];
-  VectorizedArrayType *temp1;
-  if (static_dofs_per_face < stack_array_size_threshold)
-    temp1 = &temp_data[0];
-  else
-    temp1 = this->scratch_data;
-
-  const unsigned int n_q_points_1d_actual = fe_degree > -1 ? n_q_points_1d : 0;
-  if (fe_degree > -1 &&
-      this->subface_index >= GeometryInfo<dim - 1>::max_children_per_cell &&
-      this->data->element_type <=
-        internal::MatrixFreeFunctions::tensor_symmetric)
-    internal::FEFaceEvaluationImpl<
-      true,
-      dim,
-      fe_degree,
-      n_q_points_1d_actual,
-      n_components,
-      VectorizedArrayType>::integrate_in_face(*this->data,
-                                              temp1,
-                                              this->begin_values(),
-                                              this->begin_gradients(),
-                                              this->scratch_data +
-                                                2 * n_components *
-                                                  dofs_per_face,
-                                              integrate_values,
-                                              integrate_gradients,
-                                              this->subface_index);
-  else
-    internal::FEFaceEvaluationImpl<
-      false,
-      dim,
-      fe_degree,
-      n_q_points_1d_actual,
-      n_components,
-      VectorizedArrayType>::integrate_in_face(*this->data,
-                                              temp1,
-                                              this->begin_values(),
-                                              this->begin_gradients(),
-                                              this->scratch_data +
-                                                2 * n_components *
-                                                  dofs_per_face,
-                                              integrate_values,
-                                              integrate_gradients,
-                                              this->subface_index);
-
-  internal::FEFaceNormalEvaluationImpl<dim,
-                                       fe_degree,
-                                       n_components,
-                                       VectorizedArrayType>::
-    template interpolate<false, false>(
-      *this->data, temp1, values_array, integrate_gradients, this->face_no);
+  internal::FEFaceEvaluationSelector<
+    dim,
+    fe_degree,
+    n_q_points_1d,
+    n_components,
+    Number,
+    VectorizedArrayType>::integrate(*this->data,
+                                    values_array,
+                                    this->begin_values(),
+                                    this->begin_gradients(),
+                                    this->scratch_data,
+                                    integrate_values,
+                                    integrate_gradients,
+                                    this->face_no,
+                                    this->subface_index,
+                                    this->face_orientation,
+                                    this->mapping_data
+                                      ->descriptor[this->active_fe_index]
+                                      .face_orientations);
 }
 
 
@@ -8235,502 +7474,43 @@ FEFaceEvaluation<
                                         const bool        evaluate_values,
                                         const bool        evaluate_gradients)
 {
-  const unsigned int side = this->face_no % 2;
-
-  constexpr unsigned int static_dofs_per_face =
-    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                     numbers::invalid_unsigned_int;
-  const unsigned int dofs_per_face =
-    fe_degree > -1 ? static_dofs_per_face :
-                     Utilities::pow(this->data->fe_degree + 1, dim - 1);
-
-  constexpr unsigned int stack_array_size_threshold = 100;
-
-  VectorizedArrayType
-    temp_data[static_dofs_per_face < stack_array_size_threshold ?
-                n_components_ * 2 * dofs_per_face :
-                1];
-  VectorizedArrayType *__restrict temp1;
-  if (static_dofs_per_face < stack_array_size_threshold)
-    temp1 = &temp_data[0];
-  else
-    temp1 = this->scratch_data;
-
-  internal::VectorReader<Number, VectorizedArrayType> reader;
-  std::integral_constant<bool,
-                         internal::is_vectorizable<VectorType, Number>::value>
-    vector_selector;
-
-  // case 1: contiguous and interleaved indices
-  if (((evaluate_gradients == false &&
-        this->data->nodal_at_cell_boundaries == true) ||
-       (this->data->element_type ==
-          internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-        fe_degree > 1)) &&
-      this->dof_info
-          ->index_storage_variants[this->dof_access_index][this->cell] ==
-        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-          interleaved_contiguous)
-    {
-      AssertDimension(
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
-        VectorizedArrayType::n_array_elements);
-      const unsigned int dof_index =
-        this->dof_info
-          ->dof_indices_contiguous[this->dof_access_index]
-                                  [this->cell *
-                                   VectorizedArrayType::n_array_elements] +
-        this->dof_info
-            ->component_dof_indices_offset[this->active_fe_index]
-                                          [this->first_selected_component] *
-          VectorizedArrayType::n_array_elements;
-
-      if (fe_degree > 1 && evaluate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 1 + side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 = index_array[2 * i];
-              const unsigned int ind2 = index_array[2 * i + 1];
-              AssertIndexRange(ind1, dofs_per_cell);
-              AssertIndexRange(ind2, dofs_per_cell);
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  reader.process_dofs_vectorized(
-                    1,
-                    dof_index + (ind1 + comp * static_dofs_per_component) *
-                                  VectorizedArrayType::n_array_elements,
-                    input_vector,
-                    temp1 + i + 2 * comp * dofs_per_face,
-                    vector_selector);
-                  reader.process_dofs_vectorized(
-                    1,
-                    dof_index + (ind2 + comp * static_dofs_per_component) *
-                                  VectorizedArrayType::n_array_elements,
-                    input_vector,
-                    temp1 + dofs_per_face + i + 2 * comp * dofs_per_face,
-                    vector_selector);
-                  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                    grad_weight *
-                    (temp1[i + 2 * comp * dofs_per_face] -
-                     temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind = index_array[i];
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                reader.process_dofs_vectorized(
-                  1,
-                  dof_index + (ind + comp * static_dofs_per_component) *
-                                VectorizedArrayType::n_array_elements,
-                  input_vector,
-                  temp1 + i + 2 * comp * dofs_per_face,
-                  vector_selector);
-            }
-        }
-    }
-
-  // case 2: contiguous and interleaved indices with fixed stride
-  else if (((evaluate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               interleaved_contiguous_strided)
-    {
-      AssertDimension(
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
-        VectorizedArrayType::n_array_elements);
-      const unsigned int *indices =
-        &this->dof_info
-           ->dof_indices_contiguous[this->dof_access_index]
-                                   [this->cell *
-                                    VectorizedArrayType::n_array_elements];
-      if (fe_degree > 1 && evaluate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 1 + side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 =
-                index_array[2 * i] * VectorizedArrayType::n_array_elements;
-              const unsigned int ind2 =
-                index_array[2 * i + 1] * VectorizedArrayType::n_array_elements;
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  reader.process_dof_gather(
-                    indices,
-                    input_vector,
-                    ind1 +
-                      comp * static_dofs_per_component *
-                        VectorizedArrayType::n_array_elements +
-                      this->dof_info->component_dof_indices_offset
-                          [this->active_fe_index]
-                          [this->first_selected_component] *
-                        VectorizedArrayType::n_array_elements,
-                    temp1[i + 2 * comp * dofs_per_face],
-                    vector_selector);
-                  VectorizedArrayType grad;
-                  reader.process_dof_gather(
-                    indices,
-                    input_vector,
-                    ind2 +
-                      comp * static_dofs_per_component *
-                        VectorizedArrayType::n_array_elements +
-                      this->dof_info->component_dof_indices_offset
-                          [this->active_fe_index]
-                          [this->first_selected_component] *
-                        VectorizedArrayType::n_array_elements,
-                    grad,
-                    vector_selector);
-                  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                    grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind =
-                index_array[i] * VectorizedArrayType::n_array_elements;
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                reader.process_dof_gather(
-                  indices,
-                  input_vector,
-                  ind +
-                    comp * static_dofs_per_component *
-                      VectorizedArrayType::n_array_elements +
-                    this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index]
-                        [this->first_selected_component] *
-                      VectorizedArrayType::n_array_elements,
-                  temp1[i + 2 * comp * dofs_per_face],
-                  vector_selector);
-            }
-        }
-    }
-
-  // case 3: contiguous and interleaved indices with mixed stride
-  else if (((evaluate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               interleaved_contiguous_mixed_strides)
-    {
-      const unsigned int *strides =
-        &this->dof_info->dof_indices_interleave_strides
-           [this->dof_access_index]
-           [this->cell * VectorizedArrayType::n_array_elements];
-      unsigned int indices[VectorizedArrayType::n_array_elements];
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        indices[v] =
-          this->dof_info->dof_indices_contiguous
-            [this->dof_access_index]
-            [this->cell * VectorizedArrayType::n_array_elements + v] +
-          this->dof_info
-              ->component_dof_indices_offset[this->active_fe_index]
-                                            [this->first_selected_component] *
-            strides[v];
-      const unsigned int nvec =
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
-
-      if (fe_degree > 1 && evaluate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 1 + side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          if (nvec == VectorizedArrayType::n_array_elements)
-            for (unsigned int comp = 0; comp < n_components_; ++comp)
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  unsigned int ind1[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind1[v] = indices[v] + (comp * static_dofs_per_component +
-                                            index_array[2 * i]) *
-                                             strides[v];
-                  unsigned int ind2[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind2[v] = indices[v] + (comp * static_dofs_per_component +
-                                            index_array[2 * i + 1]) *
-                                             strides[v];
-                  reader.process_dof_gather(ind1,
-                                            input_vector,
-                                            0,
-                                            temp1[i + 2 * comp * dofs_per_face],
-                                            vector_selector);
-                  VectorizedArrayType grad;
-                  reader.process_dof_gather(
-                    ind2, input_vector, 0, grad, vector_selector);
-                  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                    grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
-                }
-          else
-            {
-              for (unsigned int i = 0; i < n_components_ * 2 * dofs_per_face;
-                   ++i)
-                temp1[i] = VectorizedArrayType();
-              for (unsigned int v = 0; v < nvec; ++v)
-                for (unsigned int comp = 0; comp < n_components_; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      const unsigned int ind1 =
-                        indices[v] + (comp * static_dofs_per_component +
-                                      index_array[2 * i]) *
-                                       strides[v];
-                      const unsigned int ind2 =
-                        indices[v] + (comp * static_dofs_per_component +
-                                      index_array[2 * i + 1]) *
-                                       strides[v];
-                      reader.process_dof(
-                        ind1,
-                        const_cast<VectorType &>(input_vector),
-                        temp1[i + 2 * comp * dofs_per_face][v]);
-                      Number grad;
-                      reader.process_dof(ind2,
-                                         const_cast<VectorType &>(input_vector),
-                                         grad);
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v] =
-                        grad_weight[0] *
-                        (temp1[i + 2 * comp * dofs_per_face][v] - grad);
-                    }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          if (nvec == VectorizedArrayType::n_array_elements)
-            for (unsigned int comp = 0; comp < n_components_; ++comp)
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  unsigned int ind[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind[v] = indices[v] + (comp * static_dofs_per_component +
-                                           index_array[i]) *
-                                            strides[v];
-                  reader.process_dof_gather(ind,
-                                            input_vector,
-                                            0,
-                                            temp1[i + 2 * comp * dofs_per_face],
-                                            vector_selector);
-                }
-          else
-            {
-              for (unsigned int i = 0; i < n_components_ * dofs_per_face; ++i)
-                temp1[i] = VectorizedArrayType();
-              for (unsigned int v = 0; v < nvec; ++v)
-                for (unsigned int comp = 0; comp < n_components_; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      const unsigned int ind1 =
-                        indices[v] +
-                        (comp * static_dofs_per_component + index_array[i]) *
-                          strides[v];
-                      reader.process_dof(
-                        ind1,
-                        const_cast<VectorType &>(input_vector),
-                        temp1[i + 2 * comp * dofs_per_face][v]);
-                    }
-            }
-        }
-    }
-
-  // case 4: contiguous indices without interleaving
-  else if (((evaluate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               contiguous &&
-           this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
-                                                       [this->cell] ==
-             VectorizedArrayType::n_array_elements)
-    {
-      const unsigned int *indices =
-        &this->dof_info
-           ->dof_indices_contiguous[this->dof_access_index]
-                                   [this->cell *
-                                    VectorizedArrayType::n_array_elements];
-      if (evaluate_gradients == true &&
-          this->data->element_type ==
-            internal::MatrixFreeFunctions::tensor_symmetric_hermite)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 1 + side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 = index_array[2 * i];
-              const unsigned int ind2 = index_array[2 * i + 1];
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  reader.process_dof_gather(
-                    indices,
-                    input_vector,
-                    ind1 + comp * static_dofs_per_component +
-                      this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index][this->first_selected_component],
-                    temp1[i + 2 * comp * dofs_per_face],
-                    vector_selector);
-                  VectorizedArrayType grad;
-                  reader.process_dof_gather(
-                    indices,
-                    input_vector,
-                    ind2 + comp * static_dofs_per_component +
-                      this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index][this->first_selected_component],
-                    grad,
-                    vector_selector);
-                  temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                    grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            for (unsigned int comp = 0; comp < n_components_; ++comp)
-              {
-                const unsigned int ind = index_array[i];
-                reader.process_dof_gather(
-                  indices,
-                  input_vector,
-                  ind + comp * static_dofs_per_component +
-                    this->dof_info->component_dof_indices_offset
-                      [this->active_fe_index][this->first_selected_component],
-                  temp1[i + comp * 2 * dofs_per_face],
-                  vector_selector);
-              }
-        }
-    }
-
-  // case 5: default vector access
-  else
+  static_assert(internal::has_begin<VectorType>::value &&
+                  (std::is_same<decltype(std::declval<VectorType>().begin()),
+                                double *>::value ||
+                   std::is_same<decltype(std::declval<VectorType>().begin()),
+                                float *>::value),
+                "This function requires a vector type with begin() function "
+                "evaluating to a pointer to basic number (float,double). "
+                "Use read_dof_values() followed by evaluate() instead.");
+
+  if (!internal::FEFaceEvaluationSelector<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType>::
+        gather_evaluate(input_vector.begin(),
+                        *this->data,
+                        *this->dof_info,
+                        this->begin_values(),
+                        this->begin_gradients(),
+                        this->scratch_data,
+                        evaluate_values,
+                        evaluate_gradients,
+                        this->active_fe_index,
+                        this->first_selected_component,
+                        this->cell,
+                        this->face_no,
+                        this->subface_index,
+                        this->dof_access_index,
+                        this->face_orientation,
+                        this->mapping_data->descriptor[this->active_fe_index]
+                          .face_orientations))
     {
       this->read_dof_values(input_vector);
-      internal::FEFaceNormalEvaluationImpl<dim,
-                                           fe_degree,
-                                           n_components_,
-                                           VectorizedArrayType>::
-        template interpolate<true, false>(*this->data,
-                                          this->values_dofs[0],
-                                          temp1,
-                                          evaluate_gradients,
-                                          this->face_no);
+      this->evaluate(evaluate_values, evaluate_gradients);
     }
 
-  if (fe_degree > -1 &&
-      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-      this->data->element_type <=
-        internal::MatrixFreeFunctions::tensor_symmetric)
-    internal::FEFaceEvaluationImpl<
-      true,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      n_components_,
-      VectorizedArrayType>::evaluate_in_face(*this->data,
-                                             temp1,
-                                             this->values_quad[0],
-                                             this->gradients_quad[0][0],
-                                             this->scratch_data +
-                                               2 * n_components_ *
-                                                 dofs_per_face,
-                                             evaluate_values,
-                                             evaluate_gradients,
-                                             this->subface_index);
-  else
-    internal::FEFaceEvaluationImpl<
-      false,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      n_components_,
-      VectorizedArrayType>::evaluate_in_face(*this->data,
-                                             temp1,
-                                             this->values_quad[0],
-                                             this->gradients_quad[0][0],
-                                             this->scratch_data +
-                                               2 * n_components_ *
-                                                 dofs_per_face,
-                                             evaluate_values,
-                                             evaluate_gradients,
-                                             this->subface_index);
-
-  if (this->face_orientation)
-    adjust_for_face_orientation(false, evaluate_values, evaluate_gradients);
-
 #  ifdef DEBUG
   if (evaluate_values == true)
     this->values_quad_initialized = true;
@@ -8759,559 +7539,51 @@ FEFaceEvaluation<
                                           const bool  integrate_gradients,
                                           VectorType &destination)
 {
-  const unsigned int side = this->face_no % 2;
-  const unsigned int dofs_per_face =
-    fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                     Utilities::pow(this->data->fe_degree + 1, dim - 1);
-
-  constexpr unsigned int stack_array_size_threshold = 100;
-
-  VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ?
-                                  n_components_ * 2 * dofs_per_face :
-                                  1];
-  VectorizedArrayType *__restrict temp1;
-  if (dofs_per_face < stack_array_size_threshold)
-    temp1 = &temp_data[0];
-  else
-    temp1 = this->scratch_data;
-
-  if (this->face_orientation)
-    adjust_for_face_orientation(true, integrate_values, integrate_gradients);
-  if (fe_degree > -1 &&
-      this->subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-      this->data->element_type <=
-        internal::MatrixFreeFunctions::tensor_symmetric)
-    internal::FEFaceEvaluationImpl<
-      true,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      n_components_,
-      VectorizedArrayType>::integrate_in_face(*this->data,
-                                              temp1,
-                                              this->values_quad[0],
-                                              this->gradients_quad[0][0],
-                                              this->scratch_data +
-                                                2 * n_components_ *
-                                                  dofs_per_face,
-                                              integrate_values,
-                                              integrate_gradients,
-                                              this->subface_index);
-  else
-    internal::FEFaceEvaluationImpl<
-      false,
-      dim,
-      fe_degree,
-      n_q_points_1d,
-      n_components_,
-      VectorizedArrayType>::integrate_in_face(*this->data,
-                                              temp1,
-                                              this->values_quad[0],
-                                              this->gradients_quad[0][0],
-                                              this->scratch_data +
-                                                2 * n_components_ *
-                                                  dofs_per_face,
-                                              integrate_values,
-                                              integrate_gradients,
-                                              this->subface_index);
-
-#  ifdef DEBUG
-  this->dof_values_initialized = true;
-#  endif
-
-  internal::VectorDistributorLocalToGlobal<Number, VectorizedArrayType> writer;
-  std::integral_constant<bool,
-                         internal::is_vectorizable<VectorType, Number>::value>
-    vector_selector;
-
-  // case 1: contiguous and interleaved indices
-  if (((integrate_gradients == false &&
-        this->data->nodal_at_cell_boundaries == true) ||
-       (this->data->element_type ==
-          internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-        fe_degree > 1)) &&
-      this->dof_info
-          ->index_storage_variants[this->dof_access_index][this->cell] ==
-        internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-          interleaved_contiguous)
-    {
-      AssertDimension(
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
-        VectorizedArrayType::n_array_elements);
-      const unsigned int dof_index =
-        this->dof_info
-          ->dof_indices_contiguous[this->dof_access_index]
-                                  [this->cell *
-                                   VectorizedArrayType::n_array_elements] +
-        this->dof_info
-            ->component_dof_indices_offset[this->active_fe_index]
-                                          [this->first_selected_component] *
-          VectorizedArrayType::n_array_elements;
-
-      if (fe_degree > 1 && integrate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 2 - side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 = index_array[2 * i];
-              const unsigned int ind2 = index_array[2 * i + 1];
-              AssertIndexRange(ind1, dofs_per_cell);
-              AssertIndexRange(ind2, dofs_per_cell);
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  VectorizedArrayType val =
-                    temp1[i + 2 * comp * dofs_per_face] -
-                    grad_weight *
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  VectorizedArrayType grad =
-                    grad_weight *
-                    temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  writer.process_dofs_vectorized(
-                    1,
-                    dof_index + (ind1 + comp * static_dofs_per_component) *
-                                  VectorizedArrayType::n_array_elements,
-                    destination,
-                    &val,
-                    vector_selector);
-                  writer.process_dofs_vectorized(
-                    1,
-                    dof_index + (ind2 + comp * static_dofs_per_component) *
-                                  VectorizedArrayType::n_array_elements,
-                    destination,
-                    &grad,
-                    vector_selector);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind = index_array[i];
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                writer.process_dofs_vectorized(
-                  1,
-                  dof_index + (ind + comp * static_dofs_per_component) *
-                                VectorizedArrayType::n_array_elements,
-                  destination,
-                  temp1 + i + 2 * comp * dofs_per_face,
-                  vector_selector);
-            }
-        }
-    }
-
-  // case 2: contiguous and interleaved indices with fixed stride
-  else if (((integrate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               interleaved_contiguous_strided)
-    {
-      AssertDimension(
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell],
-        VectorizedArrayType::n_array_elements);
-      const unsigned int *indices =
-        &this->dof_info
-           ->dof_indices_contiguous[this->dof_access_index]
-                                   [this->cell *
-                                    VectorizedArrayType::n_array_elements];
-      if (fe_degree > 1 && integrate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 2 - side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 =
-                index_array[2 * i] * VectorizedArrayType::n_array_elements;
-              const unsigned int ind2 =
-                index_array[2 * i + 1] * VectorizedArrayType::n_array_elements;
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  VectorizedArrayType val =
-                    temp1[i + 2 * comp * dofs_per_face] -
-                    grad_weight *
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  VectorizedArrayType grad =
-                    grad_weight *
-                    temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  writer.process_dof_gather(
-                    indices,
-                    destination,
-                    ind1 +
-                      comp * static_dofs_per_component *
-                        VectorizedArrayType::n_array_elements +
-                      this->dof_info->component_dof_indices_offset
-                          [this->active_fe_index]
-                          [this->first_selected_component] *
-                        VectorizedArrayType::n_array_elements,
-                    val,
-                    vector_selector);
-                  writer.process_dof_gather(
-                    indices,
-                    destination,
-                    ind2 +
-                      comp * static_dofs_per_component *
-                        VectorizedArrayType::n_array_elements +
-                      this->dof_info->component_dof_indices_offset
-                          [this->active_fe_index]
-                          [this->first_selected_component] *
-                        VectorizedArrayType::n_array_elements,
-                    grad,
-                    vector_selector);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind =
-                index_array[i] * VectorizedArrayType::n_array_elements;
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                writer.process_dof_gather(
-                  indices,
-                  destination,
-                  ind +
-                    comp * static_dofs_per_component *
-                      VectorizedArrayType::n_array_elements +
-                    this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index]
-                        [this->first_selected_component] *
-                      VectorizedArrayType::n_array_elements,
-                  temp1[i + 2 * comp * dofs_per_face],
-                  vector_selector);
-            }
-        }
-    }
-
-  // case 3: contiguous and interleaved indices with mixed stride
-  else if (((integrate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               interleaved_contiguous_mixed_strides)
-    {
-      const unsigned int *strides =
-        &this->dof_info->dof_indices_interleave_strides
-           [this->dof_access_index]
-           [this->cell * VectorizedArrayType::n_array_elements];
-      unsigned int indices[VectorizedArrayType::n_array_elements];
-      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
-        indices[v] =
-          this->dof_info->dof_indices_contiguous
-            [this->dof_access_index]
-            [this->cell * VectorizedArrayType::n_array_elements + v] +
-          this->dof_info
-              ->component_dof_indices_offset[this->active_fe_index]
-                                            [this->first_selected_component] *
-            strides[v];
-      const unsigned int nvec =
-        this->dof_info
-          ->n_vectorization_lanes_filled[this->dof_access_index][this->cell];
-
-      if (fe_degree > 1 && integrate_gradients == true)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 2 - side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          if (nvec == VectorizedArrayType::n_array_elements)
-            for (unsigned int comp = 0; comp < n_components_; ++comp)
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  unsigned int ind1[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind1[v] = indices[v] + (comp * static_dofs_per_component +
-                                            index_array[2 * i]) *
-                                             strides[v];
-                  unsigned int ind2[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind2[v] = indices[v] + (comp * static_dofs_per_component +
-                                            index_array[2 * i + 1]) *
-                                             strides[v];
-                  VectorizedArrayType val =
-                    temp1[i + 2 * comp * dofs_per_face] -
-                    grad_weight *
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  VectorizedArrayType grad =
-                    grad_weight *
-                    temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  writer.process_dof_gather(
-                    ind1, destination, 0, val, vector_selector);
-                  writer.process_dof_gather(
-                    ind2, destination, 0, grad, vector_selector);
-                }
-          else
-            {
-              for (unsigned int v = 0; v < nvec; ++v)
-                for (unsigned int comp = 0; comp < n_components_; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      const unsigned int ind1 =
-                        indices[v] + (comp * static_dofs_per_component +
-                                      index_array[2 * i]) *
-                                       strides[v];
-                      const unsigned int ind2 =
-                        indices[v] + (comp * static_dofs_per_component +
-                                      index_array[2 * i + 1]) *
-                                       strides[v];
-                      Number val =
-                        temp1[i + 2 * comp * dofs_per_face][v] -
-                        grad_weight[0] * temp1[i + dofs_per_face +
-                                               2 * comp * dofs_per_face][v];
-                      Number grad =
-                        grad_weight[0] *
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face][v];
-                      writer.process_dof(ind1, destination, val);
-                      writer.process_dof(ind2, destination, grad);
-                    }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          if (nvec == VectorizedArrayType::n_array_elements)
-            for (unsigned int comp = 0; comp < n_components_; ++comp)
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  unsigned int ind[VectorizedArrayType::n_array_elements];
-                  DEAL_II_OPENMP_SIMD_PRAGMA
-                  for (unsigned int v = 0;
-                       v < VectorizedArrayType::n_array_elements;
-                       ++v)
-                    ind[v] = indices[v] + (comp * static_dofs_per_component +
-                                           index_array[i]) *
-                                            strides[v];
-                  writer.process_dof_gather(ind,
-                                            destination,
-                                            0,
-                                            temp1[i + 2 * comp * dofs_per_face],
-                                            vector_selector);
-                }
-          else
-            {
-              for (unsigned int v = 0; v < nvec; ++v)
-                for (unsigned int comp = 0; comp < n_components_; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      const unsigned int ind1 =
-                        indices[v] +
-                        (comp * static_dofs_per_component + index_array[i]) *
-                          strides[v];
-                      writer.process_dof(
-                        ind1,
-                        destination,
-                        temp1[i + 2 * comp * dofs_per_face][v]);
-                    }
-            }
-        }
-    }
-
-  // case 4: contiguous indices without interleaving
-  else if (((integrate_gradients == false &&
-             this->data->nodal_at_cell_boundaries == true) ||
-            (this->data->element_type ==
-               internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
-             fe_degree > 1)) &&
-           this->dof_info
-               ->index_storage_variants[this->dof_access_index][this->cell] ==
-             internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-               contiguous &&
-           this->dof_info->n_vectorization_lanes_filled[this->dof_access_index]
-                                                       [this->cell] ==
-             VectorizedArrayType::n_array_elements)
-    {
-      const unsigned int *indices =
-        &this->dof_info
-           ->dof_indices_contiguous[this->dof_access_index]
-                                   [this->cell *
-                                    VectorizedArrayType::n_array_elements];
-
-      if (integrate_gradients == true &&
-          this->data->element_type ==
-            internal::MatrixFreeFunctions::tensor_symmetric_hermite)
-        {
-          // we know that the gradient weights for the Hermite case on the
-          // right (side==1) are the negative from the value at the left
-          // (side==0), so we only read out one of them.
-          const VectorizedArrayType grad_weight =
-            this->data->shape_data_on_face[0][fe_degree + 2 - side];
-          AssertDimension(this->data->face_to_cell_index_hermite.size(1),
-                          2 * dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_hermite(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind1 = index_array[2 * i];
-              const unsigned int ind2 = index_array[2 * i + 1];
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                {
-                  VectorizedArrayType val =
-                    temp1[i + 2 * comp * dofs_per_face] -
-                    grad_weight *
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  VectorizedArrayType grad =
-                    grad_weight *
-                    temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                  writer.process_dof_gather(
-                    indices,
-                    destination,
-                    comp * static_dofs_per_component + ind1 +
-                      this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index][this->first_selected_component],
-                    val,
-                    vector_selector);
-                  writer.process_dof_gather(
-                    indices,
-                    destination,
-                    comp * static_dofs_per_component + ind2 +
-                      this->dof_info->component_dof_indices_offset
-                        [this->active_fe_index][this->first_selected_component],
-                    grad,
-                    vector_selector);
-                }
-            }
-        }
-      else
-        {
-          AssertDimension(this->data->face_to_cell_index_nodal.size(1),
-                          dofs_per_face);
-          const unsigned int *index_array =
-            &this->data->face_to_cell_index_nodal(this->face_no, 0);
-          for (unsigned int i = 0; i < dofs_per_face; ++i)
-            {
-              const unsigned int ind = index_array[i];
-              for (unsigned int comp = 0; comp < n_components_; ++comp)
-                writer.process_dof_gather(
-                  indices,
-                  destination,
-                  comp * static_dofs_per_component + ind +
-                    this->dof_info->component_dof_indices_offset
-                      [this->active_fe_index][this->first_selected_component],
-                  temp1[i + 2 * comp * dofs_per_face],
-                  vector_selector);
-            }
-        }
-    }
-
-  // case 5: default vector access
-  else
+  static_assert(internal::has_begin<VectorType>::value &&
+                  (std::is_same<decltype(std::declval<VectorType>().begin()),
+                                double *>::value ||
+                   std::is_same<decltype(std::declval<VectorType>().begin()),
+                                float *>::value),
+                "This function requires a vector type with begin() function "
+                "evaluating to a pointer to basic number (float,double). "
+                "Use integrate() followed by distribute_local_to_global() "
+                "instead.");
+
+  if (!internal::FEFaceEvaluationSelector<dim,
+                                          fe_degree,
+                                          n_q_points_1d,
+                                          n_components,
+                                          Number,
+                                          VectorizedArrayType>::
+        integrate_scatter(destination.begin(),
+                          *this->data,
+                          *this->dof_info,
+                          this->begin_dof_values(),
+                          this->begin_values(),
+                          this->begin_gradients(),
+                          this->scratch_data,
+                          integrate_values,
+                          integrate_gradients,
+                          this->active_fe_index,
+                          this->first_selected_component,
+                          this->cell,
+                          this->face_no,
+                          this->subface_index,
+                          this->dof_access_index,
+                          this->face_orientation,
+                          this->mapping_data->descriptor[this->active_fe_index]
+                            .face_orientations))
     {
-      internal::FEFaceNormalEvaluationImpl<dim,
-                                           fe_degree,
-                                           n_components_,
-                                           VectorizedArrayType>::
-        template interpolate<false, false>(*this->data,
-                                           temp1,
-                                           this->values_dofs[0],
-                                           integrate_gradients,
-                                           this->face_no);
+      // if we arrive here, writing into the destination vector did not succeed
+      // because some of the assumptions in integrate_scatter were not
+      // fulfilled (e.g. an element or degree that does not support direct
+      // writing), so we must do it here
       this->distribute_local_to_global(destination);
     }
 }
 
 
 
-template <int dim,
-          int fe_degree,
-          int n_q_points_1d,
-          int n_components,
-          typename Number,
-          typename VectorizedArrayType>
-inline void
-FEFaceEvaluation<
-  dim,
-  fe_degree,
-  n_q_points_1d,
-  n_components,
-  Number,
-  VectorizedArrayType>::adjust_for_face_orientation(const bool integrate,
-                                                    const bool values,
-                                                    const bool gradients)
-{
-  VectorizedArrayType *tmp_values = this->scratch_data;
-  const unsigned int * orientations =
-    &this->mapping_data->descriptor[this->active_fe_index]
-       .face_orientations[this->face_orientation][0];
-  for (unsigned int c = 0; c < n_components; ++c)
-    {
-      if (values == true)
-        {
-          if (integrate)
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              tmp_values[q] = this->values_quad[c][orientations[q]];
-          else
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              tmp_values[orientations[q]] = this->values_quad[c][q];
-          for (unsigned int q = 0; q < n_q_points; ++q)
-            this->values_quad[c][q] = tmp_values[q];
-        }
-      if (gradients == true)
-        for (unsigned int d = 0; d < dim; ++d)
-          {
-            if (integrate)
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[q] = this->gradients_quad[c][d][orientations[q]];
-            else
-              for (unsigned int q = 0; q < n_q_points; ++q)
-                tmp_values[orientations[q]] = this->gradients_quad[c][d][q];
-            for (unsigned int q = 0; q < n_q_points; ++q)
-              this->gradients_quad[c][d][q] = tmp_values[q];
-          }
-    }
-}
-
-
-
 template <int dim,
           int fe_degree,
           int n_q_points_1d,
index bd4491fbb14ff804c139125230f91ab9cd6b159b..190681797818da9ebdc39a7365d6e9867ef65236 100644 (file)
@@ -1014,64 +1014,15 @@ namespace MatrixFreeOperators
           const VectorizedArrayType *               in_array,
           VectorizedArrayType *                     out_array) const
   {
-    constexpr unsigned int dofs_per_component =
-      Utilities::pow(fe_degree + 1, dim);
-    Assert(inverse_coefficients.size() > 0 &&
-             inverse_coefficients.size() % dofs_per_component == 0,
-           ExcMessage(
-             "Expected diagonal to be a multiple of scalar dof per cells"));
-    if (inverse_coefficients.size() != dofs_per_component)
-      AssertDimension(n_actual_components * dofs_per_component,
-                      inverse_coefficients.size());
-
-    Assert(dim >= 1 || dim <= 3, ExcNotImplemented());
-
-    internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
-                                     dim,
-                                     fe_degree + 1,
-                                     fe_degree + 1,
-                                     VectorizedArrayType>
-      evaluator(inverse_shape, inverse_shape, inverse_shape);
-
-    const unsigned int shift_coefficient =
-      inverse_coefficients.size() > dofs_per_component ? dofs_per_component : 0;
-    const VectorizedArrayType *inv_coefficient = inverse_coefficients.data();
-    VectorizedArrayType        temp_data_field[dofs_per_component];
-    for (unsigned int d = 0; d < n_actual_components; ++d)
-      {
-        const VectorizedArrayType *in  = in_array + d * dofs_per_component;
-        VectorizedArrayType *      out = out_array + d * dofs_per_component;
-        // Need to select 'apply' method with hessian slot because values
-        // assume symmetries that do not exist in the inverse shapes
-        evaluator.template hessians<0, false, false>(in, temp_data_field);
-        if (dim > 1)
-          {
-            evaluator.template hessians<1, false, false>(temp_data_field, out);
-
-            if (dim == 3)
-              {
-                evaluator.template hessians<2, false, false>(out,
-                                                             temp_data_field);
-                for (unsigned int q = 0; q < dofs_per_component; ++q)
-                  temp_data_field[q] *= inv_coefficient[q];
-                evaluator.template hessians<2, true, false>(temp_data_field,
-                                                            out);
-              }
-            else if (dim == 2)
-              for (unsigned int q = 0; q < dofs_per_component; ++q)
-                out[q] *= inv_coefficient[q];
-
-            evaluator.template hessians<1, true, false>(out, temp_data_field);
-          }
-        else
-          {
-            for (unsigned int q = 0; q < dofs_per_component; ++q)
-              temp_data_field[q] *= inv_coefficient[q];
-          }
-        evaluator.template hessians<0, true, false>(temp_data_field, out);
-
-        inv_coefficient += shift_coefficient;
-      }
+    internal::CellwiseInverseMassMatrixImpl<
+      dim,
+      fe_degree,
+      n_components,
+      VectorizedArrayType>::apply(inverse_shape,
+                                  inverse_coefficients,
+                                  n_actual_components,
+                                  in_array,
+                                  out_array);
   }
 
 
@@ -1091,35 +1042,14 @@ namespace MatrixFreeOperators
                                      const VectorizedArrayType *in_array,
                                      VectorizedArrayType *      out_array) const
   {
-    constexpr unsigned int dofs_per_cell = Utilities::pow(fe_degree + 1, dim);
-    internal::EvaluatorTensorProduct<internal::evaluate_evenodd,
-                                     dim,
-                                     fe_degree + 1,
-                                     fe_degree + 1,
-                                     VectorizedArrayType>
-      evaluator(AlignedVector<VectorizedArrayType>(),
-                AlignedVector<VectorizedArrayType>(),
-                inverse_shape);
-
-    for (unsigned int d = 0; d < n_actual_components; ++d)
-      {
-        const VectorizedArrayType *in  = in_array + d * dofs_per_cell;
-        VectorizedArrayType *      out = out_array + d * dofs_per_cell;
-
-        if (dim == 3)
-          {
-            evaluator.template hessians<2, true, false>(in, out);
-            evaluator.template hessians<1, true, false>(out, out);
-            evaluator.template hessians<0, true, false>(out, out);
-          }
-        if (dim == 2)
-          {
-            evaluator.template hessians<1, true, false>(in, out);
-            evaluator.template hessians<0, true, false>(out, out);
-          }
-        if (dim == 1)
-          evaluator.template hessians<0, true, false>(in, out);
-      }
+    internal::CellwiseInverseMassMatrixImpl<dim,
+                                            fe_degree,
+                                            n_components,
+                                            VectorizedArrayType>::
+      transform_from_q_points_to_basis(inverse_shape,
+                                       n_actual_components,
+                                       in_array,
+                                       out_array);
   }
 
 
diff --git a/include/deal.II/matrix_free/vector_access_internal.h b/include/deal.II/matrix_free/vector_access_internal.h
new file mode 100644 (file)
index 0000000..5275f90
--- /dev/null
@@ -0,0 +1,683 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#ifndef dealii_matrix_free_vector_access_internal_h
+#define dealii_matrix_free_vector_access_internal_h
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/base/vectorization.h>
+
+#include <deal.II/matrix_free/dof_info.h>
+#include <deal.II/matrix_free/type_traits.h>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  // below we use type-traits from matrix-free/type_traits.h
+
+  // access to generic const vectors that have operator ().
+  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
+  // where we should first do Partitioner::local_to_global()
+  template <typename VectorType,
+            typename std::enable_if<!has_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline typename VectorType::value_type
+  vector_access(const VectorType &vec, const unsigned int entry)
+  {
+    return vec(entry);
+  }
+
+
+
+  // access to generic non-const vectors that have operator ().
+  // FIXME: this is wrong for Trilinos/Petsc MPI vectors
+  // where we should first do Partitioner::local_to_global()
+  template <typename VectorType,
+            typename std::enable_if<!has_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline typename VectorType::value_type &
+  vector_access(VectorType &vec, const unsigned int entry)
+  {
+    return vec(entry);
+  }
+
+
+
+  // access to distributed MPI vectors that have a local_element(uint)
+  // method to access data in local index space, which is what we use in
+  // DoFInfo and hence in read_dof_values etc.
+  template <typename VectorType,
+            typename std::enable_if<has_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline typename VectorType::value_type &
+  vector_access(VectorType &vec, const unsigned int entry)
+  {
+    return vec.local_element(entry);
+  }
+
+
+
+  // same for const access
+  template <typename VectorType,
+            typename std::enable_if<has_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline typename VectorType::value_type
+  vector_access(const VectorType &vec, const unsigned int entry)
+  {
+    return vec.local_element(entry);
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<has_add_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_add(VectorType &                           vec,
+                    const unsigned int                     entry,
+                    const typename VectorType::value_type &val)
+  {
+    vec.add_local_element(entry, val);
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<!has_add_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_add(VectorType &                           vec,
+                    const unsigned int                     entry,
+                    const typename VectorType::value_type &val)
+  {
+    vector_access(vec, entry) += val;
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<has_add_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_add_global(VectorType &                           vec,
+                           const types::global_dof_index          entry,
+                           const typename VectorType::value_type &val)
+  {
+    vec.add(entry, val);
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<!has_add_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_add_global(VectorType &                           vec,
+                           const types::global_dof_index          entry,
+                           const typename VectorType::value_type &val)
+  {
+    vec(entry) += val;
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<has_set_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_set(VectorType &                           vec,
+                    const unsigned int                     entry,
+                    const typename VectorType::value_type &val)
+  {
+    vec.set_local_element(entry, val);
+  }
+
+
+
+  template <typename VectorType,
+            typename std::enable_if<!has_set_local_element<VectorType>::value,
+                                    VectorType>::type * = nullptr>
+  inline void
+  vector_access_set(VectorType &                           vec,
+                    const unsigned int                     entry,
+                    const typename VectorType::value_type &val)
+  {
+    vector_access(vec, entry) = val;
+  }
+
+
+
+  // this is to make sure that the parallel partitioning in VectorType
+  // is really the same as stored in MatrixFree.
+  // version below is when has_partitioners_are_compatible == false
+  // FIXME: this is incorrect for PETSc/Trilinos MPI vectors
+  template <
+    typename VectorType,
+    typename std::enable_if<!has_partitioners_are_compatible<VectorType>::value,
+                            VectorType>::type * = nullptr>
+  inline void
+  check_vector_compatibility(
+    const VectorType &                            vec,
+    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
+  {
+    (void)vec;
+    (void)dof_info;
+
+    AssertDimension(vec.size(), dof_info.vector_partitioner->size());
+  }
+
+
+
+  // same as above for has_partitioners_are_compatible == true
+  template <
+    typename VectorType,
+    typename std::enable_if<has_partitioners_are_compatible<VectorType>::value,
+                            VectorType>::type * = nullptr>
+  inline void
+  check_vector_compatibility(
+    const VectorType &                            vec,
+    const internal::MatrixFreeFunctions::DoFInfo &dof_info)
+  {
+    (void)vec;
+    (void)dof_info;
+    Assert(vec.partitioners_are_compatible(*dof_info.vector_partitioner),
+           ExcMessage(
+             "The parallel layout of the given vector is not "
+             "compatible with the parallel partitioning in MatrixFree. "
+             "Use MatrixFree::initialize_dof_vector to get a "
+             "compatible vector."));
+  }
+
+
+
+  // Below, three classes (VectorReader, VectorSetter,
+  // VectorDistributorLocalToGlobal) implement the same interface and can be
+  // used to to read from vector, set elements of a vector and add to elements
+  // of the vector.
+
+  // 1. A class to read data from vector
+  template <typename Number, typename VectorizedArrayType>
+  struct VectorReader
+  {
+    template <typename VectorType>
+    void
+    process_dof(const unsigned int index,
+                const VectorType & vec,
+                Number &           res) const
+    {
+      res = vector_access(vec, index);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            VectorType &         vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, true>) const
+    {
+      const Number *vec_ptr = vec.begin() + dof_index;
+      for (unsigned int i = 0; i < dofs_per_cell;
+           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
+        dof_values[i].load(vec_ptr);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            const VectorType &   vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          dof_values[i][v] = vector_access(
+            vec, dof_index + v + i * VectorizedArrayType::n_array_elements);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      VectorType &         vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, true>) const
+    {
+      dealii::vectorized_load_and_transpose(dofs_per_cell,
+                                            vec.begin(),
+                                            dof_indices,
+                                            dof_values);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      const VectorType &   vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, false>) const
+    {
+      for (unsigned int d = 0; d < dofs_per_cell; ++d)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          dof_values[d][v] = vector_access(vec, dof_indices[v] + d);
+    }
+
+
+
+    // variant where VectorType::value_type is the same as Number -> can call
+    // gather
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       VectorType &         vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, true>) const
+    {
+      res.gather(vec.begin() + constant_offset, indices);
+    }
+
+
+
+    // variant where VectorType::value_type is not the same as Number -> must
+    // manually load the data
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       const VectorType &   vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, false>) const
+    {
+      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+        res[v] = vector_access(vec, indices[v] + constant_offset);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dof_global(const types::global_dof_index index,
+                       const VectorType &            vec,
+                       Number &                      res) const
+    {
+      res = vec(index);
+    }
+
+
+
+    void
+    pre_constraints(const Number &, Number &res) const
+    {
+      res = Number();
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_constraint(const unsigned int index,
+                       const Number       weight,
+                       const VectorType & vec,
+                       Number &           res) const
+    {
+      res += weight * vector_access(vec, index);
+    }
+
+
+
+    void
+    post_constraints(const Number &sum, Number &write_pos) const
+    {
+      write_pos = sum;
+    }
+
+
+
+    void
+    process_empty(VectorizedArrayType &res) const
+    {
+      res = VectorizedArrayType();
+    }
+  };
+
+
+
+  // 2. A class to add values to the vector during
+  // FEEvaluation::distribute_local_to_global() call
+  template <typename Number, typename VectorizedArrayType>
+  struct VectorDistributorLocalToGlobal
+  {
+    template <typename VectorType>
+    void
+    process_dof(const unsigned int index, VectorType &vec, Number &res) const
+    {
+      vector_access_add(vec, index, res);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            VectorType &         vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, true>) const
+    {
+      Number *vec_ptr = vec.begin() + dof_index;
+      for (unsigned int i = 0; i < dofs_per_cell;
+           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
+        {
+          VectorizedArrayType tmp;
+          tmp.load(vec_ptr);
+          tmp += dof_values[i];
+          tmp.store(vec_ptr);
+        }
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            VectorType &         vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          vector_access_add(vec,
+                            dof_index + v +
+                              i * VectorizedArrayType::n_array_elements,
+                            dof_values[i][v]);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      VectorType &         vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, true>) const
+    {
+      vectorized_transpose_and_store(
+        true, dofs_per_cell, dof_values, dof_indices, vec.begin());
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      VectorType &         vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, false>) const
+    {
+      for (unsigned int d = 0; d < dofs_per_cell; ++d)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          vector_access_add(vec, dof_indices[v] + d, dof_values[d][v]);
+    }
+
+
+
+    // variant where VectorType::value_type is the same as Number -> can call
+    // scatter
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       VectorType &         vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, true>) const
+    {
+#if DEAL_II_COMPILER_VECTORIZATION_LEVEL < 3
+      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+        vector_access(vec, indices[v] + constant_offset) += res[v];
+#else
+      // only use gather in case there is also scatter.
+      VectorizedArrayType tmp;
+      tmp.gather(vec.begin() + constant_offset, indices);
+      tmp += res;
+      tmp.scatter(indices, vec.begin() + constant_offset);
+#endif
+    }
+
+
+
+    // variant where VectorType::value_type is not the same as Number -> must
+    // manually append all data
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       VectorType &         vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, false>) const
+    {
+      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+        vector_access_add(vec, indices[v] + constant_offset, res[v]);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dof_global(const types::global_dof_index index,
+                       VectorType &                  vec,
+                       Number &                      res) const
+    {
+      vector_access_add_global(vec, index, res);
+    }
+
+
+
+    void
+    pre_constraints(const Number &input, Number &res) const
+    {
+      res = input;
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_constraint(const unsigned int index,
+                       const Number       weight,
+                       VectorType &       vec,
+                       Number &           res) const
+    {
+      vector_access_add(vec, index, weight * res);
+    }
+
+
+
+    void
+    post_constraints(const Number &, Number &) const
+    {}
+
+
+
+    void
+    process_empty(VectorizedArrayType &) const
+    {}
+  };
+
+
+
+  // 3. A class to set elements of the vector
+  template <typename Number, typename VectorizedArrayType>
+  struct VectorSetter
+  {
+    template <typename VectorType>
+    void
+    process_dof(const unsigned int index, VectorType &vec, Number &res) const
+    {
+      vector_access(vec, index) = res;
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            VectorType &         vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, true>) const
+    {
+      Number *vec_ptr = vec.begin() + dof_index;
+      for (unsigned int i = 0; i < dofs_per_cell;
+           ++i, vec_ptr += VectorizedArrayType::n_array_elements)
+        dof_values[i].store(vec_ptr);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized(const unsigned int   dofs_per_cell,
+                            const unsigned int   dof_index,
+                            VectorType &         vec,
+                            VectorizedArrayType *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          vector_access(vec,
+                        dof_index + v +
+                          i * VectorizedArrayType::n_array_elements) =
+            dof_values[i][v];
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      VectorType &         vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, true>) const
+    {
+      vectorized_transpose_and_store(
+        false, dofs_per_cell, dof_values, dof_indices, vec.begin());
+    }
+
+
+
+    template <typename VectorType, bool booltype>
+    void
+    process_dofs_vectorized_transpose(const unsigned int   dofs_per_cell,
+                                      const unsigned int * dof_indices,
+                                      VectorType &         vec,
+                                      VectorizedArrayType *dof_values,
+                                      std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+          vector_access(vec, dof_indices[v] + i) = dof_values[i][v];
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       VectorType &         vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, true>) const
+    {
+      res.scatter(indices, vec.begin() + constant_offset);
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dof_gather(const unsigned int * indices,
+                       VectorType &         vec,
+                       const unsigned int   constant_offset,
+                       VectorizedArrayType &res,
+                       std::integral_constant<bool, false>) const
+    {
+      for (unsigned int v = 0; v < VectorizedArrayType::n_array_elements; ++v)
+        vector_access(vec, indices[v] + constant_offset) = res[v];
+    }
+
+
+
+    template <typename VectorType>
+    void
+    process_dof_global(const types::global_dof_index index,
+                       VectorType &                  vec,
+                       Number &                      res) const
+    {
+      vec(index) = res;
+    }
+
+
+
+    void
+    pre_constraints(const Number &, Number &) const
+    {}
+
+
+
+    template <typename VectorType>
+    void
+    process_constraint(const unsigned int,
+                       const Number,
+                       VectorType &,
+                       Number &) const
+    {}
+
+
+
+    void
+    post_constraints(const Number &, Number &) const
+    {}
+
+
+
+    void
+    process_empty(VectorizedArrayType &) const
+    {}
+  };
+} // namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.