From e0ecacfcb9b99da11ff0d0d2c540d6d16f0c945a Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Wed, 11 Jan 2023 12:05:36 -0700
Subject: [PATCH] Implement ReferenceCell::contains_points() for the remaining
 reference cells.

---
 include/deal.II/grid/reference_cell.h | 39 +++++++++++++++++++++++++--
 1 file changed, 37 insertions(+), 2 deletions(-)

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index b149d0dc64..46b34df524 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -2172,11 +2172,46 @@ ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
     }
   else if (*this == ReferenceCells::Wedge)
     {
-      Assert(false, ExcNotImplemented());
+      // The wedge we use is a triangle extruded into the third
+      // dimension by one unit. So we can use the same logic as for
+      // triangles above (i.e., for the simplex above, using dim==2)
+      // and then check the third dimension separately.
+
+      for (unsigned int d = 0; d < 2; ++d)
+        if (p[d] < -tolerance)
+          return false;
+
+      const double sum = p[0] + p[1];
+      if (sum > 1 + tolerance * std::sqrt(2.0))
+        return false;
+
+      if (p[2] < -tolerance)
+        return false;
+      if (p[2] > 1 + tolerance)
+        return false;
+
+      return true;
     }
   else if (*this == ReferenceCells::Pyramid)
     {
-      Assert(false, ExcNotImplemented());
+      // A pyramid only lives in the upper half-space:
+      if (p[2] < -tolerance)
+        return false;
+
+      // It also only lives in the space below z=1:
+      if (p[2] > 1 + tolerance)
+        return false;
+
+      // Within what's left of the space, a pyramid is a cone that tapers
+      // towards the top. First compute the distance of the point to the
+      // axis in the max norm (this is the right norm because the vertices
+      // of the pyramid are at points +/-1, +/-1):
+      const double distance_from_axis =
+        std::max(std::fabs(p[0]), std::fabs(p[1]));
+
+      // We are inside the pyramid if the distance from the axis is less than
+      // (1-z)
+      return (distance_from_axis < 1 + tolerance - p[2]);
     }
 
   Assert(false, ExcNotImplemented());
-- 
2.39.5