From ca161a5b65548f812257f642b5522779d3f6307f Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 25 Jan 2023 23:17:53 +0100 Subject: [PATCH] Fix find_all_active_cells_around_point() Co-authored-by: Magdalena Schreter --- doc/news/changes/minor/20230207MunchSchreter | 5 + source/grid/grid_tools.cc | 3 + source/grid/grid_tools_dof_handlers.cc | 26 +++-- .../find_all_active_cells_around_point_04.cc | 98 +++++++++++++++++++ ...nd_all_active_cells_around_point_04.output | 29 ++++++ 5 files changed, 155 insertions(+), 6 deletions(-) create mode 100644 doc/news/changes/minor/20230207MunchSchreter create mode 100644 tests/grid/find_all_active_cells_around_point_04.cc create mode 100644 tests/grid/find_all_active_cells_around_point_04.output diff --git a/doc/news/changes/minor/20230207MunchSchreter b/doc/news/changes/minor/20230207MunchSchreter new file mode 100644 index 0000000000..ad67a27000 --- /dev/null +++ b/doc/news/changes/minor/20230207MunchSchreter @@ -0,0 +1,5 @@ +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. +
+(Peter Munch, Magdalena Schreter, 2023/02/07) diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index cad7fc8f3e..f5d400f7b3 100644 --- a/source/grid/grid_tools.cc +++ b/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) { diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 18a0a51b49..0378fc4646 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -614,9 +614,10 @@ namespace GridTools // need to find all neighbors const Point 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 diff --git a/tests/grid/find_all_active_cells_around_point_04.cc b/tests/grid/find_all_active_cells_around_point_04.cc new file mode 100644 index 0000000000..4654d92e41 --- /dev/null +++ b/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 + +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + + +template +void +print_result(const Mapping & mapping, + const Triangulation &tria, + const Point 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 +void +test() +{ + Triangulation 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 mapping(2); + + // Point at the edge between cells of different levels + Point 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; +} diff --git a/tests/grid/find_all_active_cells_around_point_04.output b/tests/grid/find_all_active_cells_around_point_04.output new file mode 100644 index 0000000000..1652b45e37 --- /dev/null +++ b/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:: -- 2.39.5