]> https://gitweb.dealii.org/ - dealii.git/commitdiff
ECL: enable contiguous data access 10821/head
authorPeter Munch <peterrmuench@gmail.com>
Wed, 12 Aug 2020 08:14:21 +0000 (10:14 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Thu, 13 Aug 2020 15:59:15 +0000 (17:59 +0200)
include/deal.II/matrix_free/evaluation_kernels.h
include/deal.II/matrix_free/fe_evaluation.h
tests/matrix_free/ecl_02.cc [new file with mode: 0644]
tests/matrix_free/ecl_02.output [new file with mode: 0644]

index 8c3ad32cac4fda0a74d52791fdc790dc054b31ea..1a06ec4a525022b98bc36a6a787400873b71f1f2 100644 (file)
@@ -2069,11 +2069,11 @@ namespace internal
       const bool         evaluate_gradients,
       const unsigned int active_fe_index,
       const unsigned int first_selected_component,
-      const unsigned int cell,
-      const std::array<unsigned int, n_face_orientations> face_no,
+      const std::array<unsigned int, n_face_orientations> cells,
+      const std::array<unsigned int, n_face_orientations> face_nos,
       const unsigned int                                  subface_index,
       const MatrixFreeFunctions::DoFInfo::DoFAccessIndex  dof_access_index,
-      const std::array<unsigned int, n_face_orientations> face_orientation,
+      const std::array<unsigned int, n_face_orientations> face_orientations,
       const Table<2, unsigned int> &                      orientation_map)
     {
       return process_and_io( //
@@ -2088,11 +2088,11 @@ namespace internal
         evaluate_gradients,
         active_fe_index,
         first_selected_component,
-        cell,
-        face_no,
+        cells,
+        face_nos,
         subface_index,
         dof_access_index,
-        face_orientation,
+        face_orientations,
         orientation_map,
         [](auto &      temp_1,
            auto &      temp_2,
@@ -2194,11 +2194,11 @@ namespace internal
       const bool         integrate_gradients,
       const unsigned int active_fe_index,
       const unsigned int first_selected_component,
-      const unsigned int cell,
-      const std::array<unsigned int, n_face_orientations> face_no,
+      const std::array<unsigned int, n_face_orientations> cells,
+      const std::array<unsigned int, n_face_orientations> face_nos,
       const unsigned int                                  subface_index,
       const MatrixFreeFunctions::DoFInfo::DoFAccessIndex  dof_access_index,
-      const std::array<unsigned int, n_face_orientations> face_orientation,
+      const std::array<unsigned int, n_face_orientations> face_orientations,
       const Table<2, unsigned int> &                      orientation_map)
     {
       return process_and_io( //
@@ -2213,11 +2213,11 @@ namespace internal
         integrate_gradients,
         active_fe_index,
         first_selected_component,
-        cell,
-        face_no,
+        cells,
+        face_nos,
         subface_index,
         dof_access_index,
-        face_orientation,
+        face_orientations,
         orientation_map,
         [](const auto &temp_1,
            const auto &temp_2,
@@ -2270,14 +2270,14 @@ namespace internal
           // case 5: default vector access, must be handled separately, just do
           // the face-normal interpolation
 
-          AssertDimension(face_no.size(), 1);
+          AssertDimension(face_nos.size(), 1);
 
           FEFaceNormalEvaluationImpl<dim,
                                      fe_degree,
                                      n_components,
                                      VectorizedArrayType>::
             template interpolate<false, false>(
-              data, temp1, values_array, integrate_gradients, face_no[0]);
+              data, temp1, values_array, integrate_gradients, face_nos[0]);
         },
         [&](auto &temp1, const auto &dofs_per_face) {
           if (fe_degree > -1 &&
@@ -2344,11 +2344,11 @@ namespace internal
       const bool                                                 do_gradients,
       const unsigned int active_fe_index,
       const unsigned int first_selected_component,
-      const unsigned int cell,
-      const std::array<unsigned int, n_face_orientations> face_no,
+      const std::array<unsigned int, n_face_orientations> cells,
+      const std::array<unsigned int, n_face_orientations> face_nos,
       const unsigned int                                  subface_index,
       const MatrixFreeFunctions::DoFInfo::DoFAccessIndex  dof_access_index,
-      const std::array<unsigned int, n_face_orientations> face_orientation,
+      const std::array<unsigned int, n_face_orientations> face_orientations,
       const Table<2, unsigned int> &                      orientation_map,
       const Function1a &                                  function_1a,
       const Function1b &                                  function_1b,
@@ -2359,14 +2359,14 @@ namespace internal
       const Function5 &                                   function_5,
       const Function0 &                                   function_0)
     {
-      (void)subface_index;
+      const unsigned int cell = cells[0];
 
       if (integrate &&
-          (face_orientation[0] > 0 &&
+          (face_orientations[0] > 0 &&
            subface_index < GeometryInfo<dim>::max_children_per_cell))
         {
-          AssertDimension(face_orientation.size(), 1);
-          adjust_for_face_orientation(face_orientation[0],
+          AssertDimension(face_orientations.size(), 1);
+          adjust_for_face_orientation(face_orientations[0],
                                       orientation_map,
                                       true,
                                       do_values,
@@ -2385,8 +2385,8 @@ namespace internal
          data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
           data.data.front()
             .shape_data_on_face[0][fe_degree + (integrate ?
-                                                  (2 - (face_no[0] % 2)) :
-                                                  (1 + (face_no[0] % 2)))] :
+                                                  (2 - (face_nos[0] % 2)) :
+                                                  (1 + (face_nos[0] % 2)))] :
           VectorizedArrayType(0.0 /*dummy*/);
 
       constexpr unsigned int static_dofs_per_component =
@@ -2419,388 +2419,340 @@ namespace internal
 
       // re-orientation
       std::array<const unsigned int *, n_face_orientations> orientation;
-      orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
-                         &data.face_orientations[face_orientation[0]][0] :
-                         &dummy;
+      if (n_face_orientations == 1)
+        orientation[0] = (data.data.front().nodal_at_cell_boundaries == true) ?
+                           &data.face_orientations[face_orientations[0]][0] :
+                           &dummy;
+      else
+        {
+          for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+            {
+              // the loop breaks once an invalid_unsigned_int is hit for
+              // all cases except the exterior faces in the ECL loop (where
+              // some faces might be at the boundaries but others not)
+              if (cells[v] == numbers::invalid_unsigned_int)
+                continue;
+
+              orientation[v] =
+                (data.data.front().nodal_at_cell_boundaries == true) ?
+                  &data.face_orientations[face_orientations[v]][0] :
+                  &dummy;
+            }
+        }
 
       // face_to_cell_index_hermite
       std::array<const unsigned int *, n_face_orientations> index_array_hermite;
 
-      index_array_hermite[0] =
-        (data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
-         data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
-          &data.face_to_cell_index_hermite(face_no[0], 0) :
-          &dummy;
+      if (n_face_orientations == 1)
+        index_array_hermite[0] =
+          (data.data.front().nodal_at_cell_boundaries == true &&
+           fe_degree > 1 &&
+           data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite) ?
+            &data.face_to_cell_index_hermite(face_nos[0], 0) :
+            &dummy;
 
       if (n_face_orientations > 1 &&
           data.data.front().nodal_at_cell_boundaries == true && fe_degree > 1 &&
           data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite)
         {
-          const unsigned int n_filled_lanes =
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-          for (unsigned int v = 1; v < n_filled_lanes; ++v)
+          for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
             {
+              if (cells[v] == numbers::invalid_unsigned_int)
+                continue;
+
               grad_weight[v] =
                 data.data.front().shape_data_on_face
-                  [0][fe_degree + (integrate ? (2 - (face_no[v] % 2)) :
-                                               (1 + (face_no[v] % 2)))][v];
+                  [0][fe_degree + (integrate ? (2 - (face_nos[v] % 2)) :
+                                               (1 + (face_nos[v] % 2)))][v];
 
               index_array_hermite[v] =
-                &data.face_to_cell_index_hermite(face_no[v], 0);
+                &data.face_to_cell_index_hermite(face_nos[v], 0);
             }
         }
 
       // face_to_cell_index_nodal
       std::array<const unsigned int *, n_face_orientations> index_array_nodal;
 
-      index_array_nodal[0] =
-        (data.data.front().nodal_at_cell_boundaries == true) ?
-          &data.face_to_cell_index_nodal(face_no[0], 0) :
-          &dummy;
+      if (n_face_orientations == 1)
+        index_array_nodal[0] =
+          (data.data.front().nodal_at_cell_boundaries == true) ?
+            &data.face_to_cell_index_nodal(face_nos[0], 0) :
+            &dummy;
 
       if (n_face_orientations > 1 &&
           (data.data.front().nodal_at_cell_boundaries == true))
         {
-          const unsigned int n_filled_lanes =
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
+          for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+            {
+              if (cells[v] == numbers::invalid_unsigned_int)
+                continue;
 
-          for (unsigned int v = 1; v < n_filled_lanes; ++v)
-            index_array_nodal[v] =
-              &data.face_to_cell_index_nodal(face_no[v], 0);
+              index_array_nodal[v] =
+                &data.face_to_cell_index_nodal(face_nos[v], 0);
+            }
         }
 
       const auto reorientate = [&](const unsigned int v, const unsigned int i) {
-        return (dim < 3 || face_orientation[0] == 0 ||
+        return (dim < 3 ||
+                face_orientations[n_face_orientations == 1 ? 0 : v] == 0 ||
                 subface_index < GeometryInfo<dim>::max_children_per_cell) ?
                  i :
                  orientation[v][i];
       };
 
-      // case 1: contiguous and interleaved indices
-      if (((do_gradients == false &&
-            data.data.front().nodal_at_cell_boundaries == true) ||
-           (data.element_type ==
-              MatrixFreeFunctions::tensor_symmetric_hermite &&
-            fe_degree > 1)) &&
-          dof_info.index_storage_variants[dof_access_index][cell] ==
-            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-              interleaved_contiguous)
+      if ((do_gradients == false &&
+           data.data.front().nodal_at_cell_boundaries == true) ||
+          (data.element_type == MatrixFreeFunctions::tensor_symmetric_hermite &&
+           fe_degree > 1))
         {
-          AssertDimension(
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-            VectorizedArrayType::size());
-          const unsigned int dof_index =
-            dof_info
-              .dof_indices_contiguous[dof_access_index]
-                                     [cell * VectorizedArrayType::size()] +
-            dof_info.component_dof_indices_offset[active_fe_index]
+          // case 1: contiguous and interleaved indices
+          if (n_face_orientations == 1 &&
+              dof_info.index_storage_variants[dof_access_index][cell] ==
+                MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                  interleaved_contiguous)
+            {
+              AssertDimension(n_face_orientations, 1);
+
+              AssertDimension(
+                dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+                VectorizedArrayType::size());
+              const unsigned int dof_index =
+                dof_info
+                  .dof_indices_contiguous[dof_access_index]
+                                         [cell * VectorizedArrayType::size()] +
+                dof_info
+                    .component_dof_indices_offset[active_fe_index]
                                                  [first_selected_component] *
-              VectorizedArrayType::size();
+                  VectorizedArrayType::size();
 
-          if (fe_degree > 1 && do_gradients == true)
-            {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
+              if (fe_degree > 1 && do_gradients == true)
                 {
-                  if (n_face_orientations == 1)
-                    {
-                      const unsigned int ind1 = index_array_hermite[0][2 * i];
-                      const unsigned int ind2 =
-                        index_array_hermite[0][2 * i + 1];
-                      AssertIndexRange(ind1, data.dofs_per_component_on_cell);
-                      AssertIndexRange(ind2, data.dofs_per_component_on_cell);
-                      const unsigned int i_ = reorientate(0, i);
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        function_1a(
-                          temp1[i_ + 2 * comp * dofs_per_face],
-                          temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
-                          global_vector_ptr + dof_index +
-                            (ind1 + comp * static_dofs_per_component) *
-                              VectorizedArrayType::size(),
-                          global_vector_ptr + dof_index +
-                            (ind2 + comp * static_dofs_per_component) *
-                              VectorizedArrayType::size(),
-                          grad_weight);
-                    }
-                  else
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
-                      Assert(false, ExcNotImplemented());
-
-                      const unsigned int n_filled_lanes =
-                        dof_info
-                          .n_vectorization_lanes_filled[dof_access_index][cell];
-
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          function_3a(
-                            temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
-                                 [v],
-                            temp1[reorientate(v, i) + dofs_per_face +
-                                  2 * comp * dofs_per_face][v],
-                            global_vector_ptr[dof_index +
-                                              (index_array_hermite[v][2 * i] +
-                                               comp *
-                                                 static_dofs_per_component) *
-                                                VectorizedArrayType::size() +
-                                              v],
-                            global_vector_ptr
-                              [dof_index +
-                               (index_array_hermite[v][2 * i + 1] +
-                                comp * static_dofs_per_component) *
-                                 VectorizedArrayType::size() +
-                               v],
-                            grad_weight[v]);
+                      if (n_face_orientations == 1)
+                        {
+                          const unsigned int ind1 =
+                            index_array_hermite[0][2 * i];
+                          const unsigned int ind2 =
+                            index_array_hermite[0][2 * i + 1];
+                          AssertIndexRange(ind1,
+                                           data.dofs_per_component_on_cell);
+                          AssertIndexRange(ind2,
+                                           data.dofs_per_component_on_cell);
+                          const unsigned int i_ = reorientate(0, i);
+                          for (unsigned int comp = 0; comp < n_components;
+                               ++comp)
+                            function_1a(
+                              temp1[i_ + 2 * comp * dofs_per_face],
+                              temp1[i_ + dofs_per_face +
+                                    2 * comp * dofs_per_face],
+                              global_vector_ptr + dof_index +
+                                (ind1 + comp * static_dofs_per_component) *
+                                  VectorizedArrayType::size(),
+                              global_vector_ptr + dof_index +
+                                (ind2 + comp * static_dofs_per_component) *
+                                  VectorizedArrayType::size(),
+                              grad_weight);
+                        }
+                      else
+                        {
+                          Assert(false, ExcNotImplemented());
+                        }
                     }
                 }
-            }
-          else
-            {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
+              else
                 {
-                  if (n_face_orientations == 1)
-                    {
-                      const unsigned int i_  = reorientate(0, i);
-                      const unsigned int ind = index_array_nodal[0][i];
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        function_1b(temp1[i_ + 2 * comp * dofs_per_face],
-                                    global_vector_ptr + dof_index +
-                                      (ind + comp * static_dofs_per_component) *
-                                        VectorizedArrayType::size());
-                    }
-                  else
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
-                      Assert(false, ExcNotImplemented());
-
-                      const unsigned int n_filled_lanes =
-                        dof_info
-                          .n_vectorization_lanes_filled[dof_access_index][cell];
-
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          function_3b(
-                            temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
-                                 [v],
-                            global_vector_ptr[dof_index +
-                                              (index_array_nodal[v][i] +
-                                               comp *
-                                                 static_dofs_per_component) *
-                                                VectorizedArrayType::size() +
-                                              v]);
+                      if (n_face_orientations == 1)
+                        {
+                          const unsigned int i_  = reorientate(0, i);
+                          const unsigned int ind = index_array_nodal[0][i];
+                          for (unsigned int comp = 0; comp < n_components;
+                               ++comp)
+                            function_1b(temp1[i_ + 2 * comp * dofs_per_face],
+                                        global_vector_ptr + dof_index +
+                                          (ind +
+                                           comp * static_dofs_per_component) *
+                                            VectorizedArrayType::size());
+                        }
+                      else
+                        {
+                          Assert(false, ExcNotImplemented());
+                        }
                     }
                 }
             }
-        }
 
-      // case 2: contiguous and interleaved indices with fixed stride
-      else if (((do_gradients == false &&
-                 data.data.front().nodal_at_cell_boundaries == true) ||
-                (data.element_type ==
-                   MatrixFreeFunctions::tensor_symmetric_hermite &&
-                 fe_degree > 1)) &&
-               dof_info.index_storage_variants[dof_access_index][cell] ==
-                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                   interleaved_contiguous_strided)
-        {
-          AssertDimension(
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
-            VectorizedArrayType::size());
-          const unsigned int *indices =
-            &dof_info
-               .dof_indices_contiguous[dof_access_index]
-                                      [cell * VectorizedArrayType::size()];
-          if (fe_degree > 1 && do_gradients == true)
+          // case 2: contiguous and interleaved indices with fixed stride
+          else if (n_face_orientations == 1 &&
+                   dof_info.index_storage_variants[dof_access_index][cell] ==
+                     MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                       interleaved_contiguous_strided)
             {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
+              AssertDimension(n_face_orientations, 1);
+
+              AssertDimension(
+                dof_info.n_vectorization_lanes_filled[dof_access_index][cell],
+                VectorizedArrayType::size());
+              const unsigned int *indices =
+                &dof_info
+                   .dof_indices_contiguous[dof_access_index]
+                                          [cell * VectorizedArrayType::size()];
+              if (fe_degree > 1 && do_gradients == true)
                 {
-                  if (n_face_orientations == 1)
-                    {
-                      const unsigned int i_ = reorientate(0, i);
-
-                      const unsigned int ind1 = index_array_hermite[0][2 * i] *
-                                                VectorizedArrayType::size();
-                      const unsigned int ind2 =
-                        index_array_hermite[0][2 * i + 1] *
-                        VectorizedArrayType::size();
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        function_2a(
-                          temp1[i_ + 2 * comp * dofs_per_face],
-                          temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
-                          global_vector_ptr + ind1 +
-                            comp * static_dofs_per_component *
-                              VectorizedArrayType::size() +
-                            dof_info.component_dof_indices_offset
-                                [active_fe_index][first_selected_component] *
-                              VectorizedArrayType::size(),
-                          global_vector_ptr + ind2 +
-                            comp * static_dofs_per_component *
-                              VectorizedArrayType::size() +
-                            dof_info.component_dof_indices_offset
-                                [active_fe_index][first_selected_component] *
-                              VectorizedArrayType::size(),
-                          grad_weight,
-                          indices,
-                          indices);
-                    }
-                  else
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
-                      Assert(false, ExcNotImplemented());
-
-                      const unsigned int n_filled_lanes =
-                        dof_info
-                          .n_vectorization_lanes_filled[dof_access_index][cell];
+                      if (n_face_orientations == 1)
+                        {
+                          const unsigned int i_ = reorientate(0, i);
 
-                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                        for (unsigned int comp = 0; comp < n_components; ++comp)
-                          {
-                            const unsigned int i_ = reorientate(0, i);
-                            function_3a(temp1[i_ + 2 * comp * dofs_per_face][v],
+                          const unsigned int ind1 =
+                            index_array_hermite[0][2 * i] *
+                            VectorizedArrayType::size();
+                          const unsigned int ind2 =
+                            index_array_hermite[0][2 * i + 1] *
+                            VectorizedArrayType::size();
+                          for (unsigned int comp = 0; comp < n_components;
+                               ++comp)
+                            function_2a(temp1[i_ + 2 * comp * dofs_per_face],
                                         temp1[i_ + dofs_per_face +
-                                              2 * comp * dofs_per_face][v],
-                                        global_vector_ptr
-                                          [index_array_hermite[v][2 * i] *
-                                             VectorizedArrayType::size() +
-                                           comp * static_dofs_per_component *
-                                             VectorizedArrayType::size() +
-                                           dof_info.component_dof_indices_offset
-                                               [active_fe_index]
-                                               [first_selected_component] *
-                                             VectorizedArrayType::size() +
-                                           indices[v] + v],
-                                        global_vector_ptr
-                                          [index_array_hermite[v][2 * i + 1] *
-                                             VectorizedArrayType::size() +
-                                           comp * static_dofs_per_component *
-                                             VectorizedArrayType::size() +
-                                           dof_info.component_dof_indices_offset
-                                               [active_fe_index]
-                                               [first_selected_component] *
-                                             VectorizedArrayType::size() +
-                                           indices[v] + v],
-                                        grad_weight[v]);
-                          }
+                                              2 * comp * dofs_per_face],
+                                        global_vector_ptr + ind1 +
+                                          comp * static_dofs_per_component *
+                                            VectorizedArrayType::size() +
+                                          dof_info.component_dof_indices_offset
+                                              [active_fe_index]
+                                              [first_selected_component] *
+                                            VectorizedArrayType::size(),
+                                        global_vector_ptr + ind2 +
+                                          comp * static_dofs_per_component *
+                                            VectorizedArrayType::size() +
+                                          dof_info.component_dof_indices_offset
+                                              [active_fe_index]
+                                              [first_selected_component] *
+                                            VectorizedArrayType::size(),
+                                        grad_weight,
+                                        indices,
+                                        indices);
+                        }
+                      else
+                        {
+                          Assert(false, ExcNotImplemented());
+                        }
                     }
                 }
-            }
-          else
-            {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
+              else
                 {
-                  if (n_face_orientations == 1)
-                    {
-                      const unsigned int i_ = reorientate(0, i);
-
-                      const unsigned int ind =
-                        index_array_nodal[0][i] * VectorizedArrayType::size();
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        function_2b(
-                          temp1[i_ + 2 * comp * dofs_per_face],
-                          global_vector_ptr + ind +
-                            comp * static_dofs_per_component *
-                              VectorizedArrayType::size() +
-                            dof_info.component_dof_indices_offset
-                                [active_fe_index][first_selected_component] *
-                              VectorizedArrayType::size(),
-                          indices);
-                    }
-                  else
-                    {
-                      Assert(false, ExcNotImplemented());
-
-                      const unsigned int n_filled_lanes =
-                        dof_info
-                          .n_vectorization_lanes_filled[dof_access_index][cell];
-
-                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                        for (unsigned int comp = 0; comp < n_components; ++comp)
-                          function_3b(
-                            temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
-                                 [v],
-                            global_vector_ptr
-                              [index_array_nodal[v][i] +
-                               comp * static_dofs_per_component *
-                                 VectorizedArrayType::size() +
-                               dof_info.component_dof_indices_offset
-                                   [active_fe_index][first_selected_component] *
-                                 VectorizedArrayType::size() +
-                               indices[v] + v]);
-                    }
-                }
-            }
-        }
-
-      // case 3: contiguous and interleaved indices with mixed stride
-      else if (((do_gradients == false &&
-                 data.data.front().nodal_at_cell_boundaries == true) ||
-                (data.element_type ==
-                   MatrixFreeFunctions::tensor_symmetric_hermite &&
-                 fe_degree > 1)) &&
-               dof_info.index_storage_variants[dof_access_index][cell] ==
-                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                   interleaved_contiguous_mixed_strides)
-        {
-          const unsigned int *strides =
-            &dof_info.dof_indices_interleave_strides
-               [dof_access_index][cell * VectorizedArrayType::size()];
-          unsigned int indices[VectorizedArrayType::size()];
-          for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
-            indices[v] =
-              dof_info.dof_indices_contiguous
-                [dof_access_index][cell * VectorizedArrayType::size() + v] +
-              dof_info.component_dof_indices_offset[active_fe_index]
-                                                   [first_selected_component] *
-                strides[v];
-          const unsigned int n_filled_lanes =
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-          if (fe_degree > 1 && do_gradients == true)
-            {
-              if (n_filled_lanes == VectorizedArrayType::size())
-                for (unsigned int comp = 0; comp < n_components; ++comp)
                   for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
                       if (n_face_orientations == 1)
                         {
                           const unsigned int i_ = reorientate(0, i);
-                          unsigned int       ind1[VectorizedArrayType::size()];
-                          DEAL_II_OPENMP_SIMD_PRAGMA
-                          for (unsigned int v = 0;
-                               v < VectorizedArrayType::size();
-                               ++v)
-                            ind1[v] = indices[v] +
-                                      (comp * static_dofs_per_component +
-                                       index_array_hermite[0 /*TODO*/][2 * i]) *
-                                        strides[v];
-                          unsigned int ind2[VectorizedArrayType::size()];
-                          DEAL_II_OPENMP_SIMD_PRAGMA
-                          for (unsigned int v = 0;
-                               v < VectorizedArrayType::size();
-                               ++v)
-                            ind2[v] =
-                              indices[v] +
-                              (comp * static_dofs_per_component +
-                               index_array_hermite[0 /*TODO*/][2 * i + 1]) *
-                                strides[v];
-                          function_2a(temp1[i_ + 2 * comp * dofs_per_face],
-                                      temp1[i_ + dofs_per_face +
-                                            2 * comp * dofs_per_face],
-                                      global_vector_ptr,
-                                      global_vector_ptr,
-                                      grad_weight,
-                                      ind1,
-                                      ind2);
+
+                          const unsigned int ind = index_array_nodal[0][i] *
+                                                   VectorizedArrayType::size();
+                          for (unsigned int comp = 0; comp < n_components;
+                               ++comp)
+                            function_2b(temp1[i_ + 2 * comp * dofs_per_face],
+                                        global_vector_ptr + ind +
+                                          comp * static_dofs_per_component *
+                                            VectorizedArrayType::size() +
+                                          dof_info.component_dof_indices_offset
+                                              [active_fe_index]
+                                              [first_selected_component] *
+                                            VectorizedArrayType::size(),
+                                        indices);
                         }
                       else
                         {
                           Assert(false, ExcNotImplemented());
+                        }
+                    }
+                }
+            }
 
-                          const unsigned int n_filled_lanes =
-                            dof_info
-                              .n_vectorization_lanes_filled[dof_access_index]
-                                                           [cell];
+          // case 3: contiguous and interleaved indices with mixed stride
+          else if (n_face_orientations == 1 &&
+                   dof_info.index_storage_variants[dof_access_index][cell] ==
+                     MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                       interleaved_contiguous_mixed_strides)
+            {
+              AssertDimension(n_face_orientations, 1);
+
+              const unsigned int *strides =
+                &dof_info.dof_indices_interleave_strides
+                   [dof_access_index][cell * VectorizedArrayType::size()];
+              unsigned int indices[VectorizedArrayType::size()];
+              for (unsigned int v = 0; v < VectorizedArrayType::size(); ++v)
+                indices[v] =
+                  dof_info.dof_indices_contiguous
+                    [dof_access_index][cell * VectorizedArrayType::size() + v] +
+                  dof_info
+                      .component_dof_indices_offset[active_fe_index]
+                                                   [first_selected_component] *
+                    strides[v];
+              const unsigned int n_filled_lanes =
+                dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
 
-                          for (unsigned int v = 0; v < n_filled_lanes; ++v)
+              if (fe_degree > 1 && do_gradients == true)
+                {
+                  if (n_filled_lanes == VectorizedArrayType::size())
+                    for (unsigned int comp = 0; comp < n_components; ++comp)
+                      for (unsigned int i = 0; i < dofs_per_face; ++i)
+                        {
+                          if (n_face_orientations == 1)
+                            {
+                              const unsigned int i_ = reorientate(0, i);
+                              unsigned int ind1[VectorizedArrayType::size()];
+                              DEAL_II_OPENMP_SIMD_PRAGMA
+                              for (unsigned int v = 0;
+                                   v < VectorizedArrayType::size();
+                                   ++v)
+                                ind1[v] =
+                                  indices[v] +
+                                  (comp * static_dofs_per_component +
+                                   index_array_hermite[0 /*TODO*/][2 * i]) *
+                                    strides[v];
+                              unsigned int ind2[VectorizedArrayType::size()];
+                              DEAL_II_OPENMP_SIMD_PRAGMA
+                              for (unsigned int v = 0;
+                                   v < VectorizedArrayType::size();
+                                   ++v)
+                                ind2[v] =
+                                  indices[v] +
+                                  (comp * static_dofs_per_component +
+                                   index_array_hermite[0 /*TODO*/][2 * i + 1]) *
+                                    strides[v];
+                              function_2a(temp1[i_ + 2 * comp * dofs_per_face],
+                                          temp1[i_ + dofs_per_face +
+                                                2 * comp * dofs_per_face],
+                                          global_vector_ptr,
+                                          global_vector_ptr,
+                                          grad_weight,
+                                          ind1,
+                                          ind2);
+                            }
+                          else
                             {
-                              const unsigned int i_ = reorientate(v, i);
+                              Assert(false, ExcNotImplemented());
+                            }
+                        }
+                  else
+                    {
+                      if (integrate == false)
+                        for (unsigned int i = 0;
+                             i < n_components * 2 * dofs_per_face;
+                             ++i)
+                          temp1[i] = VectorizedArrayType();
+
+                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                        for (unsigned int comp = 0; comp < n_components; ++comp)
+                          for (unsigned int i = 0; i < dofs_per_face; ++i)
+                            {
+                              const unsigned int i_ =
+                                reorientate(n_face_orientations == 1 ? 0 : v,
+                                            i);
                               function_3a(
                                 temp1[i_ + 2 * comp * dofs_per_face][v],
                                 temp1[i_ + dofs_per_face +
@@ -2808,241 +2760,280 @@ namespace internal
                                 global_vector_ptr
                                   [indices[v] +
                                    (comp * static_dofs_per_component +
-                                    index_array_hermite[v][2 * i]) *
+                                    index_array_hermite
+                                      [n_face_orientations == 1 ? 0 : v]
+                                      [2 * i]) *
                                      strides[v]],
                                 global_vector_ptr
                                   [indices[v] +
                                    (comp * static_dofs_per_component +
-                                    index_array_hermite[v][2 * i + 1]) *
+                                    index_array_hermite
+                                      [n_face_orientations == 1 ? 0 : v]
+                                      [2 * i + 1]) *
                                      strides[v]],
-                                grad_weight[v]);
+                                grad_weight[n_face_orientations == 1 ? 0 : v]);
                             }
-                        }
                     }
+                }
               else
                 {
-                  if (integrate == false)
-                    for (unsigned int i = 0;
-                         i < n_components * 2 * dofs_per_face;
-                         ++i)
-                      temp1[i] = VectorizedArrayType();
-
-                  for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                  if (n_filled_lanes == VectorizedArrayType::size())
                     for (unsigned int comp = 0; comp < n_components; ++comp)
                       for (unsigned int i = 0; i < dofs_per_face; ++i)
                         {
-                          const unsigned int i_ =
-                            reorientate(n_face_orientations == 1 ? 0 : v, i);
-                          function_3a(
-                            temp1[i_ + 2 * comp * dofs_per_face][v],
-                            temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face]
-                                 [v],
-                            global_vector_ptr
-                              [indices[v] +
-                               (comp * static_dofs_per_component +
-                                index_array_hermite
-                                  [n_face_orientations == 1 ? 0 : v][2 * i]) *
-                                 strides[v]],
-                            global_vector_ptr
-                              [indices[v] +
-                               (comp * static_dofs_per_component +
-                                index_array_hermite[n_face_orientations == 1 ?
-                                                      0 :
-                                                      v][2 * i + 1]) *
-                                 strides[v]],
-                            grad_weight[n_face_orientations == 1 ? 0 : v]);
+                          if (n_face_orientations == 1)
+                            {
+                              unsigned int ind[VectorizedArrayType::size()];
+                              DEAL_II_OPENMP_SIMD_PRAGMA
+                              for (unsigned int v = 0;
+                                   v < VectorizedArrayType::size();
+                                   ++v)
+                                ind[v] = indices[v] +
+                                         (comp * static_dofs_per_component +
+                                          index_array_nodal[0][i]) *
+                                           strides[v];
+                              const unsigned int i_ = reorientate(0, i);
+                              function_2b(temp1[i_ + 2 * comp * dofs_per_face],
+                                          global_vector_ptr,
+                                          ind);
+                            }
+                          else
+                            {
+                              Assert(false, ExcNotImplemented());
+                            }
                         }
+                  else
+                    {
+                      if (integrate == false)
+                        for (unsigned int i = 0;
+                             i < n_components * dofs_per_face;
+                             ++i)
+                          temp1[i] = VectorizedArrayType();
+
+                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                        for (unsigned int comp = 0; comp < n_components; ++comp)
+                          for (unsigned int i = 0; i < dofs_per_face; ++i)
+                            function_3b(
+                              temp1[reorientate(
+                                      n_face_orientations == 1 ? 0 : v, i) +
+                                    2 * comp * dofs_per_face][v],
+                              global_vector_ptr
+                                [indices[v] +
+                                 (comp * static_dofs_per_component +
+                                  index_array_nodal
+                                    [n_face_orientations == 1 ? 0 : v][i]) *
+                                   strides[v]]);
+                    }
                 }
             }
-          else
+
+          // case 4: contiguous indices without interleaving
+          else if (n_face_orientations > 1 ||
+                   dof_info.index_storage_variants[dof_access_index][cell] ==
+                     MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                       contiguous)
             {
-              if (n_filled_lanes == VectorizedArrayType::size())
-                for (unsigned int comp = 0; comp < n_components; ++comp)
+              const unsigned int *indices =
+                &dof_info
+                   .dof_indices_contiguous[dof_access_index]
+                                          [cell * VectorizedArrayType::size()];
+
+              if (do_gradients == true &&
+                  data.element_type ==
+                    MatrixFreeFunctions::tensor_symmetric_hermite)
+                {
                   for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
-                      if (n_face_orientations == 1)
+                      if (n_face_orientations == 1 &&
+                          dof_info
+                              .n_vectorization_lanes_filled[dof_access_index]
+                                                           [cell] ==
+                            VectorizedArrayType::size())
                         {
-                          unsigned int ind[VectorizedArrayType::size()];
-                          DEAL_II_OPENMP_SIMD_PRAGMA
-                          for (unsigned int v = 0;
-                               v < VectorizedArrayType::size();
-                               ++v)
-                            ind[v] =
-                              indices[v] + (comp * static_dofs_per_component +
-                                            index_array_nodal[0 /*TODO*/][i]) *
-                                             strides[v];
+                          const unsigned int ind1 =
+                            index_array_hermite[0][2 * i];
+                          const unsigned int ind2 =
+                            index_array_hermite[0][2 * i + 1];
                           const unsigned int i_ = reorientate(0, i);
-                          function_2b(temp1[i_ + 2 * comp * dofs_per_face],
-                                      global_vector_ptr,
-                                      ind);
+
+                          for (unsigned int comp = 0; comp < n_components;
+                               ++comp)
+                            function_2a(
+                              temp1[i_ + 2 * comp * dofs_per_face],
+                              temp1[i_ + dofs_per_face +
+                                    2 * comp * dofs_per_face],
+                              global_vector_ptr +
+                                comp * static_dofs_per_component + ind1 +
+                                dof_info.component_dof_indices_offset
+                                  [active_fe_index][first_selected_component],
+                              global_vector_ptr +
+                                comp * static_dofs_per_component + ind2 +
+                                dof_info.component_dof_indices_offset
+                                  [active_fe_index][first_selected_component],
+                              grad_weight,
+                              indices,
+                              indices);
+                        }
+                      else if (n_face_orientations == 1)
+                        {
+                          const unsigned int ind1 =
+                            index_array_hermite[0][2 * i];
+                          const unsigned int ind2 =
+                            index_array_hermite[0][2 * i + 1];
+                          const unsigned int i_ = reorientate(0, i);
+
+                          const unsigned int n_filled_lanes =
+                            dof_info
+                              .n_vectorization_lanes_filled[dof_access_index]
+                                                           [cell];
+
+                          for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                            for (unsigned int comp = 0; comp < n_components;
+                                 ++comp)
+                              function_3a(
+                                temp1[i_ + 2 * comp * dofs_per_face][v],
+                                temp1[i_ + dofs_per_face +
+                                      2 * comp * dofs_per_face][v],
+                                global_vector_ptr
+                                  [comp * static_dofs_per_component + ind1 +
+                                   dof_info.component_dof_indices_offset
+                                     [active_fe_index]
+                                     [first_selected_component] +
+                                   indices[v]],
+                                global_vector_ptr
+                                  [comp * static_dofs_per_component + ind2 +
+                                   dof_info.component_dof_indices_offset
+                                     [active_fe_index]
+                                     [first_selected_component] +
+                                   indices[v]],
+                                grad_weight[v]);
+
+                          if (integrate == false)
+                            for (unsigned int v = n_filled_lanes;
+                                 v < VectorizedArrayType::size();
+                                 ++v)
+                              for (unsigned int comp = 0; comp < n_components;
+                                   ++comp)
+                                {
+                                  temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0;
+                                  temp1[i_ + dofs_per_face +
+                                        2 * comp * dofs_per_face][v]      = 0.0;
+                                }
                         }
                       else
                         {
                           Assert(false, ExcNotImplemented());
 
-
                           const unsigned int n_filled_lanes =
                             dof_info
                               .n_vectorization_lanes_filled[dof_access_index]
                                                            [cell];
 
                           for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                            function_3b(
-                              temp1[reorientate(v, i) +
-                                    2 * comp * dofs_per_face][v],
-                              global_vector_ptr[indices[v] +
-                                                (comp *
-                                                   static_dofs_per_component +
-                                                 index_array_nodal[v][i]) *
-                                                  strides[v]]);
+                            for (unsigned int comp = 0; comp < n_components;
+                                 ++comp)
+                              function_3a(
+                                temp1[reorientate(v, i) +
+                                      2 * comp * dofs_per_face][v],
+                                temp1[reorientate(v, i) + dofs_per_face +
+                                      2 * comp * dofs_per_face][v],
+                                global_vector_ptr
+                                  [comp * static_dofs_per_component +
+                                   index_array_hermite[v][2 * i] +
+                                   dof_info.component_dof_indices_offset
+                                     [active_fe_index]
+                                     [first_selected_component] +
+                                   indices[v]],
+                                global_vector_ptr
+                                  [comp * static_dofs_per_component +
+                                   index_array_hermite[v][2 * i + 1] +
+                                   dof_info.component_dof_indices_offset
+                                     [active_fe_index]
+                                     [first_selected_component] +
+                                   indices[v]],
+                                grad_weight[v]);
                         }
                     }
+                }
               else
                 {
-                  if (integrate == false)
-                    for (unsigned int i = 0; i < n_components * dofs_per_face;
-                         ++i)
-                      temp1[i] = VectorizedArrayType();
-
-                  for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                  for (unsigned int i = 0; i < dofs_per_face; ++i)
                     for (unsigned int comp = 0; comp < n_components; ++comp)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        function_3b(
-                          temp1[reorientate(n_face_orientations == 1 ? 0 : v,
-                                            i) +
-                                2 * comp * dofs_per_face][v],
-                          global_vector_ptr
-                            [indices[v] +
-                             (comp * static_dofs_per_component +
-                              index_array_nodal
-                                [n_face_orientations == 1 ? 0 : v][i]) *
-                               strides[v]]);
-                }
-            }
-        }
-
-      // case 4: contiguous indices without interleaving
-      else if (((do_gradients == false &&
-                 data.data.front().nodal_at_cell_boundaries == true) ||
-                (data.element_type ==
-                   MatrixFreeFunctions::tensor_symmetric_hermite &&
-                 fe_degree > 1)) &&
-               dof_info.index_storage_variants[dof_access_index][cell] ==
-                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
-                   contiguous &&
-               dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
-                 VectorizedArrayType::size())
-        {
-          const unsigned int *indices =
-            &dof_info
-               .dof_indices_contiguous[dof_access_index]
-                                      [cell * VectorizedArrayType::size()];
-
-          if (do_gradients == true &&
-              data.element_type ==
-                MatrixFreeFunctions::tensor_symmetric_hermite)
-            {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  if (n_face_orientations == 1)
-                    {
-                      const unsigned int ind1 = index_array_hermite[0][2 * i];
-                      const unsigned int ind2 =
-                        index_array_hermite[0][2 * i + 1];
-                      const unsigned int i_ = reorientate(0, i);
-
-                      for (unsigned int comp = 0; comp < n_components; ++comp)
-                        function_2a(
-                          temp1[i_ + 2 * comp * dofs_per_face],
-                          temp1[i_ + dofs_per_face + 2 * comp * dofs_per_face],
-                          global_vector_ptr + comp * static_dofs_per_component +
-                            ind1 +
-                            dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component],
-                          global_vector_ptr + comp * static_dofs_per_component +
-                            ind2 +
-                            dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component],
-                          grad_weight,
-                          indices,
-                          indices);
-                    }
-                  else
-                    {
-                      Assert(false, ExcNotImplemented());
+                      {
+                        if (n_face_orientations == 1 &&
+                            dof_info
+                                .n_vectorization_lanes_filled[dof_access_index]
+                                                             [cell] ==
+                              VectorizedArrayType::size())
+                          {
+                            const unsigned int ind = index_array_nodal[0][i];
+                            const unsigned int i_  = reorientate(0, i);
+
+                            function_2b(
+                              temp1[i_ + 2 * comp * dofs_per_face],
+                              global_vector_ptr +
+                                comp * static_dofs_per_component + ind +
+                                dof_info.component_dof_indices_offset
+                                  [active_fe_index][first_selected_component],
+                              indices);
+                          }
+                        else if (n_face_orientations == 1)
+                          {
+                            const unsigned int ind = index_array_nodal[0][i];
+                            const unsigned int i_  = reorientate(0, i);
 
-                      const unsigned int n_filled_lanes =
-                        dof_info
-                          .n_vectorization_lanes_filled[dof_access_index][cell];
+                            const unsigned int n_filled_lanes =
+                              dof_info
+                                .n_vectorization_lanes_filled[dof_access_index]
+                                                             [cell];
 
-                      for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                        for (unsigned int comp = 0; comp < n_components; ++comp)
-                          function_3a(
-                            temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
-                                 [v],
-                            temp1[reorientate(v, i) + dofs_per_face +
-                                  2 * comp * dofs_per_face][v],
-                            global_vector_ptr
-                              [comp * static_dofs_per_component +
-                               index_array_hermite[v][2 * i] +
-                               dof_info.component_dof_indices_offset
-                                 [active_fe_index][first_selected_component] +
-                               indices[v]],
-                            global_vector_ptr
-                              [comp * static_dofs_per_component +
-                               index_array_hermite[v][2 * i + 1] +
-                               dof_info.component_dof_indices_offset
-                                 [active_fe_index][first_selected_component] +
-                               indices[v]],
-                            grad_weight[v]);
-                    }
+                            for (unsigned int v = 0; v < n_filled_lanes; ++v)
+                              function_3b(
+                                temp1[i_ + 2 * comp * dofs_per_face][v],
+                                global_vector_ptr
+                                  [comp * static_dofs_per_component + ind +
+                                   dof_info.component_dof_indices_offset
+                                     [active_fe_index]
+                                     [first_selected_component] +
+                                   indices[v]]);
+
+                            if (integrate == false)
+                              for (unsigned int v = n_filled_lanes;
+                                   v < VectorizedArrayType::size();
+                                   ++v)
+                                temp1[i_ + 2 * comp * dofs_per_face][v] = 0.0;
+                          }
+                        else
+                          {
+                            for (unsigned int v = 0;
+                                 v < VectorizedArrayType::size();
+                                 ++v)
+                              if (cells[v] != numbers::invalid_unsigned_int)
+                                function_3b(
+                                  temp1[reorientate(v, i) +
+                                        2 * comp * dofs_per_face][v],
+                                  global_vector_ptr
+                                    [index_array_nodal[v][i] +
+                                     comp * static_dofs_per_component +
+                                     dof_info.component_dof_indices_offset
+                                       [active_fe_index]
+                                       [first_selected_component] +
+                                     dof_info.dof_indices_contiguous
+                                       [dof_access_index][cells[v]]]);
+                          }
+                      }
                 }
             }
           else
             {
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                for (unsigned int comp = 0; comp < n_components; ++comp)
-                  {
-                    if (n_face_orientations == 1)
-                      {
-                        const unsigned int ind = index_array_nodal[0][i];
-                        const unsigned int i_  = reorientate(0, i);
-
-                        function_2b(
-                          temp1[i_ + 2 * comp * dofs_per_face],
-                          global_vector_ptr + comp * static_dofs_per_component +
-                            ind +
-                            dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component],
-                          indices);
-                      }
-                    else
-                      {
-                        Assert(false, ExcNotImplemented());
-
-                        const unsigned int n_filled_lanes =
-                          dof_info
-                            .n_vectorization_lanes_filled[dof_access_index]
-                                                         [cell];
-
-                        for (unsigned int v = 0; v < n_filled_lanes; ++v)
-                          function_3b(
-                            temp1[reorientate(v, i) + 2 * comp * dofs_per_face]
-                                 [v],
-                            global_vector_ptr
-                              [comp * static_dofs_per_component +
-                               index_array_nodal[v][i] +
-                               dof_info.component_dof_indices_offset
-                                 [active_fe_index][first_selected_component] +
-                               indices[v]]);
-                      }
-                  }
+              // case 5: default vector access
+              function_5(temp1);
+              return false;
             }
         }
-
-      // case 5: default vector access
       else
         {
+          // case 5: default vector access
           function_5(temp1);
           return false;
         }
@@ -3051,11 +3042,11 @@ namespace internal
         function_0(temp1, dofs_per_face);
 
       if (!integrate &&
-          (face_orientation[0] > 0 &&
+          (face_orientations[0] > 0 &&
            subface_index < GeometryInfo<dim>::max_children_per_cell))
         {
-          AssertDimension(face_orientation.size(), 1);
-          adjust_for_face_orientation(face_orientation[0],
+          AssertDimension(face_orientations.size(), 1);
+          adjust_for_face_orientation(face_orientations[0],
                                       orientation_map,
                                       false,
                                       do_values,
index 488637848ff0aff8f34e8744cc7e5bebe31f04fe..21136c51f8ceebf2ac33f9670ec5450f35dc7518 100644 (file)
@@ -754,6 +754,14 @@ public:
   std::array<unsigned int, VectorizedArrayType::size()>
   get_cell_ids() const;
 
+
+  /**
+   * Return the id of the cells/faces this FEEvaluation/FEFaceEvaluation is
+   * associated with.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()>
+  get_cell_or_face_ids() const;
+
   //@}
 
   /**
@@ -3142,6 +3150,20 @@ public:
    * points is inaccurate and this value must be used instead.
    */
   const unsigned int n_q_points;
+
+
+private:
+  /**
+   * Return face number of each face of the current face batch.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()>
+  compute_face_no_data();
+
+  /**
+   * Determine the orientation of each face of the current face batch.
+   */
+  std::array<unsigned int, VectorizedArrayType::size()>
+  compute_face_orientations();
 };
 
 
@@ -3849,6 +3871,65 @@ namespace internal
 
 
 
+template <int dim,
+          int n_components_,
+          typename Number,
+          bool is_face,
+          typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEEvaluationBase<dim, n_components_, Number, is_face, VectorizedArrayType>::
+  get_cell_or_face_ids() const
+{
+  const unsigned int v_len = VectorizedArrayType::size();
+  std::array<unsigned int, VectorizedArrayType::size()> cells;
+
+  // initialize array
+  for (unsigned int i = 0; i < v_len; ++i)
+    cells[i] = numbers::invalid_unsigned_int;
+
+  if (is_face &&
+      this->dof_access_index ==
+        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+      this->is_interior_face == false)
+    {
+      // cell-based face-loop: plus face
+      for (unsigned int i = 0; i < v_len; i++)
+        {
+          // compute actual (non vectorized) cell ID
+          const unsigned int cell_this = this->cell * v_len + i;
+          // compute face ID
+          unsigned int fn =
+            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+                                                                  this->face_no,
+                                                                  i);
+
+          if (fn == numbers::invalid_unsigned_int)
+            continue; // invalid face ID: no neighbor on boundary
+
+          // get cell ID on both sides of face
+          auto cell_m = this->matrix_info->get_face_info(fn / v_len)
+                          .cells_interior[fn % v_len];
+          auto cell_p = this->matrix_info->get_face_info(fn / v_len)
+                          .cells_exterior[fn % v_len];
+
+          // compare the IDs with the given cell ID
+          if (cell_m == cell_this)
+            cells[i] = cell_p; // neighbor has the other ID
+          else if (cell_p == cell_this)
+            cells[i] = cell_m;
+        }
+    }
+  else
+    {
+      for (unsigned int i = 0; i < v_len; ++i)
+        cells[i] = cell * v_len + i;
+    }
+
+  return cells;
+}
+
+
+
 template <int dim,
           int n_components_,
           typename Number,
@@ -8142,29 +8223,68 @@ FEFaceEvaluation<dim,
     ExcMessage(
       "Only EvaluationFlags::values and EvaluationFlags::gradients are supported."));
 
-  if (!internal::FEFaceEvaluationSelector<dim,
-                                          fe_degree,
-                                          n_q_points_1d,
-                                          n_components,
-                                          Number,
-                                          VectorizedArrayType>::
-        gather_evaluate(input_vector.begin(),
-                        *this->data,
-                        *this->dof_info,
-                        this->begin_values(),
-                        this->begin_gradients(),
-                        this->scratch_data,
-                        evaluation_flag & EvaluationFlags::values,
-                        evaluation_flag & EvaluationFlags::gradients,
-                        this->active_fe_index,
-                        this->first_selected_component,
-                        this->cell,
-                        std::array<unsigned int, 1>{{this->face_no}},
-                        this->subface_index,
-                        this->dof_access_index,
-                        std::array<unsigned int, 1>{{this->face_orientation}},
-                        this->mapping_data->descriptor[this->active_fe_index]
-                          .face_orientations))
+  const auto fu = [&]() {
+    if (this->dof_access_index ==
+          internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+        this->is_interior_face == false)
+      {
+        const auto cells             = this->get_cell_or_face_ids();
+        const auto face_nos          = this->compute_face_no_data();
+        const auto face_orientations = this->compute_face_orientations();
+
+        return internal::FEFaceEvaluationSelector<dim,
+                                                  fe_degree,
+                                                  n_q_points_1d,
+                                                  n_components,
+                                                  Number,
+                                                  VectorizedArrayType>::
+          gather_evaluate(input_vector.begin(),
+                          *this->data,
+                          *this->dof_info,
+                          this->begin_values(),
+                          this->begin_gradients(),
+                          this->scratch_data,
+                          evaluation_flag & EvaluationFlags::values,
+                          evaluation_flag & EvaluationFlags::gradients,
+                          this->active_fe_index,
+                          this->first_selected_component,
+                          cells,
+                          face_nos,
+                          this->subface_index,
+                          this->dof_access_index,
+                          face_orientations,
+                          this->mapping_data->descriptor[this->active_fe_index]
+                            .face_orientations);
+      }
+    else
+      {
+        return internal::FEFaceEvaluationSelector<dim,
+                                                  fe_degree,
+                                                  n_q_points_1d,
+                                                  n_components,
+                                                  Number,
+                                                  VectorizedArrayType>::
+          gather_evaluate(input_vector.begin(),
+                          *this->data,
+                          *this->dof_info,
+                          this->begin_values(),
+                          this->begin_gradients(),
+                          this->scratch_data,
+                          evaluation_flag & EvaluationFlags::values,
+                          evaluation_flag & EvaluationFlags::gradients,
+                          this->active_fe_index,
+                          this->first_selected_component,
+                          std::array<unsigned int, 1>{{this->cell}},
+                          std::array<unsigned int, 1>{{this->face_no}},
+                          this->subface_index,
+                          this->dof_access_index,
+                          std::array<unsigned int, 1>{{this->face_orientation}},
+                          this->mapping_data->descriptor[this->active_fe_index]
+                            .face_orientations);
+      }
+  };
+
+  if (!fu())
     {
       this->read_dof_values(input_vector);
       this->evaluate(evaluation_flag);
@@ -8235,6 +8355,11 @@ FEFaceEvaluation<dim,
                 "Use integrate() followed by distribute_local_to_global() "
                 "instead.");
 
+  Assert((this->dof_access_index ==
+            internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
+          this->is_interior_face == false) == false,
+         ExcNotImplemented());
+
   if (!internal::FEFaceEvaluationSelector<dim,
                                           fe_degree,
                                           n_q_points_1d,
@@ -8252,7 +8377,7 @@ FEFaceEvaluation<dim,
                           evaluation_flag & EvaluationFlags::gradients,
                           this->active_fe_index,
                           this->first_selected_component,
-                          this->cell,
+                          std::array<unsigned int, 1>{{this->cell}},
                           std::array<unsigned int, 1>{{this->face_no}},
                           this->subface_index,
                           this->dof_access_index,
@@ -8318,6 +8443,161 @@ FEFaceEvaluation<dim,
 
 
 
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::compute_face_no_data()
+{
+  std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
+
+  if (this->dof_access_index !=
+        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
+      this->is_interior_face == true)
+    {
+      std::fill(face_no_data.begin(),
+                face_no_data.begin() +
+                  this->dof_info->n_vectorization_lanes_filled
+                    [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+                    [this->cell],
+                this->face_no);
+    }
+  else
+    {
+      std::fill(face_no_data.begin(),
+                face_no_data.end(),
+                numbers::invalid_unsigned_int);
+
+      for (unsigned int i = 0;
+           i < this->dof_info->n_vectorization_lanes_filled
+                 [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+                 [this->cell];
+           i++)
+        {
+          // compute actual (non vectorized) cell ID
+          const unsigned int cell_this =
+            this->cell * VectorizedArrayType::size() + i;
+          // compute face ID
+          const unsigned int face_index =
+            this->matrix_info->get_cell_and_face_to_plain_faces()(this->cell,
+                                                                  this->face_no,
+                                                                  i);
+
+          Assert(face_index != numbers::invalid_unsigned_int,
+                 ExcNotInitialized());
+
+          // get cell ID on both sides of face
+          auto cell_m =
+            this->matrix_info
+              ->get_face_info(face_index / VectorizedArrayType::size())
+              .cells_interior[face_index % VectorizedArrayType::size()];
+
+          // compare the IDs with the given cell ID
+          face_no_data[i] =
+            (cell_m == cell_this) ?
+              this->matrix_info
+                ->get_face_info(face_index / VectorizedArrayType::size())
+                .exterior_face_no :
+              this->matrix_info
+                ->get_face_info(face_index / VectorizedArrayType::size())
+                .interior_face_no;
+        }
+    }
+
+  return face_no_data;
+}
+
+
+
+template <int dim,
+          int fe_degree,
+          int n_q_points_1d,
+          int n_components_,
+          typename Number,
+          typename VectorizedArrayType>
+std::array<unsigned int, VectorizedArrayType::size()>
+FEFaceEvaluation<dim,
+                 fe_degree,
+                 n_q_points_1d,
+                 n_components_,
+                 Number,
+                 VectorizedArrayType>::compute_face_orientations()
+{
+  std::array<unsigned int, VectorizedArrayType::size()> face_no_data;
+
+  if (this->dof_access_index !=
+        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell ||
+      this->is_interior_face == true)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+  else
+    {
+      std::fill(face_no_data.begin(), face_no_data.end(), 0);
+
+      if (dim == 3)
+        {
+          for (unsigned int i = 0;
+               i < this->dof_info->n_vectorization_lanes_filled
+                     [internal::MatrixFreeFunctions::DoFInfo::dof_access_cell]
+                     [this->cell];
+               i++)
+            {
+              // compute actual (non vectorized) cell ID
+              const unsigned int cell_this =
+                this->cell * VectorizedArrayType::size() + i;
+              // compute face ID
+              const unsigned int face_index =
+                this->matrix_info->get_cell_and_face_to_plain_faces()(
+                  this->cell, this->face_no, i);
+
+              Assert(face_index != numbers::invalid_unsigned_int,
+                     ExcNotInitialized());
+
+              const unsigned int macro =
+                face_index / VectorizedArrayType::size();
+              const unsigned int lane =
+                face_index % VectorizedArrayType::size();
+
+              const auto &faces = this->matrix_info->get_face_info(macro);
+
+              // get cell ID on both sides of face
+              auto cell_m = faces.cells_interior[lane];
+
+              const bool is_interior_face = cell_m != cell_this;
+              const bool fo_interior_face = faces.face_orientation >= 8;
+
+              unsigned int face_orientation = faces.face_orientation % 8;
+
+              if (is_interior_face != fo_interior_face)
+                {
+                  // invert (see also:
+                  // Triangulation::update_periodic_face_map())
+                  static const std::array<unsigned int, 8> table{
+                    {0, 1, 0, 3, 6, 5, 4, 7}};
+
+                  face_orientation = table[face_orientation];
+                }
+
+              // compare the IDs with the given cell ID
+              face_no_data[i] = face_orientation;
+            }
+        }
+    }
+
+  return face_no_data;
+}
+
+
+
 /*------------------------- end FEFaceEvaluation ------------------------- */
 
 
diff --git a/tests/matrix_free/ecl_02.cc b/tests/matrix_free/ecl_02.cc
new file mode 100644 (file)
index 0000000..8a0df9a
--- /dev/null
@@ -0,0 +1,153 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/matrix_free/fe_evaluation.h>
+#include <deal.II/matrix_free/matrix_free.h>
+
+#include <deal.II/multigrid/mg_constrained_dofs.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include "../tests.h"
+
+
+// tests matrix-free read_cell_data for locally refined mesh
+
+template <int dim,
+          int fe_degree,
+          int n_points                 = fe_degree + 1,
+          typename Number              = double,
+          typename VectorizedArrayType = VectorizedArray<Number>>
+void
+test(const unsigned int n_refinements = 1)
+{
+  using VectorType = LinearAlgebra::distributed::Vector<Number>;
+
+  Triangulation<dim> tria;
+  GridGenerator::subdivided_hyper_rectangle(
+    tria, {2, 1}, Point<dim>(0.0, 0.0), Point<dim>(2.0, 1.0), true);
+
+  std::vector<dealii::GridTools::PeriodicFacePair<
+    typename dealii::Triangulation<dim>::cell_iterator>>
+    periodic_faces;
+
+  if (dim >= 1)
+    dealii::GridTools::collect_periodic_faces(tria, 0, 1, 0, periodic_faces);
+
+  if (dim >= 2)
+    dealii::GridTools::collect_periodic_faces(tria, 2, 3, 1, periodic_faces);
+
+  tria.add_periodicity(periodic_faces);
+
+  FE_DGQ<dim>     fe(fe_degree);
+  DoFHandler<dim> dof_handler(tria);
+  dof_handler.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  QGauss<1>     quad(n_points);
+
+  AffineConstraints<Number> constraint;
+
+  using MF = MatrixFree<dim, Number, VectorizedArrayType>;
+
+  typename MF::AdditionalData additional_data;
+  additional_data.mapping_update_flags                = update_values;
+  additional_data.mapping_update_flags_inner_faces    = update_values;
+  additional_data.mapping_update_flags_boundary_faces = update_values;
+  additional_data.mapping_update_flags_faces_by_cells = update_values;
+  additional_data.hold_all_faces_to_owned_cells       = true;
+
+  MF matrix_free;
+  matrix_free.reinit(mapping, dof_handler, constraint, quad, additional_data);
+
+  VectorType src, dst;
+
+  matrix_free.initialize_dof_vector(src);
+  matrix_free.initialize_dof_vector(dst);
+
+  FEEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType> phi(
+    matrix_free);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_m(matrix_free, true);
+  FEFaceEvaluation<dim, fe_degree, n_points, 1, Number, VectorizedArrayType>
+    phi_p(matrix_free, false);
+
+
+  for (unsigned int i = 0; i < src.size() / 2; i++)
+    src[i] = 1;
+
+  for (unsigned int i = src.size() / 2; i < src.size(); i++)
+    src[i] = 2;
+
+  dst = 0.0;
+
+  /**
+   * Element-centric loop
+   */
+  matrix_free.template loop_cell_centric<VectorType, VectorType>(
+    [&](const auto &, auto &dst, const auto &src, const auto range) {
+      for (unsigned int cell = range.first; cell < range.second; ++cell)
+        {
+          for (unsigned int face = 0; face < GeometryInfo<dim>::faces_per_cell;
+               face++)
+            {
+              phi_m.reinit(cell, face);
+              phi_p.reinit(cell, face);
+
+              phi_m.read_dof_values(src);
+
+              for (unsigned int i = 0; i < phi_m.static_dofs_per_component; i++)
+                deallog << static_cast<int>(phi_m.begin_dof_values()[i][0])
+                        << " ";
+              deallog << std::endl;
+              phi_m.gather_evaluate(src, true, false);
+              for (unsigned int i = 0; i < phi_m.static_n_q_points; i++)
+                deallog << static_cast<int>(phi_m.begin_values()[i][0]) << " ";
+              deallog << std::endl;
+
+              phi_p.read_dof_values(src);
+              for (unsigned int i = 0; i < phi_p.static_dofs_per_component; i++)
+                deallog << static_cast<int>(phi_p.begin_dof_values()[i][0])
+                        << " ";
+              deallog << std::endl;
+
+              phi_p.gather_evaluate(src, true, false);
+              for (unsigned int i = 0; i < phi_p.static_n_q_points; i++)
+                deallog << static_cast<int>(phi_p.begin_values()[i][0]) << " ";
+              deallog << std::endl << std::endl;
+            }
+        }
+    },
+    dst,
+    src,
+    false,
+    MF::DataAccessOnFaces::gradients);
+}
+
+int
+main()
+{
+  initlog();
+  test<2, 1, 2, double, VectorizedArray<double, 2>>();
+}
diff --git a/tests/matrix_free/ecl_02.output b/tests/matrix_free/ecl_02.output
new file mode 100644 (file)
index 0000000..d584d53
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::1 1 1 1 
+DEAL::2 2 
+DEAL::
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::1 1 1 1 
+DEAL::2 2 
+DEAL::
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::1 1 1 1 
+DEAL::1 1 
+DEAL::

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.