Add tolerance to box in compute_point_locations_try_all

include/deal.II/base/bounding_box.h

@@ -256,6 +256,22 @@ class BoundingBox
   BoundingBox<spacedim, Number>
   create_extended(const Number amount) const;
 
+  /**
+   * Increase (or decrease) each side of the bounding box by the given
+   * @p relative_amount.
+   *
+   * After calling this method, the lower left corner of the bounding box will
+   * have each coordinate decreased by @p relative_amount * side_length(direction),
+   * and the upper right corner of the bounding box will have each coordinate
+   * increased by @p relative_amount * side_length(direction).
+   *
+   * If you call this method with a negative number, and one of the axes of the
+   * original bounding box is smaller than relative_amount *
+   * side_length(direction) / 2, the method will trigger an assertion.
+   */
+  BoundingBox<spacedim, Number>
+  create_extended_relative(const Number relative_amount) const;
+
   /**
    * Compute the volume (i.e. the dim-dimensional measure) of the BoundingBox.
    */
@@ -574,6 +590,29 @@ BoundingBox<spacedim, Number>::create_extended(const Number amount) const
 }
 
 
+
+template <int spacedim, typename Number>
+inline BoundingBox<spacedim, Number>
+BoundingBox<spacedim, Number>::create_extended_relative(
+  const Number relative_amount) const
+{
+  // create and modify copy
+  auto bb = *this;
+
+  for (unsigned int d = 0; d < spacedim; ++d)
+    {
+      bb.boundary_points.first[d] -= relative_amount * side_length(d);
+      bb.boundary_points.second[d] += relative_amount * side_length(d);
+      Assert(bb.boundary_points.first[d] <= bb.boundary_points.second[d],
+             ExcMessage("Bounding Box can't be shrunk this much: the points' "
+                        "order should remain bottom left, top right."));
+    }
+
+  return bb;
+}
+
+
+
 template <int spacedim, typename Number>
 template <class Archive>
 void

source/grid/grid_tools.cc

@@ -5829,7 +5829,9 @@ namespace GridTools
 
     // Check all points within a given pair of box and cell
     const auto check_all_points_within_box = [&](const auto &leaf) {
-      const auto &box       = leaf.first;
+      const double                relative_tolerance = 1e-12;
+      const BoundingBox<spacedim> box =
+        leaf.first.create_extended_relative(relative_tolerance);
       const auto &cell_hint = leaf.second;
 
       for (const auto &point_and_id :

tests/grid/compute_point_locations_try_all_02.cc

@@ -0,0 +1,96 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2019 - 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.
+//
+// ---------------------------------------------------------------------
+
+// Test GridTools::compute_point_locations_try_all for the case where
+// some points only lie outside a cell by some roundoff, and another case
+// where they are farther away
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+check_found(const Triangulation<dim> &     tria,
+            const GridTools::Cache<dim> &  cache,
+            const std::vector<Point<dim>> &points)
+{
+  auto cell_qpoint_map =
+    GridTools::compute_point_locations_try_all(cache, points);
+  const auto &cells   = std::get<0>(cell_qpoint_map);
+  const auto &qpoints = std::get<1>(cell_qpoint_map);
+  const auto &indices = std::get<2>(cell_qpoint_map);
+  const auto &other_p = std::get<3>(cell_qpoint_map);
+  size_t      n_cells = cells.size();
+
+  deallog << "Points found in " << n_cells << " cells" << std::endl;
+
+  for (unsigned int i = 0; i < n_cells; ++i)
+    {
+      deallog << "On cell " << cells[i]->id() << " found:" << std::endl;
+      unsigned int j = 0;
+      for (const Point<dim> &p : qpoints[i])
+        deallog << "real " << points[indices[i][j++]] << " unit " << p
+                << std::endl;
+    }
+  deallog << "Points not found: ";
+  for (const auto i : other_p)
+    deallog << points[i] << "   ";
+  deallog << std::endl << std::endl;
+}
+
+
+int
+main()
+{
+  initlog();
+
+  deallog << std::setprecision(18);
+
+  const int          dim = 2;
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -0.3, 0.7);
+  triangulation.refine_global(2);
+
+  MappingQ<dim>         mapping(1);
+  GridTools::Cache<dim> cache(triangulation, mapping);
+
+  // Point is outside by 1e-16
+  Point<dim> p0(-0.299999999999999989, -0.247168783648703205),
+    p1(-0.300000000000000044, -0.102831216351296786);
+
+  check_found(triangulation, cache, {{p0, p1}});
+
+  // Shift point by 1e-8
+  p0[0] -= 1e-8;
+
+  check_found(triangulation, cache, {{p0, p1}});
+
+  // Shift point by more
+  p0[0] += 0.5;
+  p1[0] += 1e-15;
+
+  check_found(triangulation, cache, {{p0, p1}});
+
+  // Shift outside
+  p0[0] += 1.;
+
+  check_found(triangulation, cache, {{p0, p1}});
+}

tests/grid/compute_point_locations_try_all_02.output

@@ -0,0 +1,24 @@
+
+DEAL::Points found in 1 cells
+DEAL::On cell 0_2:00 found:
+DEAL::real -0.299999999999999989 -0.247168783648703205 unit 0.00000000000000000 0.211324865405187134
+DEAL::real -0.300000000000000044 -0.102831216351296786 unit -2.22044604925031308e-16 0.788675134594812866
+DEAL::Points not found: 
+DEAL::
+DEAL::Points found in 1 cells
+DEAL::On cell 0_2:00 found:
+DEAL::real -0.300000000000000044 -0.102831216351296786 unit -2.22044604925031308e-16 0.788675134594812866
+DEAL::Points not found: -0.300000009999999984 -0.247168783648703205   
+DEAL::
+DEAL::Points found in 2 cells
+DEAL::On cell 0_2:01 found:
+DEAL::real 0.199999990000000016 -0.247168783648703205 unit 0.999999960000000132 0.211324865405187134
+DEAL::On cell 0_2:00 found:
+DEAL::real -0.299999999999999045 -0.102831216351296786 unit 3.77475828372553224e-15 0.788675134594812866
+DEAL::Points not found: 
+DEAL::
+DEAL::Points found in 1 cells
+DEAL::On cell 0_2:00 found:
+DEAL::real -0.299999999999999045 -0.102831216351296786 unit 3.77475828372553224e-15 0.788675134594812866
+DEAL::Points not found: 1.19999999000000002 -0.247168783648703205   
+DEAL::