--- /dev/null
+Fixed: Fix GridTools::find_all_active_cells_around_point() to return the correct cells
+if the requested point is located on a cell edge adjacent to neighbor cells at different
+refinement levels.
+<br>
+(Peter Munch, Magdalena Schreter, 2023/02/07)
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)
{
// 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);
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)
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+}
--- /dev/null
+
+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::
projects / dealii.git / commitdiff