From ad1b7e826675f8f2abdc8704bcde973eebbd1164 Mon Sep 17 00:00:00 2001
From: Martin Kronbichler <martin.kronbichler@uni-a.de>
Date: Fri, 14 Apr 2023 14:09:35 +0200
Subject: [PATCH] Test case

---
 .../compute_point_locations_try_all_02.cc     | 96 +++++++++++++++++++
 .../compute_point_locations_try_all_02.output | 24 +++++
 2 files changed, 120 insertions(+)
 create mode 100644 tests/grid/compute_point_locations_try_all_02.cc
 create mode 100644 tests/grid/compute_point_locations_try_all_02.output

diff --git a/tests/grid/compute_point_locations_try_all_02.cc b/tests/grid/compute_point_locations_try_all_02.cc
new file mode 100644
index 0000000000..3519200747
--- /dev/null
+++ b/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}});
+}
diff --git a/tests/grid/compute_point_locations_try_all_02.output b/tests/grid/compute_point_locations_try_all_02.output
new file mode 100644
index 0000000000..2b11a01782
--- /dev/null
+++ b/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::
-- 
2.39.5