From: David Wells Date: Sat, 7 Jun 2025 19:11:21 +0000 (-0400) Subject: Move some more code from QProjector to ReferenceCell. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9f9e25e30e6a4c22e7f248dd3c2bd8fcf9beacc3;p=dealii.git Move some more code from QProjector to ReferenceCell. We should implement all such things in ReferenceCell. We only need this for subfaces in 3d (as anisotropically refined faces in 2d are equivalent to isotropically refined faces). --- diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h index 6a45af7aa4..9546a4930f 100644 --- a/include/deal.II/grid/reference_cell.h +++ b/include/deal.II/grid/reference_cell.h @@ -410,6 +410,18 @@ public: std_cxx20::ranges::iota_view isotropic_child_indices() const; + /** + * Convert an internal::SubfaceCase into its equivalent RefinementCase. For + * example, SubfacePossibilities<3>::Possibilities::y1x2x is is geometrically + * equivalent to RefinementPossibilities<3>::cut_xy. + */ + template + std::pair> + equivalent_refinement_case( + const types::geometric_orientation combined_face_orientation, + const internal::SubfaceCase subface_case, + const unsigned int subface_no) const; + /** * Return the reference-cell type of face @p face_no of the current * object. For example, if the current object is diff --git a/source/base/qprojector.cc b/source/base/qprojector.cc index dd2bfdbbe7..11c17f0465 100644 --- a/source/base/qprojector.cc +++ b/source/base/qprojector.cc @@ -188,146 +188,6 @@ namespace internal - std::pair> - select_subface_no_and_refinement_case( - const unsigned int subface_no, - const bool face_orientation, - const bool face_flip, - const bool face_rotation, - const internal::SubfaceCase<3> ref_case) - { - constexpr int dim = 3; - // for each subface of a given FaceRefineCase - // there is a corresponding equivalent - // subface number of one of the "standard" - // RefineCases (cut_x, cut_y, cut_xy). Map - // the given values to those equivalent - // ones. - - // first, define an invalid number - static const unsigned int e = numbers::invalid_unsigned_int; - - static const RefinementCase - equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] - [GeometryInfo<3>::max_children_per_face] = { - // case_none. there should be only - // invalid values here. However, as - // this function is also called (in - // tests) for cells which have no - // refined faces, use isotropic - // refinement instead - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_x - {RefinementCase::cut_x, - RefinementCase::cut_x, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_x1y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_x, - RefinementCase::no_refinement}, - // case_x2y - {RefinementCase::cut_x, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_x1y2y - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_y - {RefinementCase::cut_y, - RefinementCase::cut_y, - RefinementCase::no_refinement, - RefinementCase::no_refinement}, - // case_y1x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_y, - RefinementCase::no_refinement}, - // case_y2x - {RefinementCase::cut_y, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::no_refinement}, - // case_y1x2x - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}, - // case_xy (case_isotropic) - {RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy, - RefinementCase::cut_xy}}; - - static const unsigned int - equivalent_subface_number[internal::SubfaceCase::case_isotropic + - 1][GeometryInfo<3>::max_children_per_face] = - {// case_none, see above - {0, 1, 2, 3}, - // case_x - {0, 1, e, e}, - // case_x1y - {0, 2, 1, e}, - // case_x2y - {0, 1, 3, e}, - // case_x1y2y - {0, 2, 1, 3}, - // case_y - {0, 1, e, e}, - // case_y1x - {0, 1, 1, e}, - // case_y2x - {0, 2, 3, e}, - // case_y1x2x - {0, 1, 2, 3}, - // case_xy (case_isotropic) - {0, 1, 2, 3}}; - - // If face-orientation or face_rotation are - // non-standard, cut_x and cut_y have to be - // exchanged. - static const RefinementCase ref_case_permutation[4] = { - RefinementCase::no_refinement, - RefinementCase::cut_y, - RefinementCase::cut_x, - RefinementCase::cut_xy}; - - // set a corresponding (equivalent) - // RefineCase and subface number - const RefinementCase equ_ref_case = - equivalent_refine_case[ref_case][subface_no]; - const unsigned int equ_subface_no = - equivalent_subface_number[ref_case][subface_no]; - // make sure, that we got a valid subface and RefineCase - Assert(equ_ref_case != RefinementCase::no_refinement, - ExcInternalError()); - Assert(equ_subface_no != e, ExcInternalError()); - // now, finally respect non-standard faces - const RefinementCase final_ref_case = - (face_orientation == face_rotation ? - ref_case_permutation[equ_ref_case] : - equ_ref_case); - - const unsigned int final_subface_no = - GeometryInfo::child_cell_on_face(RefinementCase( - final_ref_case), - 4, - equ_subface_no, - face_orientation, - face_flip, - face_rotation, - equ_ref_case); - - return std::make_pair(final_subface_no, final_ref_case); - } - /** * Append the points and weights of a quadrature rule projected onto a * subobject (either face or subface) to the end of two vectors. These @@ -1187,9 +1047,6 @@ QProjector<3>::DataSetDescriptor::subface( { const unsigned int dim = 3; - const auto [face_orientation, face_rotation, face_flip] = - internal::split_face_orientation(combined_orientation); - Assert(reference_cell == ReferenceCells::Hexahedron, ExcNotImplemented()); (void)reference_cell; @@ -1228,14 +1085,13 @@ QProjector<3>::DataSetDescriptor::subface( 0 // cut_xy }; - 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); + const auto [final_subface_no, refinement_case] = + reference_cell.equivalent_refinement_case(combined_orientation, + ref_case, + subface_no); return ((face_no * total_subfaces_per_face + - ref_case_offset[final_subface_no_and_ref_case.second - 1] + - final_subface_no_and_ref_case.first) + + ref_case_offset[refinement_case - 1] + final_subface_no) + GeometryInfo::faces_per_cell * total_subfaces_per_face * combined_orientation) * n_quadrature_points; diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc index b1bf43b63e..264355627e 100644 --- a/source/grid/reference_cell.cc +++ b/source/grid/reference_cell.cc @@ -116,6 +116,146 @@ ReferenceCell::to_string() const +template +std::pair> +ReferenceCell::equivalent_refinement_case( + const types::geometric_orientation combined_face_orientation, + const internal::SubfaceCase subface_case, + const unsigned int subface_no) const +{ + if constexpr (dim == 3) + { + static const RefinementCase + equivalent_refine_case[internal::SubfaceCase::case_isotropic + 1] + [GeometryInfo<3>::max_children_per_face] = { + // case_none. there should be only + // invalid values here. However, as + // this function is also called (in + // tests) for cells which have no + // refined faces, use isotropic + // refinement instead + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_x + {RefinementCase::cut_x, + RefinementCase::cut_x, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_x1y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_x, + RefinementCase::no_refinement}, + // case_x2y + {RefinementCase::cut_x, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_x1y2y + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_y + {RefinementCase::cut_y, + RefinementCase::cut_y, + RefinementCase::no_refinement, + RefinementCase::no_refinement}, + // case_y1x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_y, + RefinementCase::no_refinement}, + // case_y2x + {RefinementCase::cut_y, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::no_refinement}, + // case_y1x2x + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}, + // case_xy (case_isotropic) + {RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy, + RefinementCase::cut_xy}}; + + constexpr unsigned int X = numbers::invalid_unsigned_int; + static const unsigned int + equivalent_subface_number[internal::SubfaceCase::case_isotropic + + 1][GeometryInfo<3>::max_children_per_face] = { + // case_none, see above + {0, 1, 2, 3}, + // case_x + {0, 1, X, X}, + // case_x1y + {0, 2, 1, X}, + // case_x2y + {0, 1, 3, X}, + // case_x1y2y + {0, 2, 1, 3}, + // case_y + {0, 1, X, X}, + // case_y1x + {0, 1, 1, X}, + // case_y2x + {0, 2, 3, X}, + // case_y1x2x + {0, 1, 2, 3}, + // case_xy (case_isotropic) + {0, 1, 2, 3}}; + + static const RefinementCase rotated_refinement_case[4] = { + RefinementCase::no_refinement, + RefinementCase::cut_y, + RefinementCase::cut_x, + RefinementCase::cut_xy}; + const auto [face_orientation, face_rotation, face_flip] = + internal::split_face_orientation(combined_face_orientation); + + const auto equivalent_refinement_case = + equivalent_refine_case[subface_case][subface_no]; + const unsigned int equivalent_subface_no = + equivalent_subface_number[subface_case][subface_no]; + // make sure, that we got a valid subface and RefineCase + Assert(equivalent_refinement_case != RefinementCase::no_refinement, + ExcInternalError()); + Assert(equivalent_subface_no != X, ExcInternalError()); + // now, finally respect non-standard faces + const RefinementCase final_refinement_case = + (face_orientation == face_rotation ? + rotated_refinement_case[equivalent_refinement_case] : + equivalent_refinement_case); + + const unsigned int final_subface_no = + GeometryInfo::child_cell_on_face(RefinementCase( + final_refinement_case), + /*face_no = */ 4, + equivalent_subface_no, + face_orientation, + face_flip, + face_rotation, + equivalent_refinement_case); + + return std::make_pair(final_subface_no, final_refinement_case); + } + else + { + (void)combined_face_orientation; + (void)subface_case; + (void)subface_no; + + DEAL_II_NOT_IMPLEMENTED(); + return {}; + } +} + + + template std::unique_ptr> ReferenceCell::get_default_mapping(const unsigned int degree) const diff --git a/source/grid/reference_cell.inst.in b/source/grid/reference_cell.inst.in index 7693fb3ab9..62697913fe 100644 --- a/source/grid/reference_cell.inst.in +++ b/source/grid/reference_cell.inst.in @@ -32,6 +32,12 @@ for (deal_II_dimension : DIMENSIONS) template const Quadrature & ReferenceCell::get_nodal_type_quadrature() const; + template std::pair> + ReferenceCell::equivalent_refinement_case( + const types::geometric_orientation, + const internal::SubfaceCase, + const unsigned int) const; + template Point ReferenceCell::closest_point( const Point &) const; }