From: Martin Kronbichler <martin.kronbichler@it.uu.se>
Date: Mon, 20 Dec 2021 09:20:13 +0000 (+0100)
Subject: Bug fix of MF face eval for hanging nodes and non-standard orientation
X-Git-Tag: v9.4.0-rc1~722^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e1dd4e9a25420bbf3c545611c883a422973e5ead;p=dealii.git

Bug fix of MF face eval for hanging nodes and non-standard orientation
---

diff --git a/include/deal.II/matrix_free/evaluation_kernels.h b/include/deal.II/matrix_free/evaluation_kernels.h
index d4101092fa..d26e9ce2b4 100644
--- a/include/deal.II/matrix_free/evaluation_kernels.h
+++ b/include/deal.II/matrix_free/evaluation_kernels.h
@@ -3783,19 +3783,18 @@ namespace internal
                 }
             }
           else if (fe_eval.get_face_orientation() != 0)
-            for (unsigned int c = 0; c < n_components; ++c)
-              adjust_for_face_orientation(
-                dim,
-                n_components,
-                evaluation_flag,
-                &fe_eval.get_shape_info().face_orientations_quad(
-                  fe_eval.get_face_orientation(), 0),
-                false,
-                Utilities::pow(n_q_points_1d, dim - 1),
-                temp,
-                fe_eval.begin_values(),
-                fe_eval.begin_gradients(),
-                fe_eval.begin_hessians());
+            adjust_for_face_orientation(
+              dim,
+              n_components,
+              evaluation_flag,
+              &fe_eval.get_shape_info().face_orientations_quad(
+                fe_eval.get_face_orientation(), 0),
+              false,
+              Utilities::pow(n_q_points_1d, dim - 1),
+              temp,
+              fe_eval.begin_values(),
+              fe_eval.begin_gradients(),
+              fe_eval.begin_hessians());
         }
 
       return false;
diff --git a/include/deal.II/matrix_free/face_setup_internal.h b/include/deal.II/matrix_free/face_setup_internal.h
index bcfd049961..bdcc6516e6 100644
--- a/include/deal.II/matrix_free/face_setup_internal.h
+++ b/include/deal.II/matrix_free/face_setup_internal.h
@@ -28,6 +28,7 @@
 #include <deal.II/grid/tria_accessor.h>
 
 #include <deal.II/matrix_free/face_info.h>
+#include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/task_info.h>
 
 #include <fstream>
@@ -1041,6 +1042,20 @@ namespace internal
         }
       else
         info.face_orientation = exterior_face_orientation;
+
+      // make sure to select correct subface index in case of non-standard
+      // orientation of the coarser neighbor face
+      if (cell->level() > neighbor->level() && exterior_face_orientation > 0)
+        {
+          const Table<2, unsigned int> orientation =
+            ShapeInfo<double>::compute_orientation_table(2);
+          const std::array<unsigned int, 8> inverted_orientations{
+            {0, 1, 2, 3, 6, 5, 4, 7}};
+          info.subface_index =
+            orientation[inverted_orientations[exterior_face_orientation]]
+                       [info.subface_index];
+        }
+
       return info;
     }
 
diff --git a/include/deal.II/matrix_free/shape_info.h b/include/deal.II/matrix_free/shape_info.h
index d7269738da..1187d0708c 100644
--- a/include/deal.II/matrix_free/shape_info.h
+++ b/include/deal.II/matrix_free/shape_info.h
@@ -375,6 +375,13 @@ namespace internal
       static bool
       is_supported(const FiniteElement<dim, spacedim> &fe);
 
+      /**
+       * Compute a table with numbers of re-orientation for all versions of
+       * face flips, orientation, and rotation (relating only to 3D elements).
+       */
+      static Table<2, unsigned int>
+      compute_orientation_table(const unsigned int n_points_per_dim);
+
       /**
        * Return data of univariate shape functions which defines the
        * dimension @p dimension of tensor product shape functions
diff --git a/include/deal.II/matrix_free/shape_info.templates.h b/include/deal.II/matrix_free/shape_info.templates.h
index fd13f5fb0f..9f1a084a81 100644
--- a/include/deal.II/matrix_free/shape_info.templates.h
+++ b/include/deal.II/matrix_free/shape_info.templates.h
@@ -854,41 +854,8 @@ namespace internal
           // (similar to MappingInfoStorage::QuadratureDescriptor::initialize)
           if (dim == 3)
             {
-              const auto compute_orientations =
-                [](const unsigned int      n,
-                   Table<2, unsigned int> &face_orientations) {
-                  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;
-                      }
-                };
-              compute_orientations(fe_degree + 1, face_orientations_dofs);
-              compute_orientations(n_q_points_1d, face_orientations_quad);
+              face_orientations_dofs = compute_orientation_table(fe_degree + 1);
+              face_orientations_quad = compute_orientation_table(n_q_points_1d);
             }
           else
             {
@@ -1163,6 +1130,36 @@ namespace internal
 
 
 
+    template <typename Number>
+    Table<2, unsigned int>
+    ShapeInfo<Number>::compute_orientation_table(const unsigned int n)
+    {
+      Table<2, unsigned int> face_orientations(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;
+          }
+      return face_orientations;
+    }
+
+
+
     template <typename Number>
     std::size_t
     ShapeInfo<Number>::memory_consumption() const