]> https://gitweb.dealii.org/ - dealii.git/commitdiff
QProjector: add another project_to_subface() overload.
authorDavid Wells <drwells@email.unc.edu>
Sat, 15 Feb 2025 14:07:41 +0000 (09:07 -0500)
committerDavid Wells <drwells@email.unc.edu>
Sat, 15 Feb 2025 20:29:45 +0000 (15:29 -0500)
include/deal.II/base/qprojector.h
source/base/qprojector.cc

index 07643e71ab266e7e1add4f867c3ce6d1150f22b7..7aa720229ebeb8eb6b893b8ca68c762f2283d4bc 100644 (file)
@@ -209,6 +209,23 @@ public:
                               const bool                       face_rotation,
                               const internal::SubfaceCase<dim> ref_case);
 
+  /**
+   * Compute the cell quadrature formula corresponding to using
+   * <tt>quadrature</tt> on subface <tt>subface_no</tt> of face
+   * <tt>face_no</tt> with RefinementCase<dim-1> <tt>ref_case</tt>. 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<dim>
+  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<dim - 1>     &ref_case);
+
   /**
    * Take a collection of face quadrature formulas and generate a cell
    * quadrature formula from it where the quadrature points of the given
index 5a8b0ee74630efe15a7b677942745a90d0c6a7e3..1af849c4b420f094536e32a70c798759d2085c43 100644 (file)
@@ -658,37 +658,122 @@ QProjector<dim>::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 <int dim>
+Quadrature<dim>
+QProjector<dim>::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<dim - 1>     &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<dim - 1>(ref_case));
 
-  const Quadrature<2> mutation = internal::QProjector::mutate_quadrature(
-    quadrature,
-    internal::combined_face_orientation(face_orientation,
-                                        face_rotation,
-                                        face_flip));
+  std::vector<Point<dim>> q_points;
+  std::vector<double>     q_weights = quadrature.get_weights();
+  q_points.reserve(quadrature.size());
 
-  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);
+  if constexpr (dim == 1)
+    {
+      AssertDimension(quadrature.size(), 1);
+      q_points.emplace_back(static_cast<double>(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<dim>(std::move(q_points), std::move(q_weights));
 }
 
 

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.