]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move some more code from QProjector to ReferenceCell. 18551/head
authorDavid Wells <drwells@email.unc.edu>
Sat, 7 Jun 2025 19:11:21 +0000 (15:11 -0400)
committerDavid Wells <drwells@email.unc.edu>
Thu, 12 Jun 2025 13:18:21 +0000 (09:18 -0400)
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).

include/deal.II/grid/reference_cell.h
source/base/qprojector.cc
source/grid/reference_cell.cc
source/grid/reference_cell.inst.in

index 6a45af7aa408e30ee98ab9d4c6afdd4dac35687d..9546a4930fd30069e80b990e79d75989206dc7d1 100644 (file)
@@ -410,6 +410,18 @@ public:
   std_cxx20::ranges::iota_view<unsigned int, unsigned int>
   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 <int dim>
+  std::pair<unsigned int, RefinementCase<dim - 1>>
+  equivalent_refinement_case(
+    const types::geometric_orientation combined_face_orientation,
+    const internal::SubfaceCase<dim>   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
index dd2bfdbbe7c1c2138ed06f7a1302aac5d055a5b4..11c17f0465d8d7419c6c4c4976114027c26a406b 100644 (file)
@@ -188,146 +188,6 @@ namespace internal
 
 
 
-      std::pair<unsigned int, RefinementCase<2>>
-      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<dim - 1>
-          equivalent_refine_case[internal::SubfaceCase<dim>::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<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy},
-                                  // case_x
-                                  {RefinementCase<dim - 1>::cut_x,
-                                   RefinementCase<dim - 1>::cut_x,
-                                   RefinementCase<dim - 1>::no_refinement,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_x1y
-                                  {RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_x,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_x2y
-                                  {RefinementCase<dim - 1>::cut_x,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_x1y2y
-                                  {RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy},
-                                  // case_y
-                                  {RefinementCase<dim - 1>::cut_y,
-                                   RefinementCase<dim - 1>::cut_y,
-                                   RefinementCase<dim - 1>::no_refinement,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_y1x
-                                  {RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_y,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_y2x
-                                  {RefinementCase<dim - 1>::cut_y,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::no_refinement},
-                                  // case_y1x2x
-                                  {RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy},
-                                  // case_xy (case_isotropic)
-                                  {RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy,
-                                   RefinementCase<dim - 1>::cut_xy}};
-
-        static const unsigned int
-          equivalent_subface_number[internal::SubfaceCase<dim>::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<dim - 1> ref_case_permutation[4] = {
-          RefinementCase<dim - 1>::no_refinement,
-          RefinementCase<dim - 1>::cut_y,
-          RefinementCase<dim - 1>::cut_x,
-          RefinementCase<dim - 1>::cut_xy};
-
-        // set a corresponding (equivalent)
-        // RefineCase and subface number
-        const RefinementCase<dim - 1> 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<dim>::no_refinement,
-               ExcInternalError());
-        Assert(equ_subface_no != e, ExcInternalError());
-        // now, finally respect non-standard faces
-        const RefinementCase<dim - 1> final_ref_case =
-          (face_orientation == face_rotation ?
-             ref_case_permutation[equ_ref_case] :
-             equ_ref_case);
-
-        const unsigned int final_subface_no =
-          GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
-                                                  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<unsigned int, RefinementCase<2>>
-    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<dim>::faces_per_cell * total_subfaces_per_face *
             combined_orientation) *
          n_quadrature_points;
index b1bf43b63edc4509aa20899cafd2db582a01cec7..264355627e19637676faef9fc167dce5ba6e88ff 100644 (file)
@@ -116,6 +116,146 @@ ReferenceCell::to_string() const
 
 
 
+template <int dim>
+std::pair<unsigned int, RefinementCase<dim - 1>>
+ReferenceCell::equivalent_refinement_case(
+  const types::geometric_orientation combined_face_orientation,
+  const internal::SubfaceCase<dim>   subface_case,
+  const unsigned int                 subface_no) const
+{
+  if constexpr (dim == 3)
+    {
+      static const RefinementCase<dim - 1>
+        equivalent_refine_case[internal::SubfaceCase<dim>::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<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy},
+                                // case_x
+                                {RefinementCase<dim - 1>::cut_x,
+                                 RefinementCase<dim - 1>::cut_x,
+                                 RefinementCase<dim - 1>::no_refinement,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_x1y
+                                {RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_x,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_x2y
+                                {RefinementCase<dim - 1>::cut_x,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_x1y2y
+                                {RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy},
+                                // case_y
+                                {RefinementCase<dim - 1>::cut_y,
+                                 RefinementCase<dim - 1>::cut_y,
+                                 RefinementCase<dim - 1>::no_refinement,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_y1x
+                                {RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_y,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_y2x
+                                {RefinementCase<dim - 1>::cut_y,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::no_refinement},
+                                // case_y1x2x
+                                {RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy},
+                                // case_xy (case_isotropic)
+                                {RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy,
+                                 RefinementCase<dim - 1>::cut_xy}};
+
+      constexpr unsigned int X = numbers::invalid_unsigned_int;
+      static const unsigned int
+        equivalent_subface_number[internal::SubfaceCase<dim>::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<dim - 1> rotated_refinement_case[4] = {
+        RefinementCase<dim - 1>::no_refinement,
+        RefinementCase<dim - 1>::cut_y,
+        RefinementCase<dim - 1>::cut_x,
+        RefinementCase<dim - 1>::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<dim>::no_refinement,
+             ExcInternalError());
+      Assert(equivalent_subface_no != X, ExcInternalError());
+      // now, finally respect non-standard faces
+      const RefinementCase<dim - 1> final_refinement_case =
+        (face_orientation == face_rotation ?
+           rotated_refinement_case[equivalent_refinement_case] :
+           equivalent_refinement_case);
+
+      const unsigned int final_subface_no =
+        GeometryInfo<dim>::child_cell_on_face(RefinementCase<dim>(
+                                                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 <int dim, int spacedim>
 std::unique_ptr<Mapping<dim, spacedim>>
 ReferenceCell::get_default_mapping(const unsigned int degree) const
index 7693fb3ab994801f897978847013d7b51c277b9d..62697913fe0f2ced4a40da7d5c3aea69bdab4419 100644 (file)
@@ -32,6 +32,12 @@ for (deal_II_dimension : DIMENSIONS)
     template const Quadrature<deal_II_dimension> &
     ReferenceCell::get_nodal_type_quadrature() const;
 
+    template std::pair<unsigned int, RefinementCase<deal_II_dimension - 1>>
+    ReferenceCell::equivalent_refinement_case(
+      const types::geometric_orientation,
+      const internal::SubfaceCase<deal_II_dimension>,
+      const unsigned int) const;
+
     template Point<deal_II_dimension> ReferenceCell::closest_point(
       const Point<deal_II_dimension> &) const;
   }

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.