From e4913b0f5d41bbfd3c5f39f8e303056dbfd20f70 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Mon, 31 Jul 2023 15:37:38 +0200 Subject: [PATCH] Use Cache::get_vertex_to_cell_map() in GT::find_all_locally_owned_active_cells_around_point() --- include/deal.II/grid/grid_tools.h | 10 +++- source/grid/grid_tools.cc | 3 +- source/grid/grid_tools_dof_handlers.cc | 59 +++++++++++++-------- source/grid/grid_tools_dof_handlers.inst.in | 3 +- 4 files changed, 50 insertions(+), 25 deletions(-) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index c437c55aba..073cf10d68 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1810,6 +1810,11 @@ namespace GridTools * corresponding neighboring cells with points in unit coordinates are also * identified. * + * The parameter @p vertex_to_cells allows to accelerate the process of + * identifying the neighbors of a cell, by first precomputing a map from the + * vertex indices to the cells. Such data structure is, e.g., provided by + * GridTools::Cache::get_vertex_to_cell_map(). + * * This function is useful e.g. for discontinuous function spaces where, for * the case the given point `p` lies on a vertex, edge or face, several * cells might hold independent values of the solution that get combined in @@ -1843,7 +1848,10 @@ namespace GridTools const Point & p, const double tolerance, const std::pair::active_cell_iterator, - Point> & first_cell); + Point> & first_cell, + const std::vector< + std::set::active_cell_iterator>> + *vertex_to_cells = nullptr); /** * A variant of the previous function that internally calls one of the diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 6d84bbd0ed..039ffdcef1 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -6073,7 +6073,8 @@ namespace GridTools cache.get_triangulation(), point, tolerance, - first_cell); + first_cell, + &cache.get_vertex_to_cell_map()); if (enforce_unique_mapping) { diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 9ec9e0df4d..2913b3e64f 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -615,7 +615,10 @@ namespace GridTools const Point & p, const double tolerance, const std::pair::active_cell_iterator, - Point> & first_cell) + Point> & first_cell, + const std::vector< + std::set::active_cell_iterator>> + *vertex_to_cells) { std::vector< std::pair::active_cell_iterator, @@ -671,14 +674,19 @@ namespace GridTools 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; - std::vector::active_cell_iterator> - cells = find_cells_adjacent_to_vertex( - mesh, my_cell->vertex_index(local_vertex_index)); - for (const auto &cell : cells) - { + + 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); + + 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 @@ -712,21 +720,28 @@ namespace GridTools (vertex_indices[1].second << vertex_indices[1].first); for (unsigned int d = 0; d < 2; ++d) { - auto tentative_cells = find_cells_adjacent_to_vertex( - mesh, - my_cell->vertex_index(first_vertex + (d << free_direction))); - 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 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)); } } diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in index 644a19b533..9ca26ba67b 100644 --- a/source/grid/grid_tools_dof_handlers.inst.in +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -69,7 +69,8 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS; const Point &, const double, const std::pair> &); + Point> &, + const std::vector> *); template std::vector< std::pair