Explicitly use the inverse orientation instead of hard-coding.

Partially reverts #15678.

While #15678 fixed the permutation bug, it didn't address the true cause of the
problem: since face orientations are computed in the "apply this permutation to
face 1 to get face 2" direction, QProjector should use inverse orientations.
Since 2d orientations are their own inverses this only shows up in 3d.

Whenever two faces abutt, the first face is always in the default orientation
and the second face's orientation is computed relative to that (as decribed
here). Hence, when we project quadrature points onto the first face they do not
need to reoriented. However, we need to apply the *reverse* permutation on the
second face so that they end up in the same positions as the first face. More
formally: most, but not all, orientations are their own inverses. In particular,
triangle orientations 3 and 5 are each-other's inverses.

This matches the notion of inverse orientation used in #16828 for hypercubes. In
the future, we should combine the hypercube and non-hypercube implementations to
avoid these kinds of inconsistencies.

Part of #14667.

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 5d87fd625d3b4c28de53723be8c0c92323eefb61..e9d092d14cc6f2489c549f6761509ddb5f3504dc 100644 (file)
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -2942,9 +2942,9 @@ ReferenceCell::permute_by_combined_orientation(
             case 1:
               return {vertices[0], vertices[1], vertices[2]};
             case 3:
-              return {vertices[1], vertices[2], vertices[0]};
-            case 5:
               return {vertices[2], vertices[0], vertices[1]};
+            case 5:
+              return {vertices[1], vertices[2], vertices[0]};
             case 0:
               return {vertices[0], vertices[2], vertices[1]};
             case 2:

diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc
index 80ca05d5fab66620b7b09a4f0b11388437cb5f25..e7a6f2acc7007311c6ad95aec54373a9f3f4d26b 100644 (file)
--- a/source/base/qprojector.cc
+++ b/source/base/qprojector.cc
@@ -911,6 +911,8 @@ QProjector<3>::project_to_all_faces(const ReferenceCell      &reference_cell,
     // loop over all faces (triangles) ...
     for (unsigned int face_no = 0; face_no < faces.size(); ++face_no)
       {
+        const ReferenceCell face_reference_cell =
+          reference_cell.face_reference_cell(face_no);
         // We will use linear polynomials to map the reference quadrature
         // points correctly to on faces. There are as many linear shape
         // functions as there are vertices in the face.
@@ -929,9 +931,30 @@ QProjector<3>::project_to_all_faces(const ReferenceCell      &reference_cell,
           {
             const auto &face = faces[face_no];
 
+            // The goal of this function is to compute identical sets of
+            // quadrature points on the common face of two abutting cells. Our
+            // orientation convention is that, given such a pair of abutting
+            // cells:
+            //
+            // 1. The shared face, from the perspective of the first cell, is
+            //    in the default orientation.
+            // 2. The shared face, from the perspective of the second cell, has
+            //    its orientation computed relative to the first cell: i.e.,
+            //    'orientation' is the vertex permutation applied to the first
+            //    cell's face to get the second cell's face.
+            //
+            // The first case is trivial since points do not need to be
+            // oriented. However, in the second case, we need to use the
+            // *reverse* of the stored orientation (i.e., the permutation
+            // applied to the second cell's face which yields the first cell's
+            // face) so that we get identical quadrature points.
+            //
+            // For more information see connectivity.h.
             const boost::container::small_vector<Point<3>, 8> support_points =
-              reference_cell.face_reference_cell(face_no)
-                .permute_by_combined_orientation<Point<3>>(face, orientation);
+              face_reference_cell.permute_by_combined_orientation<Point<3>>(
+                face,
+                face_reference_cell.get_inverse_combined_orientation(
+                  orientation));
 
             // the quadrature rule to be projected ...
             const auto &sub_quadrature_points =

diff --git a/tests/simplex/orientation_02.cc b/tests/simplex/orientation_02.cc
index 74ab98587c6231e4f2b63c6185bd35f56c47db8b..54d1fa7148bec2a108584192dd44464521bf3325 100644 (file)
--- a/tests/simplex/orientation_02.cc
+++ b/tests/simplex/orientation_02.cc
@@ -23,7 +23,7 @@
 #include "../tests.h"
 
 void
-test(const unsigned int orientation)
+test(const unsigned char orientation)
 {
   const unsigned int face_no = 0;
 
@@ -43,11 +43,11 @@ test(const unsigned int orientation)
   }
 
   {
-    const auto &face = dummy.begin()->face(face_no);
+    const auto &face                = dummy.begin()->face(face_no);
+    const auto  face_reference_cell = face->reference_cell();
     const auto  permuted =
-      ReferenceCell(ReferenceCells::Triangle)
-        .permute_according_orientation(std::array<unsigned int, 3>{{0, 1, 2}},
-                                       orientation);
+      ReferenceCells::Triangle.permute_according_orientation(
+        std::array<unsigned int, 3>{{0, 1, 2}}, orientation);
 
     auto direction =
       cross_product_3d(vertices[permuted[1]] - vertices[permuted[0]],
@@ -57,13 +57,10 @@ test(const unsigned int orientation)
     vertices.push_back(face->center() + direction);
 
     CellData<3> cell;
-    cell.vertices.resize(4);
-
-    cell.vertices[permuted[0]] = face->vertex_index(0);
-    cell.vertices[permuted[1]] = face->vertex_index(1);
-    cell.vertices[permuted[2]] = face->vertex_index(2);
-    cell.vertices[3]           = 4;
-
+    cell.vertices    = {face->vertex_index(permuted[0]),
+                        face->vertex_index(permuted[1]),
+                        face->vertex_index(permuted[2]),
+                        4};
     cell.material_id = 1;
     cells.push_back(cell);
   }
@@ -74,7 +71,7 @@ test(const unsigned int orientation)
   cell++;
 
   // check orientation
-  deallog << "face orientation: " << orientation << ' '
+  deallog << "face orientation: " << int(orientation) << ' '
           << int(cell->combined_face_orientation(0)) << ' ' << std::endl;
 
   // check vertices