Merge gather and adjust_for_face_orientation

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index 1a101a63d01739c79a71a2c5b8bc5c8dd2cab6d8..0cadb08e874b1d6b323b4636073cfd6422377c22 100644 (file)
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -2355,7 +2355,9 @@ namespace internal
     {
       (void)subface_index;
 
-      if (integrate && (face_orientation > 0))
+      if (integrate &&
+          (face_orientation > 0 &&
+           subface_index < GeometryInfo<dim>::max_children_per_cell))
         adjust_for_face_orientation(face_orientation,
                                     orientation_map,
                                     true,
@@ -2396,6 +2398,18 @@ namespace internal
       if (integrate)
         function_0(temp1, dofs_per_face);
 
+      const unsigned int  dummy = 0;
+      const unsigned int *orientation =
+        (data.data.front().nodal_at_cell_boundaries == true) ?
+          &data.face_orientations[face_orientation][0] :
+          &dummy;
+      const auto reorientate = [&](const unsigned int i) {
+        return (dim < 3 || face_orientation == 0 ||
+                subface_index < GeometryInfo<dim>::max_children_per_cell) ?
+                 i :
+                 orientation[i];
+      };
+
       // case 1: contiguous and interleaved indices
       if (((do_gradients == false &&
             data.data.front().nodal_at_cell_boundaries == true) ||
@@ -2434,10 +2448,11 @@ namespace internal
                   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);
+                  const unsigned int i_ = reorientate(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],
+                      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(),
@@ -2455,9 +2470,10 @@ namespace internal
                 &data.face_to_cell_index_nodal(face_no, 0);
               for (unsigned int i = 0; i < dofs_per_face; ++i)
                 {
+                  const unsigned int i_  = reorientate(i);
                   const unsigned int ind = index_array[i];
                   for (unsigned int comp = 0; comp < n_components; ++comp)
-                    function_1b(temp1[i + 2 * comp * dofs_per_face],
+                    function_1b(temp1[i_ + 2 * comp * dofs_per_face],
                                 global_vector_ptr + dof_index +
                                   (ind + comp * static_dofs_per_component) *
                                     VectorizedArrayType::size());
@@ -2496,14 +2512,16 @@ namespace internal
                 &data.face_to_cell_index_hermite(face_no, 0);
               for (unsigned int i = 0; i < dofs_per_face; ++i)
                 {
+                  const unsigned int i_ = reorientate(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)
                     function_2a(
-                      temp1[i + 2 * comp * dofs_per_face],
-                      temp1[i + dofs_per_face + 2 * comp * dofs_per_face],
+                      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() +
@@ -2529,11 +2547,13 @@ namespace internal
                 &data.face_to_cell_index_nodal(face_no, 0);
               for (unsigned int i = 0; i < dofs_per_face; ++i)
                 {
+                  const unsigned int i_ = reorientate(i);
+
                   const unsigned int ind =
                     index_array[i] * VectorizedArrayType::size();
                   for (unsigned int comp = 0; comp < n_components; ++comp)
                     function_2b(
-                      temp1[i + 2 * comp * dofs_per_face],
+                      temp1[i_ + 2 * comp * dofs_per_face],
                       global_vector_ptr + ind +
                         comp * static_dofs_per_component *
                           VectorizedArrayType::size() +
@@ -2585,7 +2605,8 @@ namespace internal
                 for (unsigned int comp = 0; comp < n_components; ++comp)
                   for (unsigned int i = 0; i < dofs_per_face; ++i)
                     {
-                      unsigned int ind1[VectorizedArrayType::size()];
+                      const unsigned int i_ = reorientate(i);
+                      unsigned int       ind1[VectorizedArrayType::size()];
                       DEAL_II_OPENMP_SIMD_PRAGMA
                       for (unsigned int v = 0; v < VectorizedArrayType::size();
                            ++v)
@@ -2602,8 +2623,8 @@ namespace internal
                                         index_array[2 * i + 1]) *
                                          strides[v];
                       function_2a(
-                        temp1[i + 2 * comp * dofs_per_face],
-                        temp1[i + dofs_per_face + 2 * comp * dofs_per_face],
+                        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,
@@ -2622,6 +2643,7 @@ namespace internal
                     for (unsigned int comp = 0; comp < n_components; ++comp)
                       for (unsigned int i = 0; i < dofs_per_face; ++i)
                         {
+                          const unsigned int i_ = reorientate(i);
                           const unsigned int ind1 =
                             indices[v] + (comp * static_dofs_per_component +
                                           index_array[2 * i]) *
@@ -2630,8 +2652,8 @@ namespace internal
                             indices[v] + (comp * static_dofs_per_component +
                                           index_array[2 * i + 1]) *
                                            strides[v];
-                          function_3a(temp1[i + 2 * comp * dofs_per_face][v],
-                                      temp1[i + dofs_per_face +
+                          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],
@@ -2657,7 +2679,8 @@ namespace internal
                           indices[v] +
                           (comp * static_dofs_per_component + index_array[i]) *
                             strides[v];
-                      function_2b(temp1[i + 2 * comp * dofs_per_face],
+                      const unsigned int i_ = reorientate(i);
+                      function_2b(temp1[i_ + 2 * comp * dofs_per_face],
                                   global_vector_ptr,
                                   ind);
                     }
@@ -2676,7 +2699,8 @@ namespace internal
                             indices[v] + (comp * static_dofs_per_component +
                                           index_array[i]) *
                                            strides[v];
-                          function_3b(temp1[i + 2 * comp * dofs_per_face][v],
+                          const unsigned int i_ = reorientate(i);
+                          function_3b(temp1[i_ + 2 * comp * dofs_per_face][v],
                                       global_vector_ptr[ind1]);
                         }
                 }
@@ -2718,10 +2742,12 @@ namespace internal
                 {
                   const unsigned int ind1 = index_array[2 * i];
                   const unsigned int ind2 = index_array[2 * i + 1];
+                  const unsigned int i_   = reorientate(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],
+                      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
@@ -2745,7 +2771,9 @@ namespace internal
                 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],
+                    const unsigned int i_  = reorientate(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
@@ -2765,7 +2793,9 @@ namespace internal
       if (!integrate)
         function_0(temp1, dofs_per_face);
 
-      if (!integrate && (face_orientation > 0))
+      if (!integrate &&
+          (face_orientation > 0 &&
+           subface_index < GeometryInfo<dim>::max_children_per_cell))
         adjust_for_face_orientation(face_orientation,
                                     orientation_map,
                                     false,

diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index ccfe4902112c3966191f5eeab71b34e20589bc14..f20abbc98dcca60152491b0f971a1a71c6210af2 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -467,6 +467,14 @@ namespace internal
        */
       dealii::Table<2, unsigned int> face_to_cell_index_hermite;
 
+      /**
+       * For degrees on faces, the basis functions are not
+       * in the correct order if a face is not in the standard orientation
+       * to a given element. This data structure is used to re-order the
+       * basis functions to represent the correct order.
+       */
+      dealii::Table<2, unsigned int> face_orientations;
+
     private:
       /**
        * Check whether we have symmetries in the shape values. In that case,

diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index 4bced728889da95bc01fab2ccf68d286f25185a0..8b12453bab636ba0d203535eb495ace17d15275f 100644 (file)
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -481,6 +481,46 @@ namespace internal
                       face_to_cell_index_nodal(f, l) = ind;
                     }
             }
+
+          // face orientation for faces in 3D
+          // (similar to MappingInfoStorage::QuadratureDescriptor::initialize)
+          if (dim == 3)
+            {
+              const unsigned int n = fe_degree + 1;
+              face_orientations.reinit(8, n * n);
+              for (unsigned int j = 0, i = 0; j < n; ++j)
+                for (unsigned int k = 0; k < n; ++k, ++i)
+                  {
+                    // face_orientation=true,  face_flip=false,
+                    // face_rotation=false
+                    face_orientations[0][i] = i;
+                    // face_orientation=false, face_flip=false,
+                    // face_rotation=false
+                    face_orientations[1][i] = j + k * n;
+                    // face_orientation=true,  face_flip=true,
+                    // face_rotation=false
+                    face_orientations[2][i] = (n - 1 - k) + (n - 1 - j) * n;
+                    // face_orientation=false, face_flip=true,
+                    // face_rotation=false
+                    face_orientations[3][i] = (n - 1 - j) + (n - 1 - k) * n;
+                    // face_orientation=true,  face_flip=false,
+                    // face_rotation=true
+                    face_orientations[4][i] = j + (n - 1 - k) * n;
+                    // face_orientation=false, face_flip=false,
+                    // face_rotation=true
+                    face_orientations[5][i] = k + (n - 1 - j) * n;
+                    // face_orientation=true,  face_flip=true,
+                    // face_rotation=true
+                    face_orientations[6][i] = (n - 1 - j) + k * n;
+                    // face_orientation=false, face_flip=true,
+                    // face_rotation=true
+                    face_orientations[7][i] = (n - 1 - k) + j * n;
+                  }
+            }
+          else
+            {
+              face_orientations.reinit(1, 1);
+            }
         }
 
       if (element_type == tensor_symmetric_hermite)