]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Merge implementation of gather_evaluate and integrate_scatter 10811/head
authorPeter Munch <peterrmuench@gmail.com>
Sun, 9 Aug 2020 17:49:39 +0000 (19:49 +0200)
committerPeter Munch <peterrmuench@gmail.com>
Mon, 10 Aug 2020 06:45:24 +0000 (08:45 +0200)
include/deal.II/matrix_free/evaluation_kernels.h

index 88cc9bfef2f4f3cdebe0624da996ea07b672abb1..1a101a63d01739c79a71a2c5b8bc5c8dd2cab6d8 100644 (file)
@@ -2075,457 +2075,108 @@ namespace internal
       const unsigned int                                 face_orientation,
       const Table<2, unsigned int> &                     orientation_map)
     {
-      const unsigned int side = face_no % 2;
-
-      constexpr unsigned int static_dofs_per_component =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
-                         numbers::invalid_unsigned_int;
-      constexpr unsigned int static_dofs_per_face =
-        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
-                         numbers::invalid_unsigned_int;
-      const unsigned int dofs_per_face =
-        fe_degree > -1 ?
-          static_dofs_per_face :
-          Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
-
-      constexpr unsigned int stack_array_size_threshold = 100;
-
-      VectorizedArrayType
-        temp_data[static_dofs_per_face < stack_array_size_threshold ?
-                    n_components * 2 * dofs_per_face :
-                    1];
-      VectorizedArrayType *__restrict temp1;
-      if (static_dofs_per_face < stack_array_size_threshold)
-        temp1 = &temp_data[0];
-      else
-        temp1 = scratch_data;
-
-      // case 1: contiguous and interleaved indices
-      if (((evaluate_gradients == false &&
-            data.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)
-        {
-          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();
-
-          if (fe_degree > 1 && evaluate_gradients == true)
-            {
-              // we know that the gradient weights for the Hermite case on the
-              // right (side==1) are the negative from the value at the left
-              // (side==0), so we only read out one of them.
-              const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
-              AssertDimension(data.face_to_cell_index_hermite.size(1),
-                              2 * dofs_per_face);
-              const unsigned int *index_array =
-                &data.face_to_cell_index_hermite(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind1 = index_array[2 * i];
-                  const unsigned int ind2 = index_array[2 * i + 1];
-                  AssertIndexRange(ind1, data.dofs_per_component_on_cell);
-                  AssertIndexRange(ind2, data.dofs_per_component_on_cell);
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      do_vectorized_read(src_ptr + dof_index +
-                                           (ind1 +
-                                            comp * static_dofs_per_component) *
-                                             VectorizedArrayType::size(),
-                                         temp1[i + 2 * comp * dofs_per_face]);
-                      do_vectorized_read(
-                        src_ptr + dof_index +
-                          (ind2 + comp * static_dofs_per_component) *
-                            VectorizedArrayType::size(),
-                        temp1[dofs_per_face + i + 2 * comp * dofs_per_face]);
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                        grad_weight *
-                        (temp1[i + 2 * comp * dofs_per_face] -
-                         temp1[i + dofs_per_face + 2 * comp * dofs_per_face]);
-                    }
-                }
-            }
-          else
-            {
-              AssertDimension(data.face_to_cell_index_nodal.size(1),
-                              dofs_per_face);
-              const unsigned int *index_array =
-                &data.face_to_cell_index_nodal(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind = index_array[i];
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    do_vectorized_read(src_ptr + dof_index +
-                                         (ind +
-                                          comp * static_dofs_per_component) *
-                                           VectorizedArrayType::size(),
-                                       temp1[i + 2 * comp * dofs_per_face]);
-                }
-            }
-        }
-
-      // case 2: contiguous and interleaved indices with fixed stride
-      else if (((evaluate_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 && evaluate_gradients == true)
-            {
-              // we know that the gradient weights for the Hermite case on the
-              // right (side==1) are the negative from the value at the left
-              // (side==0), so we only read out one of them.
-              const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
-              AssertDimension(data.face_to_cell_index_hermite.size(1),
-                              2 * dofs_per_face);
-
-              const unsigned int *index_array =
-                &data.face_to_cell_index_hermite(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind1 =
-                    index_array[2 * i] * VectorizedArrayType::size();
-                  const unsigned int ind2 =
-                    index_array[2 * i + 1] * VectorizedArrayType::size();
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      do_vectorized_gather(
-                        src_ptr + ind1 +
-                          comp * static_dofs_per_component *
-                            VectorizedArrayType::size() +
-                          dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component] *
-                            VectorizedArrayType::size(),
-                        indices,
-                        temp1[i + 2 * comp * dofs_per_face]);
-                      VectorizedArrayType grad;
-                      do_vectorized_gather(
-                        src_ptr + ind2 +
-                          comp * static_dofs_per_component *
-                            VectorizedArrayType::size() +
-                          dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component] *
-                            VectorizedArrayType::size(),
-                        indices,
-                        grad);
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                        grad_weight *
-                        (temp1[i + 2 * comp * dofs_per_face] - grad);
-                    }
-                }
-            }
-          else
-            {
-              AssertDimension(data.face_to_cell_index_nodal.size(1),
-                              dofs_per_face);
-              const unsigned int *index_array =
-                &data.face_to_cell_index_nodal(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind =
-                    index_array[i] * VectorizedArrayType::size();
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    do_vectorized_gather(
-                      src_ptr + ind +
-                        comp * static_dofs_per_component *
-                          VectorizedArrayType::size() +
-                        dof_info.component_dof_indices_offset
-                            [active_fe_index][first_selected_component] *
-                          VectorizedArrayType::size(),
-                      indices,
-                      temp1[i + 2 * comp * dofs_per_face]);
-                }
-            }
-        }
-
-      // case 3: contiguous and interleaved indices with mixed stride
-      else if (((evaluate_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 nvec =
-            dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
-
-          if (fe_degree > 1 && evaluate_gradients == true)
-            {
-              // we know that the gradient weights for the Hermite case on the
-              // right (side==1) are the negative from the value at the left
-              // (side==0), so we only read out one of them.
-              const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
-              AssertDimension(data.face_to_cell_index_hermite.size(1),
-                              2 * dofs_per_face);
-
-              const unsigned int *index_array =
-                &data.face_to_cell_index_hermite(face_no, 0);
-              if (nvec == VectorizedArrayType::size())
-                for (unsigned int comp = 0; comp < n_components; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++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[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[2 * i + 1]) *
-                                         strides[v];
-                      do_vectorized_gather(src_ptr,
-                                           ind1,
-                                           temp1[i + 2 * comp * dofs_per_face]);
-                      VectorizedArrayType grad;
-                      do_vectorized_gather(src_ptr, ind2, grad);
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                        grad_weight *
-                        (temp1[i + 2 * comp * dofs_per_face] - grad);
-                    }
-              else
-                {
-                  for (unsigned int i = 0; i < n_components * 2 * dofs_per_face;
-                       ++i)
-                    temp1[i] = VectorizedArrayType();
-                  for (unsigned int v = 0; v < nvec; ++v)
-                    for (unsigned int comp = 0; comp < n_components; ++comp)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        {
-                          const unsigned int ind1 =
-                            indices[v] + (comp * static_dofs_per_component +
-                                          index_array[2 * i]) *
-                                           strides[v];
-                          const unsigned int ind2 =
-                            indices[v] + (comp * static_dofs_per_component +
-                                          index_array[2 * i + 1]) *
-                                           strides[v];
-                          temp1[i + 2 * comp * dofs_per_face][v] =
-                            src_ptr[ind1];
-                          const Number grad = src_ptr[ind2];
-                          temp1[i + dofs_per_face +
-                                2 * comp * dofs_per_face][v] =
-                            grad_weight[0] *
-                            (temp1[i + 2 * comp * dofs_per_face][v] - grad);
-                        }
-                }
-            }
-          else
-            {
-              AssertDimension(data.face_to_cell_index_nodal.size(1),
-                              dofs_per_face);
-              const unsigned int *index_array =
-                &data.face_to_cell_index_nodal(face_no, 0);
-              if (nvec == VectorizedArrayType::size())
-                for (unsigned int comp = 0; comp < n_components; ++comp)
-                  for (unsigned int i = 0; i < dofs_per_face; ++i)
-                    {
-                      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[i]) *
-                            strides[v];
-                      do_vectorized_gather(src_ptr,
-                                           ind,
-                                           temp1[i + 2 * comp * dofs_per_face]);
-                    }
-              else
-                {
-                  for (unsigned int i = 0; i < n_components * dofs_per_face;
-                       ++i)
-                    temp1[i] = VectorizedArrayType();
-                  for (unsigned int v = 0; v < nvec; ++v)
-                    for (unsigned int comp = 0; comp < n_components; ++comp)
-                      for (unsigned int i = 0; i < dofs_per_face; ++i)
-                        {
-                          const unsigned int ind1 =
-                            indices[v] + (comp * static_dofs_per_component +
-                                          index_array[i]) *
-                                           strides[v];
-                          temp1[i + 2 * comp * dofs_per_face][v] =
-                            src_ptr[ind1];
-                        }
-                }
-            }
-        }
-
-      // case 4: contiguous indices without interleaving
-      else if (((evaluate_gradients == false &&
-                 data.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 (evaluate_gradients == true &&
-              data.element_type ==
-                MatrixFreeFunctions::tensor_symmetric_hermite)
-            {
-              // we know that the gradient weights for the Hermite case on the
-              // right (side==1) are the negative from the value at the left
-              // (side==0), so we only read out one of them.
-              const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 1 + side];
-              AssertDimension(data.face_to_cell_index_hermite.size(1),
-                              2 * dofs_per_face);
-
-              const unsigned int *index_array =
-                &data.face_to_cell_index_hermite(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind1 = index_array[2 * i];
-                  const unsigned int ind2 = index_array[2 * i + 1];
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      do_vectorized_gather(
-                        src_ptr + ind1 + comp * static_dofs_per_component +
-                          dof_info.component_dof_indices_offset
-                            [active_fe_index][first_selected_component],
-                        indices,
-                        temp1[i + 2 * comp * dofs_per_face]);
-                      VectorizedArrayType grad;
-                      do_vectorized_gather(
-                        src_ptr + ind2 + comp * static_dofs_per_component +
-                          dof_info.component_dof_indices_offset
-                            [active_fe_index][first_selected_component],
-                        indices,
-                        grad);
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face] =
-                        grad_weight *
-                        (temp1[i + 2 * comp * dofs_per_face] - grad);
-                    }
-                }
-            }
+      return process_and_io( //
+        false /*=evaluate*/,
+        src_ptr,
+        data,
+        dof_info,
+        values_quad,
+        gradients_quad,
+        scratch_data,
+        evaluate_values,
+        evaluate_gradients,
+        active_fe_index,
+        first_selected_component,
+        cell,
+        face_no,
+        subface_index,
+        dof_access_index,
+        face_orientation,
+        orientation_map,
+        [](auto &      temp_1,
+           auto &      temp_2,
+           const auto  src_ptr_1,
+           const auto  src_ptr_2,
+           const auto &grad_weight) {
+          // case 1a)
+          do_vectorized_read(src_ptr_1, temp_1);
+          do_vectorized_read(src_ptr_2, temp_2);
+          temp_2 = grad_weight * (temp_1 - temp_2);
+        },
+        [](auto &temp, const auto src_ptr) {
+          // case 1b)
+          do_vectorized_read(src_ptr, temp);
+        },
+        [](auto &      temp_1,
+           auto &      temp_2,
+           const auto  src_ptr_1,
+           const auto  src_ptr_2,
+           const auto &grad_weight,
+           const auto &indices_1,
+           const auto &indices_2) {
+          // case 2a)
+          do_vectorized_gather(src_ptr_1, indices_1, temp_1);
+          do_vectorized_gather(src_ptr_2, indices_2, temp_2);
+          temp_2 = grad_weight * (temp_1 - temp_2);
+        },
+        [](auto &temp, const auto src_ptr, const auto &indices) {
+          // case 2b)
+          do_vectorized_gather(src_ptr, indices, temp);
+        },
+        [](auto &      temp_1,
+           auto &      temp_2,
+           const auto &src_ptr_1,
+           const auto &src_ptr_2,
+           const auto &grad_weight) {
+          // case 3a)
+          temp_1 = src_ptr_1;
+          temp_2 = grad_weight * (temp_1 - src_ptr_2);
+        },
+        [](auto &temp, const auto &src_ptr) {
+          // case 3b)
+          temp = src_ptr;
+        },
+        [&](const auto &) {
+          // case 5)
+        },
+        [&](auto &temp1, const auto &dofs_per_face) {
+          if (fe_degree > -1 &&
+              subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+              data.element_type <= MatrixFreeFunctions::tensor_symmetric)
+            FEFaceEvaluationImpl<
+              true,
+              dim,
+              fe_degree,
+              n_q_points_1d,
+              n_components,
+              VectorizedArrayType>::evaluate_in_face(data,
+                                                     temp1,
+                                                     values_quad,
+                                                     gradients_quad,
+                                                     scratch_data +
+                                                       2 * n_components *
+                                                         dofs_per_face,
+                                                     evaluate_values,
+                                                     evaluate_gradients,
+                                                     subface_index);
           else
-            {
-              AssertDimension(data.face_to_cell_index_nodal.size(1),
-                              dofs_per_face);
-              const unsigned int *index_array =
-                &data.face_to_cell_index_nodal(face_no, 0);
-              for (unsigned int i = 0; i < dofs_per_face; ++i)
-                for (unsigned int comp = 0; comp < n_components; ++comp)
-                  {
-                    const unsigned int ind = index_array[i];
-                    do_vectorized_gather(
-                      src_ptr + ind + comp * static_dofs_per_component +
-                        dof_info.component_dof_indices_offset
-                          [active_fe_index][first_selected_component],
-                      indices,
-                      temp1[i + comp * 2 * dofs_per_face]);
-                  }
-            }
-        }
-
-      // case 5: default vector access
-      else
-        {
-          return false;
-        }
-
-      if (fe_degree > -1 &&
-          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-          data.element_type <= MatrixFreeFunctions::tensor_symmetric)
-        FEFaceEvaluationImpl<
-          true,
-          dim,
-          fe_degree,
-          n_q_points_1d,
-          n_components,
-          VectorizedArrayType>::evaluate_in_face(data,
-                                                 temp1,
-                                                 values_quad,
-                                                 gradients_quad,
-                                                 scratch_data + 2 *
-                                                                  n_components *
-                                                                  dofs_per_face,
-                                                 evaluate_values,
-                                                 evaluate_gradients,
-                                                 subface_index);
-      else
-        FEFaceEvaluationImpl<
-          false,
-          dim,
-          fe_degree,
-          n_q_points_1d,
-          n_components,
-          VectorizedArrayType>::evaluate_in_face(data,
-                                                 temp1,
-                                                 values_quad,
-                                                 gradients_quad,
-                                                 scratch_data + 2 *
-                                                                  n_components *
-                                                                  dofs_per_face,
-                                                 evaluate_values,
-                                                 evaluate_gradients,
-                                                 subface_index);
-
-      if (face_orientation)
-        adjust_for_face_orientation(face_orientation,
-                                    orientation_map,
-                                    false,
-                                    evaluate_values,
-                                    evaluate_gradients,
-                                    data.n_q_points_face,
-                                    scratch_data,
-                                    values_quad,
-                                    gradients_quad);
-
-      return true;
+            FEFaceEvaluationImpl<
+              false,
+              dim,
+              fe_degree,
+              n_q_points_1d,
+              n_components,
+              VectorizedArrayType>::evaluate_in_face(data,
+                                                     temp1,
+                                                     values_quad,
+                                                     gradients_quad,
+                                                     scratch_data +
+                                                       2 * n_components *
+                                                         dofs_per_face,
+                                                     evaluate_values,
+                                                     evaluate_gradients,
+                                                     subface_index);
+        });
     }
 
     static bool
@@ -2548,82 +2199,211 @@ namespace internal
       const unsigned int                                 face_orientation,
       const Table<2, unsigned int> &                     orientation_map)
     {
-      if (face_orientation)
+      return process_and_io( //
+        true /*=integrate*/,
+        dst_ptr,
+        data,
+        dof_info,
+        values_quad,
+        gradients_quad,
+        scratch_data,
+        integrate_values,
+        integrate_gradients,
+        active_fe_index,
+        first_selected_component,
+        cell,
+        face_no,
+        subface_index,
+        dof_access_index,
+        face_orientation,
+        orientation_map,
+        [](const auto &temp_1,
+           const auto &temp_2,
+           auto        dst_ptr_1,
+           auto        dst_ptr_2,
+           const auto &grad_weight) {
+          // case 1a)
+          const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
+          const VectorizedArrayType grad = grad_weight * temp_2;
+          do_vectorized_add(val, dst_ptr_1);
+          do_vectorized_add(grad, dst_ptr_2);
+        },
+        [](const auto &temp, auto dst_ptr) {
+          // case 1b)
+          do_vectorized_add(temp, dst_ptr);
+        },
+        [](const auto &temp_1,
+           const auto &temp_2,
+           auto        dst_ptr_1,
+           auto        dst_ptr_2,
+           const auto &grad_weight,
+           const auto &indices_1,
+           const auto &indices_2) {
+          // case 2a)
+          const VectorizedArrayType val  = temp_1 - grad_weight * temp_2;
+          const VectorizedArrayType grad = grad_weight * temp_2;
+          do_vectorized_scatter_add(val, indices_1, dst_ptr_1);
+          do_vectorized_scatter_add(grad, indices_2, dst_ptr_2);
+        },
+        [](const auto &temp, auto dst_ptr, const auto &indices) {
+          // case 2b)
+          do_vectorized_scatter_add(temp, indices, dst_ptr);
+        },
+        [](const auto &temp_1,
+           const auto &temp_2,
+           auto &      dst_ptr_1,
+           auto &      dst_ptr_2,
+           const auto &grad_weight) {
+          // case 3a)
+          const Number val  = temp_1 - grad_weight * temp_2;
+          const Number grad = grad_weight * temp_2;
+          dst_ptr_1 += val;
+          dst_ptr_2 += grad;
+        },
+        [](const auto &temp, auto &dst_ptr) {
+          // case 3b)
+          dst_ptr += temp;
+        },
+        [&](const auto &temp1) {
+          // case 5: default vector access, must be handled separately, just do
+          // the face-normal interpolation
+          FEFaceNormalEvaluationImpl<dim,
+                                     fe_degree,
+                                     n_components,
+                                     VectorizedArrayType>::
+            template interpolate<false, false>(
+              data, temp1, values_array, integrate_gradients, face_no);
+        },
+        [&](auto &temp1, const auto &dofs_per_face) {
+          if (fe_degree > -1 &&
+              subface_index >= GeometryInfo<dim>::max_children_per_cell &&
+              data.element_type <=
+                internal::MatrixFreeFunctions::tensor_symmetric)
+            internal::FEFaceEvaluationImpl<
+              true,
+              dim,
+              fe_degree,
+              n_q_points_1d,
+              n_components,
+              VectorizedArrayType>::integrate_in_face(data,
+                                                      temp1,
+                                                      values_quad,
+                                                      gradients_quad,
+                                                      scratch_data +
+                                                        2 * n_components *
+                                                          dofs_per_face,
+                                                      integrate_values,
+                                                      integrate_gradients,
+                                                      subface_index);
+          else
+            internal::FEFaceEvaluationImpl<
+              false,
+              dim,
+              fe_degree,
+              n_q_points_1d,
+              n_components,
+              VectorizedArrayType>::integrate_in_face(data,
+                                                      temp1,
+                                                      values_quad,
+                                                      gradients_quad,
+                                                      scratch_data +
+                                                        2 * n_components *
+                                                          dofs_per_face,
+                                                      integrate_values,
+                                                      integrate_gradients,
+                                                      subface_index);
+        });
+    }
+
+  private:
+    template <typename Number2_,
+              typename Function1a,
+              typename Function1b,
+              typename Function2a,
+              typename Function2b,
+              typename Function3a,
+              typename Function3b,
+              typename Function5,
+              typename Function0>
+    static bool
+    process_and_io(
+      const bool integrate,
+      Number2_ * global_vector_ptr,
+      const MatrixFreeFunctions::ShapeInfo<VectorizedArrayType> &data,
+      const MatrixFreeFunctions::DoFInfo &                       dof_info,
+      VectorizedArrayType *                                      values_quad,
+      VectorizedArrayType *                                      gradients_quad,
+      VectorizedArrayType *                                      scratch_data,
+      const bool                                                 do_values,
+      const bool                                                 do_gradients,
+      const unsigned int active_fe_index,
+      const unsigned int first_selected_component,
+      const unsigned int cell,
+      const unsigned int face_no,
+      const unsigned int subface_index,
+      const MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index,
+      const unsigned int                                 face_orientation,
+      const Table<2, unsigned int> &                     orientation_map,
+      const Function1a &                                 function_1a,
+      const Function1b &                                 function_1b,
+      const Function2a &                                 function_2a,
+      const Function2b &                                 function_2b,
+      const Function3a &                                 function_3a,
+      const Function3b &                                 function_3b,
+      const Function5 &                                  function_5,
+      const Function0 &                                  function_0)
+    {
+      (void)subface_index;
+
+      if (integrate && (face_orientation > 0))
         adjust_for_face_orientation(face_orientation,
                                     orientation_map,
                                     true,
-                                    integrate_values,
-                                    integrate_gradients,
+                                    do_values,
+                                    do_gradients,
                                     data.n_q_points_face,
                                     scratch_data,
                                     values_quad,
                                     gradients_quad);
 
-      const unsigned int     side = face_no % 2;
+      const unsigned int side = face_no % 2;
+
+      const unsigned int side_ = integrate ? (2 - side) : (1 + side);
+
       constexpr unsigned int static_dofs_per_component =
         fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim) :
                          numbers::invalid_unsigned_int;
+      constexpr unsigned int static_dofs_per_face =
+        fe_degree > -1 ? Utilities::pow(fe_degree + 1, dim - 1) :
+                         numbers::invalid_unsigned_int;
       const unsigned int dofs_per_face =
         fe_degree > -1 ?
-          Utilities::pow(fe_degree + 1, dim - 1) :
+          static_dofs_per_face :
           Utilities::pow(data.data.front().fe_degree + 1, dim - 1);
 
       constexpr unsigned int stack_array_size_threshold = 100;
 
-      VectorizedArrayType temp_data[dofs_per_face < stack_array_size_threshold ?
-                                      n_components * 2 * dofs_per_face :
-                                      1];
+      VectorizedArrayType
+        temp_data[static_dofs_per_face < stack_array_size_threshold ?
+                    n_components * 2 * dofs_per_face :
+                    1];
       VectorizedArrayType *__restrict temp1;
-      if (dofs_per_face < stack_array_size_threshold)
+      if (static_dofs_per_face < stack_array_size_threshold)
         temp1 = &temp_data[0];
       else
         temp1 = scratch_data;
 
-      if (fe_degree > -1 &&
-          subface_index >= GeometryInfo<dim>::max_children_per_cell &&
-          data.element_type <= internal::MatrixFreeFunctions::tensor_symmetric)
-        internal::FEFaceEvaluationImpl<
-          true,
-          dim,
-          fe_degree,
-          n_q_points_1d,
-          n_components,
-          VectorizedArrayType>::integrate_in_face(data,
-                                                  temp1,
-                                                  values_quad,
-                                                  gradients_quad,
-                                                  scratch_data +
-                                                    2 * n_components *
-                                                      dofs_per_face,
-                                                  integrate_values,
-                                                  integrate_gradients,
-                                                  subface_index);
-      else
-        internal::FEFaceEvaluationImpl<
-          false,
-          dim,
-          fe_degree,
-          n_q_points_1d,
-          n_components,
-          VectorizedArrayType>::integrate_in_face(data,
-                                                  temp1,
-                                                  values_quad,
-                                                  gradients_quad,
-                                                  scratch_data +
-                                                    2 * n_components *
-                                                      dofs_per_face,
-                                                  integrate_values,
-                                                  integrate_gradients,
-                                                  subface_index);
+      if (integrate)
+        function_0(temp1, dofs_per_face);
 
       // case 1: contiguous and interleaved indices
-      if (((integrate_gradients == false &&
+      if (((do_gradients == false &&
             data.data.front().nodal_at_cell_boundaries == true) ||
            (data.element_type ==
-              internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+              MatrixFreeFunctions::tensor_symmetric_hermite &&
             fe_degree > 1)) &&
           dof_info.index_storage_variants[dof_access_index][cell] ==
-            internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+            MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
               interleaved_contiguous)
         {
           AssertDimension(
@@ -2637,13 +2417,13 @@ namespace internal
                                                  [first_selected_component] *
               VectorizedArrayType::size();
 
-          if (fe_degree > 1 && integrate_gradients == true)
+          if (fe_degree > 1 && do_gradients == true)
             {
               // we know that the gradient weights for the Hermite case on the
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + side_];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
               const unsigned int *index_array =
@@ -2655,25 +2435,16 @@ namespace internal
                   AssertIndexRange(ind1, data.dofs_per_component_on_cell);
                   AssertIndexRange(ind2, data.dofs_per_component_on_cell);
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      VectorizedArrayType val =
-                        temp1[i + 2 * comp * dofs_per_face] -
-                        grad_weight *
-                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      VectorizedArrayType grad =
-                        grad_weight *
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      do_vectorized_add(val,
-                                        dst_ptr + dof_index +
-                                          (ind1 +
-                                           comp * static_dofs_per_component) *
-                                            VectorizedArrayType::size());
-                      do_vectorized_add(grad,
-                                        dst_ptr + dof_index +
-                                          (ind2 +
-                                           comp * static_dofs_per_component) *
-                                            VectorizedArrayType::size());
-                    }
+                    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
@@ -2686,23 +2457,22 @@ namespace internal
                 {
                   const unsigned int ind = index_array[i];
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    do_vectorized_add(temp1[i + 2 * comp * dofs_per_face],
-                                      dst_ptr + dof_index +
-                                        (ind +
-                                         comp * static_dofs_per_component) *
-                                          VectorizedArrayType::size());
+                    function_1b(temp1[i + 2 * comp * dofs_per_face],
+                                global_vector_ptr + dof_index +
+                                  (ind + comp * static_dofs_per_component) *
+                                    VectorizedArrayType::size());
                 }
             }
         }
 
       // case 2: contiguous and interleaved indices with fixed stride
-      else if (((integrate_gradients == false &&
+      else if (((do_gradients == false &&
                  data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
-                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
                dof_info.index_storage_variants[dof_access_index][cell] ==
-                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                    interleaved_contiguous_strided)
         {
           AssertDimension(
@@ -2712,13 +2482,13 @@ namespace internal
             &dof_info
                .dof_indices_contiguous[dof_access_index]
                                       [cell * VectorizedArrayType::size()];
-          if (fe_degree > 1 && integrate_gradients == true)
+          if (fe_degree > 1 && do_gradients == true)
             {
               // we know that the gradient weights for the Hermite case on the
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + side_];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2731,33 +2501,24 @@ namespace internal
                   const unsigned int ind2 =
                     index_array[2 * i + 1] * VectorizedArrayType::size();
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      VectorizedArrayType val =
-                        temp1[i + 2 * comp * dofs_per_face] -
-                        grad_weight *
-                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      VectorizedArrayType grad =
-                        grad_weight *
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      do_vectorized_scatter_add(
-                        val,
-                        indices,
-                        dst_ptr + ind1 +
-                          comp * static_dofs_per_component *
-                            VectorizedArrayType::size() +
-                          dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component] *
-                            VectorizedArrayType::size());
-                      do_vectorized_scatter_add(
-                        grad,
-                        indices,
-                        dst_ptr + ind2 +
-                          comp * static_dofs_per_component *
-                            VectorizedArrayType::size() +
-                          dof_info.component_dof_indices_offset
-                              [active_fe_index][first_selected_component] *
-                            VectorizedArrayType::size());
-                    }
+                    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
@@ -2771,27 +2532,27 @@ namespace internal
                   const unsigned int ind =
                     index_array[i] * VectorizedArrayType::size();
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    do_vectorized_scatter_add(
+                    function_2b(
                       temp1[i + 2 * comp * dofs_per_face],
-                      indices,
-                      dst_ptr + ind +
+                      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());
+                          VectorizedArrayType::size(),
+                      indices);
                 }
             }
         }
 
       // case 3: contiguous and interleaved indices with mixed stride
-      else if (((integrate_gradients == false &&
+      else if (((do_gradients == false &&
                  data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
-                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
                dof_info.index_storage_variants[dof_access_index][cell] ==
-                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                    interleaved_contiguous_mixed_strides)
         {
           const unsigned int *strides =
@@ -2808,13 +2569,13 @@ namespace internal
           const unsigned int n_filled_lanes =
             dof_info.n_vectorization_lanes_filled[dof_access_index][cell];
 
-          if (fe_degree > 1 && integrate_gradients == true)
+          if (fe_degree > 1 && do_gradients == true)
             {
               // we know that the gradient weights for the Hermite case on the
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + side_];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
 
@@ -2840,18 +2601,23 @@ namespace internal
                           indices[v] + (comp * static_dofs_per_component +
                                         index_array[2 * i + 1]) *
                                          strides[v];
-                      VectorizedArrayType val =
-                        temp1[i + 2 * comp * dofs_per_face] -
-                        grad_weight *
-                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      VectorizedArrayType grad =
-                        grad_weight *
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      do_vectorized_scatter_add(val, ind1, dst_ptr);
-                      do_vectorized_scatter_add(grad, ind2, dst_ptr);
+                      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
                 {
+                  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)
@@ -2864,15 +2630,12 @@ namespace internal
                             indices[v] + (comp * static_dofs_per_component +
                                           index_array[2 * i + 1]) *
                                            strides[v];
-                          Number val =
-                            temp1[i + 2 * comp * dofs_per_face][v] -
-                            grad_weight[0] * temp1[i + dofs_per_face +
-                                                   2 * comp * dofs_per_face][v];
-                          Number grad =
-                            grad_weight[0] * temp1[i + dofs_per_face +
-                                                   2 * comp * dofs_per_face][v];
-                          dst_ptr[ind1] += val;
-                          dst_ptr[ind2] += grad;
+                          function_3a(temp1[i + 2 * comp * dofs_per_face][v],
+                                      temp1[i + dofs_per_face +
+                                            2 * comp * dofs_per_face][v],
+                                      global_vector_ptr[ind1],
+                                      global_vector_ptr[ind2],
+                                      grad_weight[0]);
                         }
                 }
             }
@@ -2894,11 +2657,17 @@ namespace internal
                           indices[v] +
                           (comp * static_dofs_per_component + index_array[i]) *
                             strides[v];
-                      do_vectorized_scatter_add(
-                        temp1[i + 2 * comp * dofs_per_face], ind, dst_ptr);
+                      function_2b(temp1[i + 2 * comp * dofs_per_face],
+                                  global_vector_ptr,
+                                  ind);
                     }
               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)
@@ -2907,21 +2676,21 @@ namespace internal
                             indices[v] + (comp * static_dofs_per_component +
                                           index_array[i]) *
                                            strides[v];
-                          dst_ptr[ind1] +=
-                            temp1[i + 2 * comp * dofs_per_face][v];
+                          function_3b(temp1[i + 2 * comp * dofs_per_face][v],
+                                      global_vector_ptr[ind1]);
                         }
                 }
             }
         }
 
       // case 4: contiguous indices without interleaving
-      else if (((integrate_gradients == false &&
+      else if (((do_gradients == false &&
                  data.data.front().nodal_at_cell_boundaries == true) ||
                 (data.element_type ==
-                   internal::MatrixFreeFunctions::tensor_symmetric_hermite &&
+                   MatrixFreeFunctions::tensor_symmetric_hermite &&
                  fe_degree > 1)) &&
                dof_info.index_storage_variants[dof_access_index][cell] ==
-                 internal::MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
+                 MatrixFreeFunctions::DoFInfo::IndexStorageVariants::
                    contiguous &&
                dof_info.n_vectorization_lanes_filled[dof_access_index][cell] ==
                  VectorizedArrayType::size())
@@ -2931,17 +2700,18 @@ namespace internal
                .dof_indices_contiguous[dof_access_index]
                                       [cell * VectorizedArrayType::size()];
 
-          if (integrate_gradients == true &&
+          if (do_gradients == true &&
               data.element_type ==
-                internal::MatrixFreeFunctions::tensor_symmetric_hermite)
+                MatrixFreeFunctions::tensor_symmetric_hermite)
             {
               // we know that the gradient weights for the Hermite case on the
               // right (side==1) are the negative from the value at the left
               // (side==0), so we only read out one of them.
               const VectorizedArrayType grad_weight =
-                data.data.front().shape_data_on_face[0][fe_degree + 2 - side];
+                data.data.front().shape_data_on_face[0][fe_degree + side_];
               AssertDimension(data.face_to_cell_index_hermite.size(1),
                               2 * dofs_per_face);
+
               const unsigned int *index_array =
                 &data.face_to_cell_index_hermite(face_no, 0);
               for (unsigned int i = 0; i < dofs_per_face; ++i)
@@ -2949,27 +2719,20 @@ namespace internal
                   const unsigned int ind1 = index_array[2 * i];
                   const unsigned int ind2 = index_array[2 * i + 1];
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    {
-                      VectorizedArrayType val =
-                        temp1[i + 2 * comp * dofs_per_face] -
-                        grad_weight *
-                          temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      VectorizedArrayType grad =
-                        grad_weight *
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face];
-                      do_vectorized_scatter_add(
-                        val,
-                        indices,
-                        dst_ptr + comp * static_dofs_per_component + ind1 +
-                          dof_info.component_dof_indices_offset
-                            [active_fe_index][first_selected_component]);
-                      do_vectorized_scatter_add(
-                        grad,
-                        indices,
-                        dst_ptr + comp * static_dofs_per_component + ind2 +
-                          dof_info.component_dof_indices_offset
-                            [active_fe_index][first_selected_component]);
-                    }
+                    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
@@ -2979,32 +2742,40 @@ namespace internal
               const unsigned int *index_array =
                 &data.face_to_cell_index_nodal(face_no, 0);
               for (unsigned int i = 0; i < dofs_per_face; ++i)
-                {
-                  const unsigned int ind = index_array[i];
-                  for (unsigned int comp = 0; comp < n_components; ++comp)
-                    do_vectorized_scatter_add(
-                      temp1[i + 2 * comp * dofs_per_face],
-                      indices,
-                      dst_ptr + comp * static_dofs_per_component + ind +
-                        dof_info.component_dof_indices_offset
-                          [active_fe_index][first_selected_component]);
-                }
+                for (unsigned int comp = 0; comp < n_components; ++comp)
+                  {
+                    const unsigned int ind = index_array[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);
+                  }
             }
         }
 
-      // case 5: default vector access, must be handled separately, just do
-      // the face-normal interpolation
+      // case 5: default vector access
       else
         {
-          FEFaceNormalEvaluationImpl<dim,
-                                     fe_degree,
-                                     n_components,
-                                     VectorizedArrayType>::
-            template interpolate<false, false>(
-              data, temp1, values_array, integrate_gradients, face_no);
+          function_5(temp1);
           return false;
         }
 
+      if (!integrate)
+        function_0(temp1, dofs_per_face);
+
+      if (!integrate && (face_orientation > 0))
+        adjust_for_face_orientation(face_orientation,
+                                    orientation_map,
+                                    false,
+                                    do_values,
+                                    do_gradients,
+                                    data.n_q_points_face,
+                                    scratch_data,
+                                    values_quad,
+                                    gradients_quad);
+
       return true;
     }
 
@@ -3051,7 +2822,7 @@ namespace internal
               }
         }
     }
-  };
+  }; // namespace internal
 
 
 

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.