Fix find_all_active_cells_around_point()

dealii · Jan 25, 2023 · 2cf85c6 · 2cf85c6
1 parent 059bd97
commit 2cf85c6
Show file tree

Hide file tree

Showing 2 changed files with 33 additions and 6 deletions.
diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc
@@ -6036,6 +6036,8 @@ namespace GridTools
       cell_hint = first_cell.first;
       if (cell_hint.state() == IteratorState::valid)
         {
+          AssertThrow(cell_hint->is_active(), ExcInternalError());
+
           const auto active_cells_around_point =
             GridTools::find_all_active_cells_around_point(
               cache.get_mapping(),
@@ -6179,6 +6181,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)
                   {

diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc
@@ -614,9 +614,12 @@ 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;
+
+    AssertThrow(my_cell->is_active(), ExcNotImplemented());
+
+    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 +639,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 +662,12 @@ 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);
+          {
+            AssertThrow(cell->is_active(), ExcNotImplemented());
+
+            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
@@ -688,6 +706,8 @@ namespace GridTools
               my_cell->vertex_index(first_vertex + (d << free_direction)));
             for (const auto &cell : tentative_cells)
               {
+                AssertThrow(cell->is_active(), ExcNotImplemented());
+
                 bool cell_not_yet_present = true;
                 for (const auto &other_cell : cells_to_add)
                   if (cell == other_cell)
@@ -705,6 +725,8 @@ namespace GridTools
       GeometryInfo<dim>::distance_to_unit_cell(unit_point);
     for (const auto &cell : cells_to_add)
       {
+        AssertThrow(cell->is_active(), ExcNotImplemented());
+
         if (cell != my_cell)
           try
             {