From bed40b1a149cc3439705865ef31d9cc9926eff1b Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 8 Feb 2025 20:51:39 -0500 Subject: [PATCH] QProjector: add a new project_to_face() overload. This one will replace the other ones in the next step. --- include/deal.II/base/qprojector.h | 11 +++++ source/base/qprojector.cc | 76 +++++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+) diff --git a/include/deal.II/base/qprojector.h b/include/deal.II/base/qprojector.h index 400c3ebf20..37e16f9d4d 100644 --- a/include/deal.II/base/qprojector.h +++ b/include/deal.II/base/qprojector.h @@ -122,6 +122,17 @@ public: const bool face_flip, const bool face_rotation); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on face face_no. For further details, see + * the general doc for this class. + */ + static Quadrature + project_to_face(const ReferenceCell &reference_cell, + const SubQuadrature &quadrature, + const unsigned int face_no, + const types::geometric_orientation orientation); + /** * Compute the quadrature points on the cell if the given quadrature formula * is used on face face_no, subface number subface_no diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 6c34ab2d45..a6a897d566 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -462,6 +462,82 @@ QProjector<3>::project_to_oriented_face(const ReferenceCell &reference_cell, +template +Quadrature +QProjector::project_to_face(const ReferenceCell &reference_cell, + const Quadrature &quadrature, + const unsigned int face_no, + const types::geometric_orientation orientation) +{ + AssertIndexRange(face_no, reference_cell.n_faces()); + AssertIndexRange(orientation, reference_cell.n_face_orientations(face_no)); + AssertDimension(reference_cell.get_dimension(), dim); + + std::vector> q_points; + std::vector q_weights = quadrature.get_weights(); + q_points.reserve(quadrature.size()); + if constexpr (dim == 1) + { + AssertDimension(quadrature.size(), 1); + q_points.emplace_back(static_cast(face_no)); + } + else if constexpr (dim == 2) + { + if (reference_cell == ReferenceCells::Triangle) + // use linear polynomial to map the reference quadrature points + // correctly on faces, i.e., BarycentricPolynomials<1>(1) + for (unsigned int p = 0; p < quadrature.size(); ++p) + { + if (face_no == 0) + q_points.emplace_back(quadrature.point(p)[0], 0); + else if (face_no == 1) + q_points.emplace_back(1.0 - quadrature.point(p)[0], + quadrature.point(p)[0]); + else if (face_no == 2) + q_points.emplace_back(0, 1.0 - quadrature.point(p)[0]); + else + DEAL_II_ASSERT_UNREACHABLE(); + } + else if (reference_cell == ReferenceCells::Quadrilateral) + for (unsigned int p = 0; p < quadrature.size(); ++p) + { + if (face_no == 0) + q_points.emplace_back(0, quadrature.point(p)[0]); + else if (face_no == 1) + q_points.emplace_back(1, quadrature.point(p)[0]); + else if (face_no == 2) + q_points.emplace_back(quadrature.point(p)[0], 0); + else if (face_no == 3) + q_points.emplace_back(quadrature.point(p)[0], 1); + else + DEAL_II_ASSERT_UNREACHABLE(); + } + else + DEAL_II_ASSERT_UNREACHABLE(); + + if (orientation == numbers::reverse_line_orientation) + { + std::reverse(q_points.begin(), q_points.end()); + std::reverse(q_weights.begin(), q_weights.end()); + } + } + else if constexpr (dim == 3) + { + Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); + const Quadrature<2> mutation = + internal::QProjector::mutate_quadrature(quadrature, orientation); + return QProjector<3>::project_to_face(reference_cell, mutation, face_no); + } + else + { + DEAL_II_ASSERT_UNREACHABLE(); + } + + return Quadrature(std::move(q_points), std::move(q_weights)); +} + + + template <> void QProjector<1>::project_to_subface(const ReferenceCell &reference_cell, -- 2.39.5