From c95b4e227dbdc5027858b81121a163cbe1300738 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 15 Feb 2025 09:07:41 -0500 Subject: [PATCH] QProjector: add another project_to_subface() overload. --- include/deal.II/base/qprojector.h | 17 ++++ source/base/qprojector.cc | 139 ++++++++++++++++++++++++------ 2 files changed, 129 insertions(+), 27 deletions(-) diff --git a/include/deal.II/base/qprojector.h b/include/deal.II/base/qprojector.h index 07643e71ab..7aa720229e 100644 --- a/include/deal.II/base/qprojector.h +++ b/include/deal.II/base/qprojector.h @@ -209,6 +209,23 @@ public: const bool face_rotation, const internal::SubfaceCase ref_case); + /** + * Compute the cell quadrature formula corresponding to using + * quadrature on subface subface_no of face + * face_no with RefinementCase ref_case. The last + * argument is only used in 3d. + * + * @note Only the points are transformed. The quadrature weights are the + * same as those of the original rule. + */ + static Quadrature + project_to_subface(const ReferenceCell &reference_cell, + const SubQuadrature &quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const types::geometric_orientation orientation, + const RefinementCase &ref_case); + /** * Take a collection of face quadrature formulas and generate a cell * quadrature formula from it where the quadrature points of the given diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index 5a8b0ee746..1af849c4b4 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -658,37 +658,122 @@ QProjector::project_to_oriented_subface( -template <> -Quadrature<3> -QProjector<3>::project_to_oriented_subface( - const ReferenceCell &reference_cell, - const Quadrature<2> &quadrature, - const unsigned int face_no, - const unsigned int subface_no, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, - const internal::SubfaceCase<3> ref_case) +template +Quadrature +QProjector::project_to_subface( + const ReferenceCell &reference_cell, + const SubQuadrature &quadrature, + const unsigned int face_no, + const unsigned int subface_no, + const types::geometric_orientation orientation, + const RefinementCase &ref_case) { - Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); + AssertIndexRange(face_no, reference_cell.n_faces()); + // TODO once the default orientation is zero we can remove this special case + AssertIndexRange(dim == 1 ? 1 - orientation : orientation, + reference_cell.n_face_orientations(face_no)); + AssertDimension(reference_cell.get_dimension(), dim); + AssertIndexRange(subface_no, + reference_cell.face_reference_cell(face_no) + .template n_children(ref_case)); - const Quadrature<2> mutation = internal::QProjector::mutate_quadrature( - quadrature, - internal::combined_face_orientation(face_orientation, - face_rotation, - face_flip)); + std::vector> q_points; + std::vector q_weights = quadrature.get_weights(); + q_points.reserve(quadrature.size()); - const std::pair> - final_subface_no_and_ref_case = - internal::QProjector::select_subface_no_and_refinement_case( - subface_no, face_orientation, face_flip, face_rotation, ref_case); + 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) + { + if (subface_no == 0) + q_points.emplace_back(quadrature.point(p)[0] / 2, 0); + else + q_points.emplace_back(0.5 + quadrature.point(p)[0] / 2, 0); + } + else if (face_no == 1) + { + if (subface_no == 0) + q_points.emplace_back(1 - quadrature.point(p)[0] / 2, + quadrature.point(p)[0] / 2); + else + q_points.emplace_back(0.5 - quadrature.point(p)[0] / 2, + 0.5 + quadrature.point(p)[0] / 2); + } + else if (face_no == 2) + { + if (subface_no == 0) + q_points.emplace_back(0, 1 - quadrature.point(p)[0] / 2); + else + q_points.emplace_back(0, 0.5 - quadrature.point(p)[0] / 2); + } + else + DEAL_II_ASSERT_UNREACHABLE(); + } + else if (reference_cell == ReferenceCells::Quadrilateral) + for (unsigned int p = 0; p < quadrature.size(); ++p) + { + if (face_no == 0) + { + if (subface_no == 0) + q_points.emplace_back(0, quadrature.point(p)[0] / 2); + else + q_points.emplace_back(0, quadrature.point(p)[0] / 2 + 0.5); + } + else if (face_no == 1) + { + if (subface_no == 0) + q_points.emplace_back(1, quadrature.point(p)[0] / 2); + else + q_points.emplace_back(1, quadrature.point(p)[0] / 2 + 0.5); + } + else if (face_no == 2) + { + if (subface_no == 0) + q_points.emplace_back(quadrature.point(p)[0] / 2, 0); + else + q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 0); + } + else if (face_no == 3) + { + if (subface_no == 0) + q_points.emplace_back(quadrature.point(p)[0] / 2, 1); + else + q_points.emplace_back(quadrature.point(p)[0] / 2 + 0.5, 1); + } + else + DEAL_II_ASSERT_UNREACHABLE(); + } + else + DEAL_II_ASSERT_UNREACHABLE(); - return QProjector<3>::project_to_subface( - reference_cell, - mutation, - face_no, - final_subface_no_and_ref_case.first, - final_subface_no_and_ref_case.second); + 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()); + internal::QProjector::project_to_hex_face_and_append( + quadrature.get_points(), face_no, q_points, ref_case, subface_no); + } + else + { + DEAL_II_ASSERT_UNREACHABLE(); + } + + return Quadrature(std::move(q_points), std::move(q_weights)); } -- 2.39.5