Merge pull request #14670 from bangerth/contains-point

Implement ReferenceCell::contains_points() for the remaining reference cells.

doc/news/changes/minor/20230111Bangerth

@@ -0,0 +1,5 @@
+New: The function ReferenceCell::contains_points() is now implemented
+for all reference cell kinds (it had been missing for wedges and
+pyramids before).
+<br>
+(Wolfgang Bangerth, 2023/01/11)

include/deal.II/grid/reference_cell.h

@@ -2211,6 +2211,13 @@ ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
 {
   AssertDimension(dim, get_dimension());
 
+  // Introduce abbreviations to silence compiler warnings about invalid
+  // array accesses (that can't happen, of course, but the compiler
+  // doesn't know that).
+  constexpr unsigned int x_coordinate = 0;
+  constexpr unsigned int y_coordinate = (dim >= 2 ? 1 : 0);
+  constexpr unsigned int z_coordinate = (dim >= 3 ? 2 : 0);
+
   if (*this == ReferenceCells::Vertex)
     {
       // Vertices are special cases in that they do not actually
@@ -2253,11 +2260,45 @@ 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.
+
+      if ((p[x_coordinate] < -tolerance) || (p[y_coordinate] < -tolerance))
+        return false;
+
+      const double sum = p[x_coordinate] + p[y_coordinate];
+      if (sum > 1 + tolerance * std::sqrt(2.0))
+        return false;
+
+      if (p[z_coordinate] < -tolerance)
+        return false;
+      if (p[z_coordinate] > 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[z_coordinate] < -tolerance)
+        return false;
+
+      // It also only lives in the space below z=1:
+      if (p[z_coordinate] > 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[x_coordinate]), std::fabs(p[y_coordinate]));
+
+      // We are inside the pyramid if the distance from the axis is less than
+      // (1-z)
+      return (distance_from_axis < 1 + tolerance - p[z_coordinate]);
     }
 
   Assert(false, ExcNotImplemented());

tests/grid/reference_cell_type_04.cc

@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2020 - 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+// Test consistency of ReferenceCell::volume() and
+// ReferenceCell::contains_point() by computing a Monte Carlo integral
+// of the volume based on the characteristic function implied by
+// contains_point().
+
+#include <deal.II/grid/reference_cell.h>
+
+#include <random>
+
+#include "../tests.h"
+
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell)
+{
+  unsigned int       n_samples_inside = 0;
+  const unsigned int n_samples        = 200000;
+
+  for (unsigned int n = 0; n < n_samples; ++n)
+    {
+      // Choose a random point in the box [-1,1]^d that contains all
+      // of our reference cells:
+      Point<dim> p;
+      for (unsigned int d = 0; d < dim; ++d)
+        p[d] = random_value<double>(-1, 1);
+
+      if (reference_cell.contains_point(p))
+        ++n_samples_inside;
+    }
+
+  // The approximate volume of the reference cell is then the fraction
+  // of points that were found to be within the reference cell times
+  // the volume of the box within which we have sampled:
+  const double volume = 1. * n_samples_inside / n_samples * std::pow(2.0, dim);
+
+  deallog << "ReferenceCell: " << reference_cell.to_string() << std::endl;
+  deallog << "  self-reported volume = " << reference_cell.volume()
+          << std::endl;
+  deallog << "  computed approximate volume = " << volume << std::endl;
+
+  Assert(std::fabs(volume - reference_cell.volume()) < 1e-2,
+         ExcInternalError());
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("1D");
+    test<1>(ReferenceCells::Line);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("2D");
+    test<2>(ReferenceCells::Quadrilateral);
+    test<2>(ReferenceCells::Triangle);
+    deallog.pop();
+  }
+
+  {
+    deallog.push("3D");
+    test<3>(ReferenceCells::Tetrahedron);
+    test<3>(ReferenceCells::Pyramid);
+    test<3>(ReferenceCells::Wedge);
+    test<3>(ReferenceCells::Hexahedron);
+    deallog.pop();
+  }
+}

tests/grid/reference_cell_type_04.output

@@ -0,0 +1,22 @@
+
+DEAL:1D::ReferenceCell: Line
+DEAL:1D::  self-reported volume = 1.00000
+DEAL:1D::  computed approximate volume = 1.00076
+DEAL:2D::ReferenceCell: Quad
+DEAL:2D::  self-reported volume = 1.00000
+DEAL:2D::  computed approximate volume = 1.00358
+DEAL:2D::ReferenceCell: Tri
+DEAL:2D::  self-reported volume = 0.500000
+DEAL:2D::  computed approximate volume = 0.497740
+DEAL:3D::ReferenceCell: Tet
+DEAL:3D::  self-reported volume = 0.166667
+DEAL:3D::  computed approximate volume = 0.169240
+DEAL:3D::ReferenceCell: Pyramid
+DEAL:3D::  self-reported volume = 1.33333
+DEAL:3D::  computed approximate volume = 1.34152
+DEAL:3D::ReferenceCell: Wedge
+DEAL:3D::  self-reported volume = 0.500000
+DEAL:3D::  computed approximate volume = 0.499080
+DEAL:3D::ReferenceCell: Hex
+DEAL:3D::  self-reported volume = 1.00000
+DEAL:3D::  computed approximate volume = 1.00820