Merge pull request #13496 from bangerth/point_is_inside

Implement ReferenceCell::point_is_inside.

doc/news/changes/minor/20220305Bangerth

@@ -0,0 +1,6 @@
+New: There is now a function ReferenceCell::contains_point() that
+provides the generalized equivalent of a function in the GeometryInfo
+class that will eventually be deprecated.
+type.
+<br>
+(Wolfgang Bangerth, 2022/03/05)

include/deal.II/grid/reference_cell.h

@@ -457,6 +457,29 @@ class ReferenceCell
    * @{
    */
 
+  /**
+   * Return true if the given point is inside the reference cell of the present
+   * space dimension up to some tolerance. This function accepts an additional
+   * parameter (which defaults to zero) which specifies by how much the point
+   * position may actually be outside the true reference cell. This is useful
+   * because in practice we may often not be able to compute the coordinates of
+   * a point in reference coordinates exactly, but only up to numerical
+   * roundoff. For example, strictly speaking one would expect that for points
+   * on the boundary of the reference cell, the function would return `true` if
+   * the tolerance was zero. But in practice, this may or may not actually be
+   * true; for example, the point $(1/3, 2/3)$ is on the boundary of the
+   * reference triangle because $1/3+2/3 \le 1$, but since neither of its
+   * coordinates are exactly representable in floating point arithmetic, the
+   * floating point representations of $1/3$ and $2/3$ may or may not add up to
+   * anything that is less than or equal to one.
+   *
+   * The tolerance parameter may be less than zero, indicating that the point
+   * should be safely inside the cell.
+   */
+  template <int dim>
+  bool
+  contains_point(const Point<dim> &p, const double tolerance = 0) const;
+
   /*
    * Return $i$-th unit tangential vector of a face of the reference cell.
    * The vectors are arranged such that the
@@ -1898,6 +1921,69 @@ ReferenceCell::d_linear_shape_function_gradient(const Point<dim> & xi,
 }
 
 
+
+template <int dim>
+inline bool
+ReferenceCell::contains_point(const Point<dim> &p, const double tolerance) const
+{
+  AssertDimension(dim, get_dimension());
+
+  if (*this == ReferenceCells::Vertex)
+    {
+      // Vertices are special cases in that they do not actually
+      // have coordinates. Error out if this function is called
+      // with a vertex:
+      Assert(false,
+             ExcMessage("Vertices are zero-dimensional objects and "
+                        "as a consequence have no coordinates. You "
+                        "cannot meaningfully ask whether a point is "
+                        "inside a vertex (within a certain tolerance) "
+                        "without coordinate values."));
+      return false;
+    }
+  else if (*this == ReferenceCells::get_hypercube<dim>())
+    {
+      for (unsigned int d = 0; d < dim; ++d)
+        if ((p[d] < -tolerance) || (p[d] > 1 + tolerance))
+          return false;
+      return true;
+    }
+  else if (*this == ReferenceCells::get_simplex<dim>())
+    {
+      // First make sure that we are in the first quadrant or octant
+      for (unsigned int d = 0; d < dim; ++d)
+        if (p[d] < -tolerance)
+          return false;
+
+      // Now we also need to make sure that we are below the diagonal line
+      // or plane that delineates the simplex. This diagonal is given by
+      // sum(p[d])<=1, and a diagonal a distance eps away is given by
+      // sum(p[d])<=1+eps*sqrt(d). (For example, the point at (1,1) is a
+      // distance of 1/sqrt(2) away from the diagonal. That is, its
+      // sum satisfies
+      //   sum(p[d]) = 2 <= 1 + (1/sqrt(2)) * sqrt(2)
+      // in other words, it satisfies the predicate with eps=1/sqrt(2).)
+      double sum = 0;
+      for (unsigned int d = 0; d < dim; ++d)
+        sum += p[d];
+      return (sum <= 1 + tolerance * std::sqrt(1. * dim));
+    }
+  else if (*this == ReferenceCells::Wedge)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+  else if (*this == ReferenceCells::Pyramid)
+    {
+      Assert(false, ExcNotImplemented());
+    }
+
+  Assert(false, ExcNotImplemented());
+
+  return false;
+}
+
+
+
 template <int dim>
 inline Tensor<1, dim>
 ReferenceCell::unit_tangential_vectors(const unsigned int face_no,

tests/simplex/reference_cell_point_is_inside_01.cc

@@ -0,0 +1,72 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2003 - 2020 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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check ReferenceCell::contains_point()
+
+#include <deal.II/grid/reference_cell.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test(const ReferenceCell &reference_cell)
+{
+  deallog << "Reference cell is: " << reference_cell.to_string() << std::endl;
+
+  deallog << std::boolalpha;
+
+  // Generate equispaced points and check whether they are inside the box
+  for (double x = -0.1; x <= 1.1; x += 0.1)
+    for (double y = -0.1; y <= (dim > 1 ? 1.1 : -0.1); y += 0.1)
+      for (double z = -0.1; z <= (dim > 2 ? 1.1 : -0.1); z += 0.1)
+        {
+          Point<dim> p =
+            (dim == 1 ? Point<dim>(x) :
+                        (dim == 2 ? Point<dim>(x, y) : Point<dim>(x, y, z)));
+
+          deallog << p << ' ' << reference_cell.contains_point(p) << std::endl;
+        }
+
+  // Make sure that all vertices are inside:
+  for (unsigned int v = 0; v < reference_cell.n_vertices(); ++v)
+    Assert(reference_cell.contains_point(reference_cell.vertex<dim>(v)),
+           ExcInternalError());
+
+  // Make sure that all vertices are outside with a negative tolerance:
+  for (unsigned int v = 0; v < reference_cell.n_vertices(); ++v)
+    Assert(reference_cell.contains_point(reference_cell.vertex<dim>(v),
+                                         -1e-12) == false,
+           ExcInternalError());
+}
+
+
+int
+main()
+{
+  initlog();
+
+  test<1>(ReferenceCells::Line);
+  test<2>(ReferenceCells::Quadrilateral);
+  test<2>(ReferenceCells::Triangle);
+  test<3>(ReferenceCells::Hexahedron);
+  test<3>(ReferenceCells::Tetrahedron);
+
+  // TODO: wedge and pyramid are not currently implemented
+
+  return 0;
+}