From 8446074c6e913d078cfeadefdd4b7921e41d49a4 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sun, 13 Apr 2025 13:26:16 -0400 Subject: [PATCH] QProjector: refactor project_to_all_faces(). --- source/base/qprojector.cc | 107 ++++++++------------------------------ 1 file changed, 22 insertions(+), 85 deletions(-) diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 29d68eace3..66fc46bd82 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -770,95 +770,32 @@ QProjector::project_to_all_faces( const ReferenceCell &reference_cell, const hp::QCollection &quadrature) { - const auto process = [&](const std::vector>> &faces) { - // new (projected) quadrature points and weights - std::vector> points; - std::vector weights; + std::vector> points; + std::vector weights; - // 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); - - // ... and over all possible orientations - for (types::geometric_orientation orientation = 0; - orientation < reference_cell.n_face_orientations(face_no); - ++orientation) - { - 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 = - 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 = - quadrature[quadrature.size() == 1 ? 0 : face_no].get_points(); - const auto &sub_quadrature_weights = - quadrature[quadrature.size() == 1 ? 0 : face_no].get_weights(); - - // loop over all quadrature points - for (unsigned int j = 0; j < sub_quadrature_points.size(); ++j) - { - Point mapped_point; - - // map reference quadrature point - for (const unsigned int i : - face_reference_cell.vertex_indices()) - mapped_point += support_points[i] * - face_reference_cell.d_linear_shape_function( - sub_quadrature_points[j], i); - - points.push_back(mapped_point); - - // rescale quadrature weights so that the sum of the weights on - // each face equals the measure of that face. - const double scaling = reference_cell.face_measure(face_no) / - face_reference_cell.volume(); - weights.push_back(sub_quadrature_weights[j] * scaling); - } - } - } - - // construct new quadrature rule - return Quadrature(std::move(points), std::move(weights)); - }; - - std::vector>> face_vertex_locations( - reference_cell.n_faces()); - for (const unsigned int f : reference_cell.face_indices()) + for (const unsigned int face_no : reference_cell.face_indices()) { - face_vertex_locations[f].resize( - reference_cell.face_reference_cell(f).n_vertices()); - for (const unsigned int v : - reference_cell.face_reference_cell(f).vertex_indices()) - face_vertex_locations[f][v] = - reference_cell.face_vertex_location(f, v); + const ReferenceCell face_reference_cell = + reference_cell.face_reference_cell(face_no); + std::vector> face_vertices(face_reference_cell.n_vertices()); + for (const unsigned int vertex_no : face_reference_cell.vertex_indices()) + face_vertices[vertex_no] = + reference_cell.face_vertex_location(face_no, vertex_no); + + for (types::geometric_orientation combined_orientation = 0; + combined_orientation < reference_cell.n_face_orientations(face_no); + ++combined_orientation) + internal::QProjector::append_subobject_rule( + face_reference_cell, + quadrature[quadrature.size() == 1 ? 0 : face_no], + face_vertices, + reference_cell.face_measure(face_no), + combined_orientation, + points, + weights); } - return process(face_vertex_locations); + return Quadrature(std::move(points), std::move(weights)); } -- 2.39.5