]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add contiguous-interleaved storage option in MatrixFree.
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Mon, 6 Aug 2018 09:01:28 +0000 (11:01 +0200)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Tue, 7 Aug 2018 14:19:49 +0000 (16:19 +0200)
include/deal.II/matrix_free/dof_info.h
include/deal.II/matrix_free/dof_info.templates.h
include/deal.II/matrix_free/evaluation_selector.h
include/deal.II/matrix_free/fe_evaluation.h
include/deal.II/matrix_free/matrix_free.templates.h
source/matrix_free/evaluation_selector.inst.in

index 3c558f06d0fb78fe9256a7d3b1375aeea653d931..6839cd1d353da7ae290bd55930ed7cbe6a1830b1 100644 (file)
@@ -265,7 +265,56 @@ namespace internal
          * is directed to the array `dof_indices_contiguous` with the index
          * `cell_index*n_vectorization*n_components`.
          */
-        contiguous
+        contiguous,
+        /**
+         * This value indicates that the indices with a cell are contiguous and
+         * interleaved for vectorization, i.e., the first DoF index on a cell
+         * to the four or eight cells in the vectorization batch come first,
+         * than the second DoF index, and so on. Furthermore, the interleaving
+         * between cells implies that only the batches for vectorization can be
+         * accessed efficiently, whereas there is a strided access for getting
+         * only some of the entries.
+         *
+         * The two additional categories `interleaved_contiguous_strided` and
+         * `interleaved_contiguous_mixed_strides` are a consequence of this
+         * storage type. The former is for faces where at least one of the two
+         * adjacent sides will break with the interleaved storage. We then have
+         * to make a strided access as described in the next category. The last
+         * category `interleaved_contiguous_mixed_strides` appears in the ghost
+         * layer, see the more detailed description of that category below.
+         * Again, this is something that cannot be avoided in general once we
+         * interleave the indices between cells.
+         *
+         * For a cell/face of this index type, the data access in
+         * FEEvaluationBase is directed to the array `dof_indices_contiguous`
+         * with the index `cell_index*n_vectorization*n_components`.
+         */
+        interleaved_contiguous,
+        /**
+         * Similar to interleaved_contiguous storage, but for the case when the
+         * interleaved indices are only contiguous within the degrees of
+         * freedom, but not also over the components of a vectorized array.
+         * This happens typically on faces with DG where the cells have
+         * `interleaved_contiguous` storage but the faces' numbering is not the
+         * same as the cell's numbering. For a
+         * cell/face of this index type, the data access in FEEvaluationBase
+         * is directed to the array `dof_indices_contiguous` with the index
+         * `cell_index*n_vectorization*n_components`.
+         */
+        interleaved_contiguous_strided,
+        /**
+         * Similar to interleaved_contiguous_separate storage, but for the case
+         * when the interleaved indices are not `n_vectorization apart`. This
+         * happens typically within the ghost layer of DG where the remote
+         * owner has applied an interleaved storage and the current processor
+         * only sees some of the cells. For a
+         * cell/face of this index type, the data access in FEEvaluationBase
+         * is directed to the array `dof_indices_contiguous` with the index
+         * `cell_index*n_vectorization*n_components`, including the array
+         * `dof_indices_interleave_strides` for the information about the
+         * actual stride.
+         */
+        interleaved_contiguous_mixed_strides
       };
 
       /**
@@ -366,6 +415,16 @@ namespace internal
        */
       std::vector<unsigned int> dof_indices_contiguous[3];
 
+      /**
+       * Compressed index storage for faster access than through @p
+       * dof_indices used according to the description in IndexStorageVariants.
+       *
+       * The three arrays given here address the types for the faces decorated
+       * as minus (0), the faces decorated with as plus (1), and the cells
+       * (2).
+       */
+      std::vector<unsigned int> dof_indices_interleave_strides[3];
+
       /**
        * Caches the number of indices filled when vectorizing. This
        * information can implicitly deduced from the row_starts data fields,
index 844994c95b36bd2b9487c94d483520faf07d8031..caaba870dde0542462eb84b37ad56ac48dcb94c7 100644 (file)
@@ -167,6 +167,14 @@ namespace internal
       start_components.clear();
       row_starts_plain_indices.clear();
       plain_dof_indices.clear();
+      dof_indices_interleaved.clear();
+      for (unsigned int i = 0; i < 3; ++i)
+        {
+          index_storage_variants[i].clear();
+          dof_indices_contiguous[i].clear();
+          dof_indices_interleave_strides[i].clear();
+          n_vectorization_lanes_filled[i].clear();
+        }
       store_plain_indices = false;
       cell_active_fe_index.clear();
       max_fe_index = 0;
@@ -699,9 +707,14 @@ namespace internal
         numbers::invalid_unsigned_int);
       dof_indices_interleaved.resize(dof_indices.size(),
                                      numbers::invalid_unsigned_int);
+      dof_indices_interleave_strides[dof_access_cell].resize(
+        irregular_cells.size() * vectorization_length,
+        numbers::invalid_unsigned_int);
 
       std::vector<unsigned int> index_kinds(
-        static_cast<unsigned int>(IndexStorageVariants::contiguous) + 1);
+        static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous_mixed_strides) +
+        1);
       std::vector<unsigned int> offsets(vectorization_length);
       for (unsigned int i = 0; i < irregular_cells.size(); ++i)
         {
@@ -745,8 +758,10 @@ namespace internal
                         break;
                       }
                 }
+
               bool indices_are_interleaved_and_contiguous =
                 (ndofs > 1 && n_comp == vectorization_length);
+
               {
                 const unsigned int *dof_indices =
                   this->dof_indices.data() +
@@ -760,6 +775,7 @@ namespace internal
                         break;
                       }
               }
+
               if (indices_are_contiguous ||
                   indices_are_interleaved_and_contiguous)
                 {
@@ -772,52 +788,209 @@ namespace internal
                                             .first];
                 }
 
-              if (indices_are_contiguous)
+              if (indices_are_interleaved_and_contiguous)
+                {
+                  Assert(n_comp == vectorization_length, ExcInternalError());
+                  index_storage_variants[dof_access_cell][i] =
+                    IndexStorageVariants::interleaved_contiguous;
+                  for (unsigned int j = 0; j < n_comp; ++j)
+                    dof_indices_interleave_strides[2][i * vectorization_length +
+                                                      j] = n_comp;
+                }
+              else if (indices_are_contiguous)
                 {
                   index_storage_variants[dof_access_cell][i] =
                     IndexStorageVariants::contiguous;
+                  for (unsigned int j = 0; j < n_comp; ++j)
+                    dof_indices_interleave_strides[2][i * vectorization_length +
+                                                      j] = 1;
                 }
               else
                 {
+                  int                 indices_are_interleaved_and_mixed = 2;
                   const unsigned int *dof_indices =
-                    this->dof_indices.data() +
-                    row_starts[i * vectorization_length * n_components].first;
-                  if (n_comp == vectorization_length)
-                    index_storage_variants[dof_access_cell][i] =
-                      IndexStorageVariants::interleaved;
-                  else
-                    index_storage_variants[dof_access_cell][i] =
-                      IndexStorageVariants::full;
-
-                  // do not use interleaved storage if two vectorized
-                  // components point to the same field (scatter not possible)
+                    &this->dof_indices[row_starts[i * vectorization_length *
+                                                  n_components]
+                                         .first];
+                  for (unsigned int j = 0; j < n_comp; ++j)
+                    offsets[j] =
+                      dof_indices[j * ndofs + 1] - dof_indices[j * ndofs];
                   for (unsigned int k = 0; k < ndofs; ++k)
-                    for (unsigned int l = 0; l < n_comp; ++l)
-                      for (unsigned int j = l + 1; j < n_comp; ++j)
-                        if (dof_indices[j * ndofs + k] ==
-                            dof_indices[l * ndofs + k])
+                    for (unsigned int j = 0; j < n_comp; ++j)
+                      if (dof_indices[j * ndofs + k] !=
+                          dof_indices[j * ndofs] + k * offsets[j])
+                        {
+                          indices_are_interleaved_and_mixed = 0;
+                          break;
+                        }
+                  if (indices_are_interleaved_and_mixed == 2)
+                    {
+                      for (unsigned int j = 0; j < n_comp; ++j)
+                        dof_indices_interleave_strides
+                          [dof_access_cell][i * vectorization_length + j] =
+                            offsets[j];
+                      for (unsigned int j = 0; j < n_comp; ++j)
+                        dof_indices_contiguous[dof_access_cell]
+                                              [i * vectorization_length + j] =
+                                                dof_indices[j * ndofs];
+                      for (unsigned int j = 0; j < n_comp; ++j)
+                        if (offsets[j] != vectorization_length)
                           {
-                            index_storage_variants[dof_access_cell][i] =
-                              IndexStorageVariants::full;
+                            indices_are_interleaved_and_mixed = 1;
                             break;
                           }
-                  if (index_storage_variants[dof_access_cell][i] !=
-                      IndexStorageVariants::full)
+                      if (indices_are_interleaved_and_mixed == 1 ||
+                          n_comp != vectorization_length)
+                        index_storage_variants[dof_access_cell][i] =
+                          IndexStorageVariants::
+                            interleaved_contiguous_mixed_strides;
+                      else
+                        index_storage_variants[dof_access_cell][i] =
+                          IndexStorageVariants::interleaved_contiguous_strided;
+                    }
+                  else
                     {
-                      unsigned int *interleaved_dof_indices =
-                        this->dof_indices_interleaved.data() +
+                      const unsigned int *dof_indices =
+                        this->dof_indices.data() +
                         row_starts[i * vectorization_length * n_components]
                           .first;
+                      if (n_comp == vectorization_length)
+                        index_storage_variants[dof_access_cell][i] =
+                          IndexStorageVariants::interleaved;
+                      else
+                        index_storage_variants[dof_access_cell][i] =
+                          IndexStorageVariants::full;
+
+                      // do not use interleaved storage if two vectorized
+                      // components point to the same field (scatter not
+                      // possible)
                       for (unsigned int k = 0; k < ndofs; ++k)
-                        for (unsigned int j = 0; j < n_comp; ++j)
-                          interleaved_dof_indices[k * n_comp + j] =
-                            dof_indices[j * ndofs + k];
+                        for (unsigned int l = 0; l < n_comp; ++l)
+                          for (unsigned int j = l + 1; j < n_comp; ++j)
+                            if (dof_indices[j * ndofs + k] ==
+                                dof_indices[l * ndofs + k])
+                              {
+                                index_storage_variants[dof_access_cell][i] =
+                                  IndexStorageVariants::full;
+                                break;
+                              }
                     }
                 }
             }
           index_kinds[static_cast<unsigned int>(
             index_storage_variants[dof_access_cell][i])]++;
         }
+
+      // Cleanup phase: we want to avoid single cells with different properties
+      // than the bulk of the domain in order to avoid extra checks in the face
+      // identification.
+
+      // Step 1: check whether the interleaved indices were only assigned to
+      // the single cell within a vectorized array.
+      auto fix_single_interleaved_indices =
+        [&](const IndexStorageVariants variant) {
+          if (index_kinds[static_cast<unsigned int>(
+                IndexStorageVariants::interleaved_contiguous_mixed_strides)] >
+                0 &&
+              index_kinds[static_cast<unsigned int>(variant)] > 0)
+            for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+              {
+                if (index_storage_variants[dof_access_cell][i] ==
+                      IndexStorageVariants::
+                        interleaved_contiguous_mixed_strides &&
+                    n_vectorization_lanes_filled[dof_access_cell][i] == 1 &&
+                    (variant != IndexStorageVariants::contiguous ||
+                     dof_indices_interleave_strides[dof_access_cell]
+                                                   [i * vectorization_length] ==
+                       1))
+                  {
+                    index_storage_variants[dof_access_cell][i] = variant;
+                    index_kinds[static_cast<unsigned int>(
+                      IndexStorageVariants::
+                        interleaved_contiguous_mixed_strides)]--;
+                    index_kinds[static_cast<unsigned int>(variant)]++;
+                  }
+              }
+        };
+
+      fix_single_interleaved_indices(IndexStorageVariants::full);
+      fix_single_interleaved_indices(IndexStorageVariants::contiguous);
+      fix_single_interleaved_indices(IndexStorageVariants::interleaved);
+
+      unsigned int n_interleaved =
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous)] +
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous_strided)] +
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous_mixed_strides)];
+
+      // Step 2: fix single contiguous cell among others with interleaved
+      // storage
+      if (n_interleaved > 0 && index_kinds[static_cast<unsigned int>(
+                                 IndexStorageVariants::contiguous)])
+        for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+          if (index_storage_variants[dof_access_cell][i] ==
+              IndexStorageVariants::contiguous)
+            {
+              index_storage_variants[dof_access_cell][i] =
+                IndexStorageVariants::interleaved_contiguous_mixed_strides;
+              index_kinds[static_cast<unsigned int>(
+                IndexStorageVariants::contiguous)]--;
+              index_kinds[static_cast<unsigned int>(
+                IndexStorageVariants::interleaved_contiguous_mixed_strides)]++;
+            }
+
+      n_interleaved =
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous)] +
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous_strided)] +
+        index_kinds[static_cast<unsigned int>(
+          IndexStorageVariants::interleaved_contiguous_mixed_strides)];
+
+      // Step 3: Interleaved cells are left but also some non-contiguous ones
+      // -> revert all to full storage
+      if (n_interleaved > 0 &&
+          index_kinds[static_cast<unsigned int>(IndexStorageVariants::full)] +
+              index_kinds[static_cast<unsigned int>(
+                IndexStorageVariants::interleaved)] >
+            0)
+        for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+          if (index_storage_variants[dof_access_cell][i] >
+              IndexStorageVariants::contiguous)
+            {
+              index_kinds[static_cast<unsigned int>(
+                index_storage_variants[2][i])]--;
+              if (n_vectorization_lanes_filled[dof_access_cell][i] ==
+                  vectorization_length)
+                index_storage_variants[dof_access_cell][i] =
+                  IndexStorageVariants::interleaved;
+              else
+                index_storage_variants[dof_access_cell][i] =
+                  IndexStorageVariants::full;
+              index_kinds[static_cast<unsigned int>(
+                index_storage_variants[dof_access_cell][i])]++;
+            }
+
+      // Step 4: Copy the interleaved indices into their own data structure
+      for (unsigned int i = 0; i < irregular_cells.size(); ++i)
+        if (index_storage_variants[dof_access_cell][i] ==
+            IndexStorageVariants::interleaved)
+          {
+            const unsigned int ndofs =
+              dofs_per_cell[have_hp ? cell_active_fe_index[i] : 0];
+            const unsigned int *dof_indices =
+              &this->dof_indices
+                 [row_starts[i * vectorization_length * n_components].first];
+            unsigned int *interleaved_dof_indices =
+              &this->dof_indices_interleaved
+                 [row_starts[i * vectorization_length * n_components].first];
+            for (unsigned int k = 0; k < ndofs; ++k)
+              for (unsigned int j = 0; j < vectorization_length; ++j)
+                interleaved_dof_indices[k * vectorization_length + j] =
+                  dof_indices[j * ndofs + k];
+          }
     }
 
 
@@ -833,6 +1006,8 @@ namespace internal
         faces.size(), IndexStorageVariants::full);
       dof_indices_contiguous[dof_access_face_interior].resize(
         faces.size() * length, numbers::invalid_unsigned_int);
+      dof_indices_interleave_strides[dof_access_face_interior].resize(
+        faces.size() * length, numbers::invalid_unsigned_int);
       n_vectorization_lanes_filled[dof_access_face_interior].resize(
         faces.size());
 
@@ -846,6 +1021,8 @@ namespace internal
         n_exterior_faces, IndexStorageVariants::full);
       dof_indices_contiguous[dof_access_face_exterior].resize(
         n_exterior_faces * length, numbers::invalid_unsigned_int);
+      dof_indices_interleave_strides[dof_access_face_exterior].resize(
+        faces.size() * length, numbers::invalid_unsigned_int);
       n_vectorization_lanes_filled[dof_access_face_exterior].resize(
         n_exterior_faces);
 
@@ -854,6 +1031,7 @@ namespace internal
           auto face_computation = [&](const DoFAccessIndex face_index,
                                       const unsigned int * cell_indices_face) {
             bool is_contiguous      = false;
+            bool is_interleaved     = false;
             bool needs_full_storage = false;
             for (unsigned int v = 0;
                  v < length &&
@@ -861,22 +1039,74 @@ namespace internal
                  ++v)
               {
                 n_vectorization_lanes_filled[face_index][face]++;
+                if (index_storage_variants[dof_access_cell]
+                                          [cell_indices_face[v] / length] >=
+                    IndexStorageVariants::interleaved_contiguous)
+                  is_interleaved = true;
                 if (index_storage_variants[dof_access_cell]
                                           [cell_indices_face[v] / length] ==
                     IndexStorageVariants::contiguous)
                   is_contiguous = true;
+                if (index_storage_variants[dof_access_cell]
+                                          [cell_indices_face[v] / length] >=
+                    IndexStorageVariants::contiguous)
+                  dof_indices_interleave_strides[face_index][face * length +
+                                                             v] =
+                    dof_indices_interleave_strides[dof_access_cell]
+                                                  [cell_indices_face[v]];
                 if (index_storage_variants[dof_access_cell]
                                           [cell_indices_face[v] / length] <
                     IndexStorageVariants::contiguous)
                   needs_full_storage = true;
               }
-            if (is_contiguous)
+            Assert(!(is_interleaved && is_contiguous),
+                   ExcMessage("Unsupported index compression found"));
+
+            if (is_interleaved || is_contiguous)
               for (unsigned int v = 0;
                    v < n_vectorization_lanes_filled[face_index][face];
                    ++v)
                 dof_indices_contiguous[face_index][face * length + v] =
                   dof_indices_contiguous[dof_access_cell][cell_indices_face[v]];
-            if (is_contiguous && !needs_full_storage)
+            if (is_interleaved)
+              {
+                bool is_also_contiguous =
+                  n_vectorization_lanes_filled[face_index][face] == length;
+                for (unsigned int v = 0;
+                     v < n_vectorization_lanes_filled[face_index][face];
+                     ++v)
+                  if (dof_indices_contiguous[face_index][face * length + v] !=
+                        dof_indices_contiguous[face_index][face * length] + v ||
+                      dof_indices_interleave_strides[dof_access_cell]
+                                                    [cell_indices_face[v]] !=
+                        length)
+                    is_also_contiguous = false;
+
+                if (is_also_contiguous)
+                  {
+                    index_storage_variants[face_index][face] =
+                      IndexStorageVariants::interleaved_contiguous;
+                  }
+                else
+                  {
+                    bool all_indices_same_offset =
+                      n_vectorization_lanes_filled[face_index][face] == length;
+                    for (unsigned int v = 0;
+                         v < n_vectorization_lanes_filled[face_index][face];
+                         ++v)
+                      if (dof_indices_interleave_strides
+                            [dof_access_cell][cell_indices_face[v]] != length)
+                        all_indices_same_offset = false;
+                    if (all_indices_same_offset)
+                      index_storage_variants[face_index][face] =
+                        IndexStorageVariants::interleaved_contiguous_strided;
+                    else
+                      index_storage_variants[face_index][face] =
+                        IndexStorageVariants::
+                          interleaved_contiguous_mixed_strides;
+                  }
+              }
+            else if (is_contiguous && !needs_full_storage)
               index_storage_variants[face_index][face] =
                 IndexStorageVariants::contiguous;
             else
index b8cd0248ccdb02987a5c7bfd9d45a7ca38826eb1..2fd382aa06478a2789915cd9400d08d35263d4ca 100644 (file)
@@ -86,7 +86,8 @@ namespace internal
         Number *   gradients_quad,
         Number *   scratch_data,
         const bool integrate_values,
-        const bool integrate_gradients)
+        const bool integrate_gradients,
+        const bool sum_into_values_array = false)
       {
         internal::FEEvaluationImpl<
           internal::MatrixFreeFunctions::tensor_general,
@@ -101,7 +102,7 @@ namespace internal
                              scratch_data,
                              integrate_values,
                              integrate_gradients,
-                             false);
+                             sum_into_values_array);
       }
     };
 
@@ -205,7 +206,8 @@ namespace internal
         Number *   gradients_quad,
         Number *   scratch_data,
         const bool integrate_values,
-        const bool integrate_gradients)
+        const bool integrate_gradients,
+        const bool sum_into_values_array = false)
       {
         const int              runtime_degree   = shape_info.fe_degree;
         constexpr unsigned int start_n_q_points = degree + 1;
@@ -217,7 +219,8 @@ namespace internal
                       gradients_quad,
                       scratch_data,
                       integrate_values,
-                      integrate_gradients);
+                      integrate_gradients,
+                      sum_into_values_array);
         else
           Factory<dim, n_components, Number, 0, degree + 1, n_q_points_1d>::
             integrate(shape_info,
@@ -226,7 +229,8 @@ namespace internal
                       gradients_quad,
                       scratch_data,
                       integrate_values,
-                      integrate_gradients);
+                      integrate_gradients,
+                      sum_into_values_array);
       }
     };
 
@@ -329,7 +333,8 @@ namespace internal
         Number *   gradients_quad,
         Number *   scratch_data,
         const bool integrate_values,
-        const bool integrate_gradients)
+        const bool integrate_gradients,
+        const bool sum_into_values_array)
       {
         const int runtime_n_q_points_1d = shape_info.n_q_points_1d;
         if (runtime_n_q_points_1d == n_q_points_1d)
@@ -346,7 +351,7 @@ namespace internal
                             scratch_data,
                             integrate_values,
                             integrate_gradients,
-                            false);
+                            sum_into_values_array);
             else if (degree < n_q_points_1d)
               internal::FEEvaluationImplTransformToCollocation<
                 dim,
@@ -360,7 +365,7 @@ namespace internal
                                    scratch_data,
                                    integrate_values,
                                    integrate_gradients,
-                                   false);
+                                   sum_into_values_array);
             else
               internal::FEEvaluationImpl<
                 internal::MatrixFreeFunctions::tensor_symmetric,
@@ -375,7 +380,7 @@ namespace internal
                                    scratch_data,
                                    integrate_values,
                                    integrate_gradients,
-                                   false);
+                                   sum_into_values_array);
           }
         else
           Factory<dim, n_components, Number, 1, degree, n_q_points_1d + 1>::
@@ -385,7 +390,8 @@ namespace internal
                       gradients_quad,
                       scratch_data,
                       integrate_values,
-                      integrate_gradients);
+                      integrate_gradients,
+                      sum_into_values_array);
       }
     };
 
@@ -437,7 +443,8 @@ namespace internal
       Number *   gradients_quad,
       Number *   scratch_data,
       const bool integrate_values,
-      const bool integrate_gradients)
+      const bool integrate_gradients,
+      const bool sum_into_values_array = false)
     {
       Assert(shape_info.element_type <=
                internal::MatrixFreeFunctions::tensor_symmetric,
@@ -448,7 +455,8 @@ namespace internal
                                                     gradients_quad,
                                                     scratch_data,
                                                     integrate_values,
-                                                    integrate_gradients);
+                                                    integrate_gradients,
+                                                    sum_into_values_array);
     }
   } // namespace EvaluationSelectorImplementation
 } // namespace internal
@@ -505,7 +513,8 @@ struct SelectEvaluator
             Number *   gradients_quad,
             Number *   scratch_data,
             const bool integrate_values,
-            const bool integrate_gradients);
+            const bool integrate_gradients,
+            const bool sum_into_values_array = false);
 };
 
 /**
@@ -555,7 +564,8 @@ struct SelectEvaluator<dim, -1, n_q_points_1d, n_components, Number>
             Number *   gradients_quad,
             Number *   scratch_data,
             const bool integrate_values,
-            const bool integrate_gradients);
+            const bool integrate_gradients,
+            const bool sum_into_values_array = false);
 };
 
 //----------------------Implementation for SelectEvaluator---------------------
@@ -709,7 +719,8 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
   Number *                                                gradients_quad,
   Number *                                                scratch_data,
   const bool                                              integrate_values,
-  const bool                                              integrate_gradients)
+  const bool                                              integrate_gradients,
+  const bool                                              sum_into_values_array)
 {
   Assert(fe_degree >= 0 && n_q_points_1d > 0, ExcInternalError());
 
@@ -726,7 +737,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                     scratch_data,
                     integrate_values,
                     integrate_gradients,
-                    false);
+                    sum_into_values_array);
     }
   else if (fe_degree < n_q_points_1d &&
            shape_info.element_type <=
@@ -744,7 +755,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::tensor_symmetric)
@@ -762,7 +773,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
@@ -780,7 +791,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::truncated_tensor)
@@ -798,7 +809,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::tensor_general)
@@ -815,7 +826,7 @@ SelectEvaluator<dim, fe_degree, n_q_points_1d, n_components, Number>::integrate(
                                                     scratch_data,
                                                     integrate_values,
                                                     integrate_gradients,
-                                                    false);
+                                                    sum_into_values_array);
     }
   else
     AssertThrow(false, ExcNotImplemented());
@@ -914,7 +925,8 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate(
   Number *                                                gradients_quad,
   Number *                                                scratch_data,
   const bool                                              integrate_values,
-  const bool                                              integrate_gradients)
+  const bool                                              integrate_gradients,
+  const bool                                              sum_into_values_array)
 {
   if (shape_info.element_type ==
       internal::MatrixFreeFunctions::tensor_symmetric_plus_dg0)
@@ -932,7 +944,7 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::truncated_tensor)
@@ -950,7 +962,7 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate(
                            scratch_data,
                            integrate_values,
                            integrate_gradients,
-                           false);
+                           sum_into_values_array);
     }
   else if (shape_info.element_type ==
            internal::MatrixFreeFunctions::tensor_general)
@@ -966,7 +978,7 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate(
                                                   scratch_data,
                                                   integrate_values,
                                                   integrate_gradients,
-                                                  false);
+                                                  sum_into_values_array);
   else
     internal::EvaluationSelectorImplementation::
       symmetric_selector_integrate<dim, n_components, Number>(
@@ -976,7 +988,8 @@ SelectEvaluator<dim, -1, dummy, n_components, Number>::integrate(
         gradients_quad,
         scratch_data,
         integrate_values,
-        integrate_gradients);
+        integrate_gradients,
+        sum_into_values_array);
 }
 #endif // DOXYGEN
 
index 1459d498d88f2766f36af2b54b79cf093b81e8f0..daafb5ad8f74082a8f84e32cc687a2e3e3f052f2 100644 (file)
@@ -3454,6 +3454,36 @@ namespace internal
       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,
+                            VectorizedArray<Number> *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 += VectorizedArray<Number>::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,
+                            VectorType &             vec,
+                            VectorizedArray<Number> *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+             ++v)
+          dof_values[i][v] = vector_access(
+            vec, dof_index + v + i * VectorizedArray<Number>::n_array_elements);
+    }
+
     template <typename VectorType>
     void
     process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
@@ -3560,6 +3590,42 @@ namespace internal
       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,
+                            VectorizedArray<Number> *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 += VectorizedArray<Number>::n_array_elements)
+        {
+          VectorizedArray<Number> 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,
+                            VectorizedArray<Number> *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+             ++v)
+          vector_access(vec,
+                        dof_index + v +
+                          i * VectorizedArray<Number>::n_array_elements) +=
+            dof_values[i][v];
+    }
+
     template <typename VectorType>
     void
     process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
@@ -3670,6 +3736,37 @@ namespace internal
       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,
+                            VectorizedArray<Number> *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 += VectorizedArray<Number>::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,
+                            VectorizedArray<Number> *dof_values,
+                            std::integral_constant<bool, false>) const
+    {
+      for (unsigned int i = 0; i < dofs_per_cell; ++i)
+        for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+             ++v)
+          vector_access(vec,
+                        dof_index + v +
+                          i * VectorizedArray<Number>::n_array_elements) =
+            dof_values[i][v];
+    }
+
     template <typename VectorType>
     void
     process_dofs_vectorized_transpose(const unsigned int       dofs_per_cell,
@@ -4294,6 +4391,37 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::
 
   const std::vector<unsigned int> &dof_indices_cont =
     dof_info->dof_indices_contiguous[ind];
+
+  // Simple case: We have contiguous storage, so we can simply copy out the
+  // data
+  if (dof_info->index_storage_variants[ind][cell] ==
+      internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+        interleaved_contiguous)
+    {
+      const unsigned int dof_index =
+        dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements] +
+        dof_info->component_dof_indices_offset[active_fe_index]
+                                              [first_selected_component] *
+          VectorizedArray<Number>::n_array_elements;
+      if (n_components == 1 || n_fe_components == 1)
+        for (unsigned int comp = 0; comp < n_components; ++comp)
+          operation.process_dofs_vectorized(data->dofs_per_component_on_cell,
+                                            dof_index,
+                                            *src[comp],
+                                            values_dofs[comp],
+                                            vector_selector);
+      else
+        operation.process_dofs_vectorized(data->dofs_per_component_on_cell *
+                                            n_components,
+                                          dof_index,
+                                          *src[0],
+                                          values_dofs[0],
+                                          vector_selector);
+      return;
+    }
+
+  // More general case: Must go through the components one by one and apply
+  // some transformations
   const unsigned int vectorization_populated =
     dof_info->n_vectorization_lanes_filled[ind][this->cell];
   unsigned int dof_indices[VectorizedArray<Number>::n_array_elements];
@@ -4301,7 +4429,9 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::
     dof_indices[v] =
       dof_indices_cont[cell * VectorizedArray<Number>::n_array_elements + v] +
       dof_info->component_dof_indices_offset[active_fe_index]
-                                            [first_selected_component];
+                                            [first_selected_component] *
+        dof_info->dof_indices_interleave_strides
+          [ind][cell * VectorizedArray<Number>::n_array_elements + v];
   for (unsigned int v = vectorization_populated;
        v < VectorizedArray<Number>::n_array_elements;
        ++v)
@@ -4311,40 +4441,149 @@ FEEvaluationBase<dim, n_components_, Number, is_face>::
   // constraints and that the indices within each element are contiguous
   if (vectorization_populated == VectorizedArray<Number>::n_array_elements)
     {
-      if (n_components == 1 || n_fe_components == 1)
-        for (unsigned int comp = 0; comp < n_components; ++comp)
-          operation.process_dofs_vectorized_transpose(
-            data->dofs_per_component_on_cell,
-            dof_indices,
-            *src[comp],
-            values_dofs[comp],
-            vector_selector);
+      if (dof_info->index_storage_variants[ind][cell] ==
+          internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+            contiguous)
+        {
+          if (n_components == 1 || n_fe_components == 1)
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              operation.process_dofs_vectorized_transpose(
+                data->dofs_per_component_on_cell,
+                dof_indices,
+                *src[comp],
+                values_dofs[comp],
+                vector_selector);
+          else
+            operation.process_dofs_vectorized_transpose(
+              data->dofs_per_component_on_cell * n_components,
+              dof_indices,
+              *src[0],
+              &values_dofs[0][0],
+              vector_selector);
+        }
+      else if (dof_info->index_storage_variants[ind][cell] ==
+               internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                 interleaved_contiguous_strided)
+        {
+          if (n_components == 1 || n_fe_components == 1)
+            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+              {
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  operation.process_dof_gather(
+                    dof_indices,
+                    *src[comp],
+                    i * VectorizedArray<Number>::n_array_elements,
+                    values_dofs[comp][i],
+                    vector_selector);
+              }
+          else
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                   ++i)
+                {
+                  operation.process_dof_gather(
+                    dof_indices,
+                    *src[0],
+                    (comp * data->dofs_per_component_on_cell + i) *
+                      VectorizedArray<Number>::n_array_elements,
+                    values_dofs[comp][i],
+                    vector_selector);
+                }
+        }
       else
-        operation.process_dofs_vectorized_transpose(
-          data->dofs_per_component_on_cell * n_components,
-          dof_indices,
-          *src[0],
-          &values_dofs[0][0],
-          vector_selector);
+        {
+          Assert(dof_info->index_storage_variants[ind][cell] ==
+                   internal::MatrixFreeFunctions::DoFInfo::
+                     IndexStorageVariants::interleaved_contiguous_mixed_strides,
+                 ExcNotImplemented());
+          const unsigned int *offsets =
+            &dof_info->dof_indices_interleave_strides
+               [ind][VectorizedArray<Number>::n_array_elements * cell];
+          if (n_components == 1 || n_fe_components == 1)
+            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
+              {
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  operation.process_dof_gather(dof_indices,
+                                               *src[comp],
+                                               0,
+                                               values_dofs[comp][i],
+                                               vector_selector);
+                DEAL_II_OPENMP_SIMD_PRAGMA
+                for (unsigned int v = 0;
+                     v < VectorizedArray<Number>::n_array_elements;
+                     ++v)
+                  dof_indices[v] += offsets[v];
+              }
+          else
+            for (unsigned int comp = 0; comp < n_components; ++comp)
+              for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                   ++i)
+                {
+                  operation.process_dof_gather(dof_indices,
+                                               *src[0],
+                                               0,
+                                               values_dofs[comp][i],
+                                               vector_selector);
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::n_array_elements;
+                       ++v)
+                    dof_indices[v] += offsets[v];
+                }
+        }
     }
   else
     for (unsigned int comp = 0; comp < n_components; ++comp)
       {
         for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
           operation.process_empty(values_dofs[comp][i]);
-        if (n_components == 1 || n_fe_components == 1)
-          for (unsigned int v = 0; v < vectorization_populated; ++v)
-            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
-              operation.process_dof(dof_indices[v] + i,
-                                    *src[comp],
-                                    values_dofs[comp][i][v]);
+        if (dof_info->index_storage_variants[ind][cell] ==
+            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+              contiguous)
+          {
+            if (n_components == 1 || n_fe_components == 1)
+              for (unsigned int v = 0; v < vectorization_populated; ++v)
+                for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                     ++i)
+                  operation.process_dof(dof_indices[v] + i,
+                                        *src[comp],
+                                        values_dofs[comp][i][v]);
+            else
+              for (unsigned int v = 0; v < vectorization_populated; ++v)
+                for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                     ++i)
+                  operation.process_dof(dof_indices[v] + i +
+                                          comp *
+                                            data->dofs_per_component_on_cell,
+                                        *src[0],
+                                        values_dofs[comp][i][v]);
+          }
         else
-          for (unsigned int v = 0; v < vectorization_populated; ++v)
-            for (unsigned int i = 0; i < data->dofs_per_component_on_cell; ++i)
-              operation.process_dof(dof_indices[v] + i +
-                                      comp * data->dofs_per_component_on_cell,
-                                    *src[0],
-                                    values_dofs[comp][i][v]);
+          {
+            const unsigned int *offsets =
+              &dof_info->dof_indices_interleave_strides
+                 [ind][VectorizedArray<Number>::n_array_elements * cell];
+            for (unsigned int v = 0; v < vectorization_populated; ++v)
+              AssertIndexRange(offsets[v],
+                               VectorizedArray<Number>::n_array_elements + 1);
+            if (n_components == 1 || n_fe_components == 1)
+              for (unsigned int v = 0; v < vectorization_populated; ++v)
+                for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                     ++i)
+                  operation.process_dof(dof_indices[v] + i * offsets[v],
+                                        *src[comp],
+                                        values_dofs[comp][i][v]);
+            else
+              for (unsigned int v = 0; v < vectorization_populated; ++v)
+                for (unsigned int i = 0; i < data->dofs_per_component_on_cell;
+                     ++i)
+                  operation.process_dof(
+                    dof_indices[v] +
+                      (i + comp * data->dofs_per_component_on_cell) *
+                        offsets[v],
+                    *src[0],
+                    values_dofs[comp][i][v]);
+          }
       }
 }
 
@@ -6646,11 +6885,42 @@ FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                   const bool        evaluate_gradients,
                   const bool        evaluate_hessians)
 {
-  this->read_dof_values(input_vector);
-  evaluate(this->begin_dof_values(),
-           evaluate_values,
-           evaluate_gradients,
-           evaluate_hessians);
+  if (this->dof_info->index_storage_variants
+          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
+                            IndexStorageVariants::interleaved_contiguous &&
+      reinterpret_cast<std::size_t>(
+        input_vector.begin() +
+        this->dof_info->dof_indices_contiguous
+          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+          [this->cell * VectorizedArray<Number>::n_array_elements]) %
+          sizeof(VectorizedArray<Number>) ==
+        0)
+    {
+      const VectorizedArray<Number> *vec_values =
+        reinterpret_cast<const VectorizedArray<Number> *>(
+          input_vector.begin() +
+          this->dof_info->dof_indices_contiguous
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+            [this->cell * VectorizedArray<Number>::n_array_elements] +
+          this->dof_info
+              ->component_dof_indices_offset[this->active_fe_index]
+                                            [this->first_selected_component] *
+            VectorizedArray<Number>::n_array_elements);
+
+      evaluate(vec_values,
+               evaluate_values,
+               evaluate_gradients,
+               evaluate_hessians);
+    }
+  else
+    {
+      this->read_dof_values(input_vector);
+      evaluate(this->begin_dof_values(),
+               evaluate_values,
+               evaluate_gradients,
+               evaluate_hessians);
+    }
 }
 
 
@@ -6706,7 +6976,8 @@ FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::integrate(
                                                         ->gradients_quad[0][0],
                                                       this->scratch_data,
                                                       integrate_values,
-                                                      integrate_gradients);
+                                                      integrate_gradients,
+                                                      false);
 
 #  ifdef DEBUG
   this->dof_values_initialized = true;
@@ -6727,8 +6998,49 @@ FEEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                     const bool  integrate_gradients,
                     VectorType &destination)
 {
-  integrate(integrate_values, integrate_gradients, this->begin_dof_values());
-  this->distribute_local_to_global(destination);
+  if (this->dof_info->index_storage_variants
+          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+          [this->cell] == internal::MatrixFreeFunctions::DoFInfo::
+                            IndexStorageVariants::interleaved_contiguous &&
+      reinterpret_cast<std::size_t>(
+        destination.begin() +
+        this->dof_info->dof_indices_contiguous
+          [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+          [this->cell * VectorizedArray<Number>::n_array_elements]) %
+          sizeof(VectorizedArray<Number>) ==
+        0)
+    {
+      VectorizedArray<Number> *vec_values =
+        reinterpret_cast<VectorizedArray<Number> *>(
+          destination.begin() +
+          this->dof_info->dof_indices_contiguous
+            [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+            [this->cell * VectorizedArray<Number>::n_array_elements] +
+          this->dof_info
+              ->component_dof_indices_offset[this->active_fe_index]
+                                            [this->first_selected_component] *
+            VectorizedArray<Number>::n_array_elements);
+      SelectEvaluator<
+        dim,
+        fe_degree,
+        n_q_points_1d,
+        n_components,
+        VectorizedArray<Number>>::integrate(*this->data,
+                                            vec_values,
+                                            this->values_quad[0],
+                                            this->gradients_quad[0][0],
+                                            this->scratch_data,
+                                            integrate_values,
+                                            integrate_gradients,
+                                            true);
+    }
+  else
+    {
+      integrate(integrate_values,
+                integrate_gradients,
+                this->begin_dof_values());
+      this->distribute_local_to_global(destination);
+    }
 }
 
 
@@ -7167,19 +7479,357 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
     temp1 = this->scratch_data;
 
   internal::VectorReader<Number> reader;
+  std::integral_constant<
+    bool,
+    std::is_same<typename VectorType::value_type, Number>::value>
+    vector_selector;
 
-  if (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] ==
-        VectorizedArray<Number>::n_array_elements &&
-      ((evaluate_gradients == false &&
+  // 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)))
+        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],
+        VectorizedArray<Number>::n_array_elements);
+      const unsigned int dof_index =
+        this->dof_info
+          ->dof_indices_contiguous[this->dof_access_index]
+                                  [this->cell *
+                                   VectorizedArray<Number>::n_array_elements] +
+        this->dof_info
+            ->component_dof_indices_offset[this->active_fe_index]
+                                          [this->first_selected_component] *
+          VectorizedArray<Number>::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 VectorizedArray<Number> 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) *
+                                  VectorizedArray<Number>::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) *
+                                  VectorizedArray<Number>::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) *
+                                VectorizedArray<Number>::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],
+        VectorizedArray<Number>::n_array_elements);
+      const unsigned int *indices =
+        &this->dof_info
+           ->dof_indices_contiguous[this->dof_access_index]
+                                   [this->cell *
+                                    VectorizedArray<Number>::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 VectorizedArray<Number> 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] * VectorizedArray<Number>::n_array_elements;
+              const unsigned int ind2 =
+                index_array[2 * i + 1] *
+                VectorizedArray<Number>::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 *
+                        VectorizedArray<Number>::n_array_elements +
+                      this->dof_info->component_dof_indices_offset
+                          [this->active_fe_index]
+                          [this->first_selected_component] *
+                        VectorizedArray<Number>::n_array_elements,
+                    temp1[i + 2 * comp * dofs_per_face],
+                    vector_selector);
+                  VectorizedArray<Number> grad;
+                  reader.process_dof_gather(
+                    indices,
+                    input_vector,
+                    ind2 +
+                      comp * static_dofs_per_component *
+                        VectorizedArray<Number>::n_array_elements +
+                      this->dof_info->component_dof_indices_offset
+                          [this->active_fe_index]
+                          [this->first_selected_component] *
+                        VectorizedArray<Number>::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] * VectorizedArray<Number>::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 *
+                      VectorizedArray<Number>::n_array_elements +
+                    this->dof_info->component_dof_indices_offset
+                        [this->active_fe_index]
+                        [this->first_selected_component] *
+                      VectorizedArray<Number>::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 * VectorizedArray<Number>::n_array_elements];
+      unsigned int indices[VectorizedArray<Number>::n_array_elements];
+      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+           ++v)
+        indices[v] =
+          this->dof_info->dof_indices_contiguous
+            [this->dof_access_index]
+            [this->cell * VectorizedArray<Number>::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 VectorizedArray<Number> 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 == VectorizedArray<Number>::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[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::n_array_elements;
+                       ++v)
+                    ind1[v] = indices[v] + (comp * static_dofs_per_component +
+                                            index_array[2 * i]) *
+                                             strides[v];
+                  unsigned int ind2[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::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);
+                  VectorizedArray<Number> 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] = VectorizedArray<Number>();
+              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 == VectorizedArray<Number>::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[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::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] = VectorizedArray<Number>();
+              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] ==
+             VectorizedArray<Number>::n_array_elements)
     {
       const unsigned int *indices =
         &this->dof_info
@@ -7193,12 +7843,8 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
           // 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 VectorizedArray<Number> grad_weight0 =
-            (side ? -1. : 1.) *
-            this->data->shape_data_on_face[0][fe_degree + 1];
-          const VectorizedArray<Number> grad_weight1 =
-            (side ? -1. : 1.) *
-            this->data->shape_data_on_face[0][fe_degree + 2];
+          const VectorizedArray<Number> 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);
 
@@ -7217,10 +7863,7 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                       this->dof_info->component_dof_indices_offset
                         [this->active_fe_index][this->first_selected_component],
                     temp1[i + 2 * comp * dofs_per_face],
-                    std::integral_constant<
-                      bool,
-                      std::is_same<typename VectorType::value_type,
-                                   Number>::value>());
+                    vector_selector);
                   VectorizedArray<Number> grad;
                   reader.process_dof_gather(
                     indices,
@@ -7229,13 +7872,9 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                       this->dof_info->component_dof_indices_offset
                         [this->active_fe_index][this->first_selected_component],
                     grad,
-                    std::integral_constant<
-                      bool,
-                      std::is_same<typename VectorType::value_type,
-                                   Number>::value>());
+                    vector_selector);
                   temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                    grad_weight0 * temp1[i + 2 * comp * dofs_per_face] +
-                    grad_weight1 * grad;
+                    grad_weight * (temp1[i + 2 * comp * dofs_per_face] - grad);
                 }
             }
         }
@@ -7256,13 +7895,12 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                     this->dof_info->component_dof_indices_offset
                       [this->active_fe_index][this->first_selected_component],
                   temp1[i + comp * 2 * dofs_per_face],
-                  std::integral_constant<
-                    bool,
-                    std::is_same<typename VectorType::value_type,
-                                 Number>::value>());
+                  vector_selector);
               }
         }
     }
+
+  // case 5: default vector access
   else
     {
       this->read_dof_values(input_vector);
@@ -7401,19 +8039,358 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
 #  endif
 
   internal::VectorDistributorLocalToGlobal<Number> writer;
+  std::integral_constant<
+    bool,
+    std::is_same<typename VectorType::value_type, Number>::value>
+    vector_selector;
 
-  if (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] ==
-        VectorizedArray<Number>::n_array_elements &&
-      ((integrate_gradients == false &&
+  // 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)))
+        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],
+        VectorizedArray<Number>::n_array_elements);
+      const unsigned int dof_index =
+        this->dof_info
+          ->dof_indices_contiguous[this->dof_access_index]
+                                  [this->cell *
+                                   VectorizedArray<Number>::n_array_elements] +
+        this->dof_info
+            ->component_dof_indices_offset[this->active_fe_index]
+                                          [this->first_selected_component] *
+          VectorizedArray<Number>::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 VectorizedArray<Number> 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)
+                {
+                  VectorizedArray<Number> val =
+                    temp1[i + 2 * comp * dofs_per_face] -
+                    grad_weight *
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                  VectorizedArray<Number> 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) *
+                                  VectorizedArray<Number>::n_array_elements,
+                    destination,
+                    &val,
+                    vector_selector);
+                  writer.process_dofs_vectorized(
+                    1,
+                    dof_index + (ind2 + comp * static_dofs_per_component) *
+                                  VectorizedArray<Number>::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) *
+                                VectorizedArray<Number>::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],
+        VectorizedArray<Number>::n_array_elements);
+      const unsigned int *indices =
+        &this->dof_info
+           ->dof_indices_contiguous[this->dof_access_index]
+                                   [this->cell *
+                                    VectorizedArray<Number>::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 VectorizedArray<Number> 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] * VectorizedArray<Number>::n_array_elements;
+              const unsigned int ind2 =
+                index_array[2 * i + 1] *
+                VectorizedArray<Number>::n_array_elements;
+              for (unsigned int comp = 0; comp < n_components_; ++comp)
+                {
+                  VectorizedArray<Number> val =
+                    temp1[i + 2 * comp * dofs_per_face] -
+                    grad_weight *
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                  VectorizedArray<Number> 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 *
+                        VectorizedArray<Number>::n_array_elements +
+                      this->dof_info->component_dof_indices_offset
+                          [this->active_fe_index]
+                          [this->first_selected_component] *
+                        VectorizedArray<Number>::n_array_elements,
+                    val,
+                    vector_selector);
+                  writer.process_dof_gather(
+                    indices,
+                    destination,
+                    ind2 +
+                      comp * static_dofs_per_component *
+                        VectorizedArray<Number>::n_array_elements +
+                      this->dof_info->component_dof_indices_offset
+                          [this->active_fe_index]
+                          [this->first_selected_component] *
+                        VectorizedArray<Number>::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] * VectorizedArray<Number>::n_array_elements;
+              for (unsigned int comp = 0; comp < n_components_; ++comp)
+                writer.process_dof_gather(
+                  indices,
+                  destination,
+                  ind +
+                    comp * static_dofs_per_component *
+                      VectorizedArray<Number>::n_array_elements +
+                    this->dof_info->component_dof_indices_offset
+                        [this->active_fe_index]
+                        [this->first_selected_component] *
+                      VectorizedArray<Number>::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 * VectorizedArray<Number>::n_array_elements];
+      unsigned int indices[VectorizedArray<Number>::n_array_elements];
+      for (unsigned int v = 0; v < VectorizedArray<Number>::n_array_elements;
+           ++v)
+        indices[v] =
+          this->dof_info->dof_indices_contiguous
+            [this->dof_access_index]
+            [this->cell * VectorizedArray<Number>::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 VectorizedArray<Number> 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 == VectorizedArray<Number>::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[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::n_array_elements;
+                       ++v)
+                    ind1[v] = indices[v] + (comp * static_dofs_per_component +
+                                            index_array[2 * i]) *
+                                             strides[v];
+                  unsigned int ind2[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::n_array_elements;
+                       ++v)
+                    ind2[v] = indices[v] + (comp * static_dofs_per_component +
+                                            index_array[2 * i + 1]) *
+                                             strides[v];
+                  VectorizedArray<Number> val =
+                    temp1[i + 2 * comp * dofs_per_face] -
+                    grad_weight *
+                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
+                  VectorizedArray<Number> 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 == VectorizedArray<Number>::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[VectorizedArray<Number>::n_array_elements];
+                  DEAL_II_OPENMP_SIMD_PRAGMA
+                  for (unsigned int v = 0;
+                       v < VectorizedArray<Number>::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] ==
+             VectorizedArray<Number>::n_array_elements)
     {
       const unsigned int *indices =
         &this->dof_info
@@ -7428,12 +8405,8 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
           // 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 VectorizedArray<Number> grad_weight0 =
-            (side ? -1. : 1.) *
-            this->data->shape_data_on_face[0][fe_degree + 1];
-          const VectorizedArray<Number> grad_weight1 =
-            (side ? -1. : 1.) *
-            this->data->shape_data_on_face[0][fe_degree + 2];
+          const VectorizedArray<Number> 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 =
@@ -7445,11 +8418,11 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
               for (unsigned int comp = 0; comp < n_components_; ++comp)
                 {
                   VectorizedArray<Number> val =
-                    temp1[i + 2 * comp * dofs_per_face] +
-                    grad_weight0 *
+                    temp1[i + 2 * comp * dofs_per_face] -
+                    grad_weight *
                       temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
                   VectorizedArray<Number> grad =
-                    grad_weight1 *
+                    grad_weight *
                     temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
                   writer.process_dof_gather(
                     indices,
@@ -7458,10 +8431,7 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                       this->dof_info->component_dof_indices_offset
                         [this->active_fe_index][this->first_selected_component],
                     val,
-                    std::integral_constant<
-                      bool,
-                      std::is_same<typename VectorType::value_type,
-                                   Number>::value>());
+                    vector_selector);
                   writer.process_dof_gather(
                     indices,
                     destination,
@@ -7469,10 +8439,7 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                       this->dof_info->component_dof_indices_offset
                         [this->active_fe_index][this->first_selected_component],
                     grad,
-                    std::integral_constant<
-                      bool,
-                      std::is_same<typename VectorType::value_type,
-                                   Number>::value>());
+                    vector_selector);
                 }
             }
         }
@@ -7493,13 +8460,12 @@ FEFaceEvaluation<dim, fe_degree, n_q_points_1d, n_components_, Number>::
                     this->dof_info->component_dof_indices_offset
                       [this->active_fe_index][this->first_selected_component],
                   temp1[i + 2 * comp * dofs_per_face],
-                  std::integral_constant<
-                    bool,
-                    std::is_same<typename VectorType::value_type,
-                                 Number>::value>());
+                  vector_selector);
             }
         }
     }
+
+  // case 5: default vector access
   else
     {
       internal::FEFaceNormalEvaluationImpl<dim,
index cf6437f9c0bf2335fb5d34f7c7d5a01c7115d14b..ff5e25152d65e8715225cfc662c418c424483c0f 100644 (file)
@@ -1542,8 +1542,9 @@ MatrixFree<dim, Number>::initialize_indices(
                         {
                           const unsigned int p =
                             face_info.faces[f].cells_exterior[v];
-                          const unsigned int stride = 1;
-                          unsigned int       i      = 0;
+                          const unsigned int stride =
+                            dof_info[no].dof_indices_interleave_strides[2][p];
+                          unsigned int i = 0;
                           for (unsigned int e = 0;
                                e < dof_info[no].n_base_elements;
                                ++e)
@@ -1657,8 +1658,9 @@ MatrixFree<dim, Number>::initialize_indices(
                         {
                           const unsigned int p =
                             face_info.faces[f].cells_exterior[v];
-                          const unsigned int stride = 1;
-                          unsigned int       i      = 0;
+                          const unsigned int stride =
+                            dof_info[no].dof_indices_interleave_strides[2][p];
+                          unsigned int i = 0;
                           for (unsigned int e = 0;
                                e < dof_info[no].n_base_elements;
                                ++e)
index 5be002385dcf17aea77e66ea313d29579555a433..e456056bb490c90494496b80f2dea49ca7adb4fc 100644 (file)
@@ -29,6 +29,7 @@ for (deal_II_dimension : DIMENSIONS; components : SPACE_DIMENSIONS;
                 VectorizedArray<scalar_type> *,
                 VectorizedArray<scalar_type> *,
                 const bool,
+                const bool,
                 const bool);
 
     template void SelectEvaluator<deal_II_dimension,

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.