From: David Wells Date: Sat, 20 Apr 2024 21:57:44 +0000 (-0400) Subject: Explicitly use the inverse orientation instead of hard-coding. X-Git-Tag: v9.6.0-rc1~351^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F16916%2Fhead;p=dealii.git 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 5d87fd625d..e9d092d14c 100644 --- 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 80ca05d5fa..e7a6f2acc7 100644 --- 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, 8> support_points = - reference_cell.face_reference_cell(face_no) - .permute_by_combined_orientation>(face, orientation); + face_reference_cell.permute_by_combined_orientation>( + 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 74ab98587c..54d1fa7148 100644 --- 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{{0, 1, 2}}, - orientation); + ReferenceCells::Triangle.permute_according_orientation( + std::array{{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