Fix find_all_active_cells_around_point()

Co-authored-by: Magdalena Schreter <schreter.magdalena@gmail.com>

source/grid/grid_tools.cc

@@ -6179,6 +6179,9 @@ namespace GridTools
                     tolerance,
                     enforce_unique_mapping);
 
+                if (cell_hint.state() != IteratorState::valid)
+                  cell_hint = cache.get_triangulation().begin_active();
+
                 for (const auto &cell_and_reference_position :
                      cells_and_reference_positions)
                   {

source/grid/grid_tools_dof_handlers.cc

@@ -614,9 +614,10 @@ namespace GridTools
     // need to find all neighbors
     const Point<dim> unit_point = cells_and_points.front().second;
     const auto       my_cell    = cells_and_points.front().first;
-    Tensor<1, dim>   distance_to_center;
-    unsigned int     n_dirs_at_threshold     = 0;
-    unsigned int     last_point_at_threshold = numbers::invalid_unsigned_int;
+
+    Tensor<1, dim> distance_to_center;
+    unsigned int   n_dirs_at_threshold     = 0;
+    unsigned int   last_point_at_threshold = numbers::invalid_unsigned_int;
     for (unsigned int d = 0; d < dim; ++d)
       {
         distance_to_center[d] = std::abs(unit_point[d] - 0.5);
@@ -636,7 +637,18 @@ namespace GridTools
           2 * last_point_at_threshold +
           (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0);
         if (!my_cell->at_boundary(neighbor_index))
-          cells_to_add.push_back(my_cell->neighbor(neighbor_index));
+          {
+            const auto neighbor_cell = my_cell->neighbor(neighbor_index);
+
+            if (neighbor_cell->is_active())
+              cells_to_add.push_back(neighbor_cell);
+            else
+              for (const auto &child_cell : neighbor_cell->child_iterators())
+                {
+                  if (child_cell->is_active())
+                    cells_to_add.push_back(child_cell);
+                }
+          }
       }
     // corner point -> use all neighbors
     else if (n_dirs_at_threshold == dim)
@@ -648,8 +660,10 @@ namespace GridTools
           cells = find_cells_adjacent_to_vertex(
             mesh, my_cell->vertex_index(local_vertex_index));
         for (const auto &cell : cells)
-          if (cell != my_cell)
-            cells_to_add.push_back(cell);
+          {
+            if (cell != my_cell)
+              cells_to_add.push_back(cell);
+          }
       }
     // point on line in 3D: We cannot simply take the intersection between
     // the two vertices of cells because of hanging nodes. So instead we

tests/grid/find_all_active_cells_around_point_04.cc

@@ -0,0 +1,98 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2001 - 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.
+//
+// ---------------------------------------------------------------------
+
+// find_all_active_cells_around_point for neighbor cells at different levels.
+
+#include <deal.II/dofs/dof_handler.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_out.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/tria.h>
+
+#include "../tests.h"
+
+
+template <int dim, int spacedim>
+void
+print_result(const Mapping<dim, spacedim> &      mapping,
+             const Triangulation<dim, spacedim> &tria,
+             const Point<dim>                    p)
+{
+  deallog << "Testing " << dim << "D with point " << p << " tolerance default "
+          << std::endl;
+  auto c_p = GridTools::find_all_active_cells_around_point(mapping, tria, p);
+  for (auto i : c_p)
+    deallog << "Cell: " << i.first->id() << " unit point " << i.second
+            << std::endl;
+  deallog << std::endl;
+}
+
+template <int dim, int spacedim>
+void
+test()
+{
+  Triangulation<dim, spacedim> tria;
+  GridGenerator::subdivided_hyper_cube(
+    tria, 2, 0 /*left*/, 1 /*right*/, true /*colorize*/);
+
+  // In every cycle cells at the boundary face x=0 are refined.
+  for (unsigned int i = 1; i < 3; ++i)
+    {
+      for (const auto &cell : tria.active_cell_iterators())
+        if (cell->at_boundary(0))
+          cell->set_refine_flag();
+
+      tria.execute_coarsening_and_refinement();
+    }
+
+  MappingQ<dim> mapping(2);
+
+  // Point at the edge between cells of different levels
+  Point<dim> point_at_edge;
+
+  point_at_edge[dim - 1] = 0.25;
+  print_result(mapping, tria, point_at_edge);
+
+  if (dim > 1)
+    {
+      for (unsigned d = 0; d < dim; ++d)
+        point_at_edge[d] = 0.25;
+
+      print_result(mapping, tria, point_at_edge);
+    }
+
+  // debug
+  if (false)
+    {
+      std::ofstream output_file("grid_" + std::to_string(dim) + ".vtu");
+      GridOut().write_vtu(tria, output_file);
+    }
+};
+
+
+int
+main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+  test<1, 1>();
+  test<2, 2>();
+  test<3, 3>();
+  return 0;
+}

tests/grid/find_all_active_cells_around_point_04.output

@@ -0,0 +1,29 @@
+
+DEAL::Testing 1D with point 0.2500000000 tolerance default 
+DEAL::Cell: 0_1:1 unit point 0.000000000
+DEAL::Cell: 0_2:01 unit point 1.000000000
+DEAL::
+DEAL::Testing 2D with point 0.000000000 0.2500000000 tolerance default 
+DEAL::Cell: 0_2:02 unit point 0.000000000 1.000000000
+DEAL::Cell: 0_2:20 unit point 0.000000000 0.000000000
+DEAL::
+DEAL::Testing 2D with point 0.2500000000 0.2500000000 tolerance default 
+DEAL::Cell: 0_1:1 unit point 0.000000000 1.000000000
+DEAL::Cell: 0_1:3 unit point 0.000000000 0.000000000
+DEAL::Cell: 0_2:03 unit point 1.000000000 1.000000000
+DEAL::Cell: 0_2:21 unit point 1.000000000 0.000000000
+DEAL::
+DEAL::Testing 3D with point 0.000000000 0.000000000 0.2500000000 tolerance default 
+DEAL::Cell: 0_2:04 unit point 0.000000000 0.000000000 1.000000000
+DEAL::Cell: 0_2:40 unit point 0.000000000 0.000000000 0.000000000
+DEAL::
+DEAL::Testing 3D with point 0.2500000000 0.2500000000 0.2500000000 tolerance default 
+DEAL::Cell: 0_1:1 unit point 0.000000000 1.000000000 1.000000000
+DEAL::Cell: 0_1:3 unit point 0.000000000 0.000000000 1.000000000
+DEAL::Cell: 0_1:5 unit point 0.000000000 1.000000000 0.000000000
+DEAL::Cell: 0_1:7 unit point 0.000000000 0.000000000 0.000000000
+DEAL::Cell: 0_2:07 unit point 1.000000000 1.000000000 1.000000000
+DEAL::Cell: 0_2:25 unit point 1.000000000 0.000000000 1.000000000
+DEAL::Cell: 0_2:43 unit point 1.000000000 1.000000000 0.000000000
+DEAL::Cell: 0_2:61 unit point 1.000000000 0.000000000 0.000000000
+DEAL::