From eb3400d2007896cd268b4b161c540f7d4a48e94a Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 30 Jul 2023 22:30:14 +0200 Subject: [PATCH] GridTools::find_active_cell_around_point() and find_all_active_cells_around_point() for simplices --- source/grid/grid_tools.cc | 8 +- source/grid/grid_tools_dof_handlers.cc | 194 +++++++----- .../find_all_active_cells_around_point_01.cc | 72 +++++ ...nd_all_active_cells_around_point_01.output | 283 ++++++++++++++++++ 4 files changed, 474 insertions(+), 83 deletions(-) create mode 100644 tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc create mode 100644 tests/remote_point_evaluation/find_all_active_cells_around_point_01.output diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 7045acce08..e840d8955d 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -2998,8 +2998,8 @@ namespace GridTools { const Point p_unit = mapping.transform_real_to_unit_cell(*cell, p); - if (GeometryInfo::is_inside_unit_cell(p_unit, - tolerance)) + if ((*cell)->reference_cell().contains_point(p_unit, + tolerance)) { cell_and_position.first = *cell; cell_and_position.second = p_unit; @@ -3011,8 +3011,8 @@ namespace GridTools // outside it is and whether we want to use this cell as a // backup if we can't find a cell within which the point // lies. - const double dist = - GeometryInfo::distance_to_unit_cell(p_unit); + const double dist = p_unit.distance( + (*cell)->reference_cell().closest_point(p_unit)); if (dist < best_distance) { best_distance = dist; diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 9a93f8668f..e4bde46a76 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -629,97 +629,137 @@ namespace GridTools // insert the fist cell and point into the vector cells_and_points.push_back(first_cell); - // check if the given point is on the surface of the unit cell. If yes, - // 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; - for (unsigned int d = 0; d < dim; ++d) + std::vector::active_cell_iterator> + cells_to_add; + + if (my_cell->reference_cell().is_hyper_cube()) { - distance_to_center[d] = std::abs(unit_point[d] - 0.5); - if (distance_to_center[d] > 0.5 - tolerance) + // check if the given point is on the surface of the unit cell. If yes, + // need to find all neighbors + + 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) { - ++n_dirs_at_threshold; - last_point_at_threshold = d; + distance_to_center[d] = std::abs(unit_point[d] - 0.5); + if (distance_to_center[d] > 0.5 - tolerance) + { + ++n_dirs_at_threshold; + last_point_at_threshold = d; + } } - } - std::vector::active_cell_iterator> - cells_to_add; - // point is within face -> only need neighbor - if (n_dirs_at_threshold == 1) - { - unsigned int neighbor_index = - 2 * last_point_at_threshold + - (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0); - if (!my_cell->at_boundary(neighbor_index)) + // point is within face -> only need neighbor + if (n_dirs_at_threshold == 1) { - const auto neighbor_cell = my_cell->neighbor(neighbor_index); + unsigned int neighbor_index = + 2 * last_point_at_threshold + + (unit_point[last_point_at_threshold] > 0.5 ? 1 : 0); + if (!my_cell->at_boundary(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); - } + 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) - { - unsigned int local_vertex_index = 0; - for (unsigned int d = 0; d < dim; ++d) - local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d; + // corner point -> use all neighbors + else if (n_dirs_at_threshold == dim) + { + unsigned int local_vertex_index = 0; + for (unsigned int d = 0; d < dim; ++d) + local_vertex_index += (unit_point[d] > 0.5 ? 1 : 0) << d; - const auto fu = [&](const auto &tentative_cells) { - for (const auto &cell : tentative_cells) - if (cell != my_cell) - cells_to_add.push_back(cell); - }; + const auto fu = [&](const auto &tentative_cells) { + for (const auto &cell : tentative_cells) + if (cell != my_cell) + cells_to_add.push_back(cell); + }; - const auto vertex_index = my_cell->vertex_index(local_vertex_index); + const auto vertex_index = my_cell->vertex_index(local_vertex_index); - if (vertex_to_cells != nullptr) - fu((*vertex_to_cells)[vertex_index]); - else - fu(find_cells_adjacent_to_vertex(mesh, vertex_index)); - } - // point on line in 3d: We cannot simply take the intersection between - // the two vertices of cells because of hanging nodes. So instead we - // list the vertices around both points and then select the - // appropriate cells according to the result of read_to_unit_cell - // below. - else if (n_dirs_at_threshold == 2) - { - std::pair vertex_indices[3]; - unsigned int count_vertex_indices = 0; - unsigned int free_direction = numbers::invalid_unsigned_int; - for (unsigned int d = 0; d < dim; ++d) + if (vertex_to_cells != nullptr) + fu((*vertex_to_cells)[vertex_index]); + else + fu(find_cells_adjacent_to_vertex(mesh, vertex_index)); + } + // point on line in 3d: We cannot simply take the intersection between + // the two vertices of cells because of hanging nodes. So instead we + // list the vertices around both points and then select the + // appropriate cells according to the result of read_to_unit_cell + // below. + else if (n_dirs_at_threshold == 2) { - if (distance_to_center[d] > 0.5 - tolerance) + std::pair vertex_indices[3]; + unsigned int count_vertex_indices = 0; + unsigned int free_direction = numbers::invalid_unsigned_int; + for (unsigned int d = 0; d < dim; ++d) { - vertex_indices[count_vertex_indices].first = d; - vertex_indices[count_vertex_indices].second = - unit_point[d] > 0.5 ? 1 : 0; - count_vertex_indices++; + if (distance_to_center[d] > 0.5 - tolerance) + { + vertex_indices[count_vertex_indices].first = d; + vertex_indices[count_vertex_indices].second = + unit_point[d] > 0.5 ? 1 : 0; + count_vertex_indices++; + } + else + free_direction = d; } - else - free_direction = d; - } - AssertDimension(count_vertex_indices, 2); - Assert(free_direction != numbers::invalid_unsigned_int, - ExcInternalError()); + AssertDimension(count_vertex_indices, 2); + Assert(free_direction != numbers::invalid_unsigned_int, + ExcInternalError()); + + const unsigned int first_vertex = + (vertex_indices[0].second << vertex_indices[0].first) + + (vertex_indices[1].second << vertex_indices[1].first); + for (unsigned int d = 0; d < 2; ++d) + { + const auto fu = [&](const auto &tentative_cells) { + for (const auto &cell : tentative_cells) + { + bool cell_not_yet_present = true; + for (const auto &other_cell : cells_to_add) + if (cell == other_cell) + { + cell_not_yet_present = false; + break; + } + if (cell_not_yet_present) + cells_to_add.push_back(cell); + } + }; + + const auto vertex_index = + my_cell->vertex_index(first_vertex + (d << free_direction)); + + if (vertex_to_cells != nullptr) + fu((*vertex_to_cells)[vertex_index]); + else + fu(find_cells_adjacent_to_vertex(mesh, vertex_index)); + } + } + } + else + { + // Note: The non-hypercube path takes a very naive approach and + // checks all possible neighbors. This can be made faster by 1) + // checking if the point is in the inner cell and 2) identifying + // the right lines/vertices so that the number of potential + // neighbors is reduced. - const unsigned int first_vertex = - (vertex_indices[0].second << vertex_indices[0].first) + - (vertex_indices[1].second << vertex_indices[1].first); - for (unsigned int d = 0; d < 2; ++d) + for (const auto v : my_cell->vertex_indices()) { const auto fu = [&](const auto &tentative_cells) { for (const auto &cell : tentative_cells) @@ -736,8 +776,7 @@ namespace GridTools } }; - const auto vertex_index = - my_cell->vertex_index(first_vertex + (d << free_direction)); + const auto vertex_index = my_cell->vertex_index(v); if (vertex_to_cells != nullptr) fu((*vertex_to_cells)[vertex_index]); @@ -746,8 +785,6 @@ namespace GridTools } } - const double original_distance_to_unit_cell = - my_cell->reference_cell().closest_point(unit_point).distance(unit_point); for (const auto &cell : cells_to_add) { if (cell != my_cell) @@ -755,8 +792,7 @@ namespace GridTools { const Point p_unit = mapping.transform_real_to_unit_cell(cell, p); - if (cell->reference_cell().closest_point(p_unit).distance( - p_unit) < original_distance_to_unit_cell + tolerance) + if (cell->reference_cell().contains_point(p_unit, tolerance)) cells_and_points.emplace_back(cell, p_unit); } catch (typename Mapping::ExcTransformationFailed &) diff --git a/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc new file mode 100644 index 0000000000..1d6046adab --- /dev/null +++ b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.cc @@ -0,0 +1,72 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 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::find_active_cell_around_point() and +// GridTools::find_all_active_cells_around_point for simplices. These +// functions are used in find_all_locally_owned_active_cells_around_point(), +// which is used by distributed_compute_point_locations(). + +#include +#include +#include +#include + +#include "../tests.h" + +template +void +test() +{ + const auto reference_cell = ReferenceCells::get_simplex(); + const auto &mapping = + reference_cell.template get_default_linear_mapping(); + + Triangulation tria; + if (false) + GridGenerator::reference_cell(tria, reference_cell); + else + GridGenerator::subdivided_hyper_cube_with_simplices(tria, 2); + + GridTools::Cache cache(tria, mapping); + + const unsigned n_subdivisions = 8; + + for (unsigned int i = 0; i <= n_subdivisions; ++i) + for (unsigned int j = 0; j <= n_subdivisions; ++j) + { + Point p(1.0 * i / n_subdivisions, 1.0 * j / n_subdivisions); + + const auto first_point = + GridTools::find_active_cell_around_point(mapping, tria, p, {}, 1e-6); + + const auto result = GridTools::find_all_active_cells_around_point( + mapping, tria, p, 1e-6, {first_point.first, first_point.second}); + + deallog << result.size() << std::endl; // should be 2 + for (const auto &cell : result) + deallog << cell.first->id() << " " << cell.second << std::endl; + deallog << std::endl; + } +} + +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + + deallog.precision(8); + + test<2>(); +} diff --git a/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output new file mode 100644 index 0000000000..70c3599ed0 --- /dev/null +++ b/tests/remote_point_evaluation/find_all_active_cells_around_point_01.output @@ -0,0 +1,283 @@ + +DEAL:0::1 +DEAL:0::0_0: 0.0000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::3 +DEAL:0::0_0: 0.0000000 1.0000000 +DEAL:0::1_0: 1.0000000 0.0000000 +DEAL:0::4_0: 0.0000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::4_0: 0.0000000 1.0000000 +DEAL:0::5_0: 1.0000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::0_0: 0.25000000 0.75000000 +DEAL:0::1_0: 0.75000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.75000000 0.0000000 +DEAL:0::4_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::4_0: 0.25000000 0.75000000 +DEAL:0::5_0: 0.75000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::0_0: 0.50000000 0.50000000 +DEAL:0::1_0: 0.50000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::1_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.50000000 0.0000000 +DEAL:0::4_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::4_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::4_0: 0.50000000 0.50000000 +DEAL:0::5_0: 0.50000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::0_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::0_0: 0.75000000 0.25000000 +DEAL:0::1_0: 0.25000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::1_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::1_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.25000000 0.0000000 +DEAL:0::4_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::4_0: 0.75000000 0.25000000 +DEAL:0::5_0: 0.25000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::5_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::3 +DEAL:0::0_0: 1.0000000 0.0000000 +DEAL:0::1_0: 0.0000000 1.0000000 +DEAL:0::2_0: 0.0000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.0000000 0.75000000 +DEAL:0::2_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.0000000 0.50000000 +DEAL:0::2_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::1_0: 0.0000000 0.25000000 +DEAL:0::2_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::6 +DEAL:0::1_0: 0.0000000 0.0000000 +DEAL:0::2_0: 0.0000000 1.0000000 +DEAL:0::3_0: 1.0000000 0.0000000 +DEAL:0::4_0: 1.0000000 0.0000000 +DEAL:0::5_0: 0.0000000 1.0000000 +DEAL:0::6_0: 0.0000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::5_0: 0.0000000 0.75000000 +DEAL:0::6_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::5_0: 0.0000000 0.50000000 +DEAL:0::6_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::5_0: 0.0000000 0.25000000 +DEAL:0::6_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::3 +DEAL:0::5_0: 0.0000000 0.0000000 +DEAL:0::6_0: 0.0000000 1.0000000 +DEAL:0::7_0: 1.0000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::2_0: 0.25000000 0.75000000 +DEAL:0::3_0: 0.75000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::3_0: 0.75000000 0.0000000 +DEAL:0::6_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::6_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::6_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::6_0: 0.25000000 0.75000000 +DEAL:0::7_0: 0.75000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::2_0: 0.50000000 0.50000000 +DEAL:0::3_0: 0.50000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::3_0: 0.50000000 0.0000000 +DEAL:0::6_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::6_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::6_0: 0.50000000 0.50000000 +DEAL:0::7_0: 0.50000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.50000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.50000000 0.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::2_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::2_0: 0.75000000 0.25000000 +DEAL:0::3_0: 0.25000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::3_0: 0.25000000 0.0000000 +DEAL:0::6_0: 0.75000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::6_0: 0.75000000 0.25000000 +DEAL:0::7_0: 0.25000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.25000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.25000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.25000000 0.0000000 +DEAL:0:: +DEAL:0::2 +DEAL:0::2_0: 1.0000000 0.0000000 +DEAL:0::3_0: 0.0000000 1.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::3_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::3 +DEAL:0::3_0: 0.0000000 0.0000000 +DEAL:0::6_0: 1.0000000 0.0000000 +DEAL:0::7_0: 0.0000000 1.0000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.0000000 0.75000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.0000000 0.50000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.0000000 0.25000000 +DEAL:0:: +DEAL:0::1 +DEAL:0::7_0: 0.0000000 0.0000000 +DEAL:0:: -- 2.39.5