From: Niklas Fehn Date: Sun, 7 Jun 2020 17:22:33 +0000 (+0200) Subject: refactoring of GridTools::find_all_active_cells_around_point() X-Git-Tag: v9.3.0-rc1~1474^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10477%2Fhead;p=dealii.git refactoring of GridTools::find_all_active_cells_around_point() --- diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index a24b8c9a96..1a19ca6c1b 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1280,19 +1280,48 @@ namespace GridTools const std::vector &marked_vertices = {}); /** - * A variant of the previous find_active_cell_around_point() function that, - * instead of returning only the first matching cell, identifies all cells + * As compared to the functions above, this function identifies all cells * around a point for a given tolerance level `tolerance` in terms of unit - * coordinates. More precisely, whenever the point returned by - * find_active_cell_around_point() is within the given tolerance from the - * surface of the unit cell, all corresponding neighbors are also - * identified, including the location of the point in unit coordinates on - * any of these cells. + * coordinates. Given a first cell with reference coordinates as parameter + * @p first_cell, e.g. obtained by one of the functions above, all + * corresponding neighboring cells with points in unit coordinates are also + * identified. * * This function is useful e.g. for discontinuous function spaces where, for - * the case the given point `p` coincides with a vertex or an edge, several + * 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 * some way in a user code. + * + * This function is used as follows + * @code + * auto first_cell = GridTools::find_active_cell_around_point(...); + * auto all_cells = GridTools::find_all_active_cells_around_point(mapping, + * mesh, p, tolerance, first_cell); + * @endcode + */ + template class MeshType, int spacedim> +# ifndef _MSC_VER + std::vector::active_cell_iterator, + Point>> +# else + std::vector>::type, + Point>> +# endif + find_all_active_cells_around_point( + const Mapping & mapping, + const MeshType &mesh, + const Point & p, + const double tolerance, + const std::pair::active_cell_iterator, + Point> & first_cell); + + /** + * A variant of the previous function that internally calls one of the + * functions find_active_cell_around_point() to obtain a first cell, and + * subsequently adds all other cells by calling the function + * find_all_active_cells_around_point() above. */ template class MeshType, int spacedim> # ifndef _MSC_VER diff --git a/source/grid/grid_tools_dof_handlers.cc b/source/grid/grid_tools_dof_handlers.cc index 059bf89b50..c2bdd9cadf 100644 --- a/source/grid/grid_tools_dof_handlers.cc +++ b/source/grid/grid_tools_dof_handlers.cc @@ -592,133 +592,156 @@ namespace GridTools const double tolerance, const std::vector &marked_vertices) { - // first use the result of the single point function as a guess. In order - // not to make the other find_all_active_cells_around_point more expensive - // and avoid some additional logic there, we first start with one cell as - // given by that other function (that possibly goes through a larger set - // of cells) and later add a list of more cells as appropriate. - std::vector< - std::pair::active_cell_iterator, - Point>> - cells_and_points; try { - cells_and_points.push_back(find_active_cell_around_point_tolerance( - mapping, mesh, p, marked_vertices, tolerance)); + const auto cell_and_point = find_active_cell_around_point_tolerance( + mapping, mesh, p, marked_vertices, tolerance); + + return find_all_active_cells_around_point( + mapping, mesh, p, tolerance, cell_and_point); } catch (ExcPointNotFound &) {} - if (!cells_and_points.empty()) + return {}; + } + + + + template class MeshType, int spacedim> +#ifndef _MSC_VER + std::vector::active_cell_iterator, + Point>> +#else + std::vector>::type, + Point>> +#endif + find_all_active_cells_around_point( + const Mapping & mapping, + const MeshType &mesh, + const Point & p, + const double tolerance, + const std::pair::active_cell_iterator, + Point> & first_cell) + { + std::vector< + std::pair::active_cell_iterator, + Point>> + cells_and_points; + + // 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) { - // 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) + distance_to_center[d] = std::abs(unit_point[d] - 0.5); + if (distance_to_center[d] > 0.5 - tolerance) { - 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; - } + ++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)) + cells_to_add.push_back(my_cell->neighbor(neighbor_index)); + } + // 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; 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)) - cells_to_add.push_back(my_cell->neighbor(neighbor_index)); - } - // 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; - std::vector::active_cell_iterator> - 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); - } - // point on line in 3D: We cannot simply take the intersection between - // the two vertices of cels 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) + 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); + } + // 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) { - 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 (distance_to_center[d] > 0.5 - tolerance) { - 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; + 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; + } - 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 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) + { + auto tentative_cells = find_cells_adjacent_to_vertex( + mesh, + my_cell->vertex_index(first_vertex + (d << free_direction))); + for (const auto &cell : tentative_cells) { - 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); - } + 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 double original_distance_to_unit_cell = - GeometryInfo::distance_to_unit_cell(unit_point); - for (const auto &cell : cells_to_add) - { - if (cell != my_cell) - try - { - const Point p_unit = - mapping.transform_real_to_unit_cell(cell, p); - if (GeometryInfo::distance_to_unit_cell(p_unit) < - original_distance_to_unit_cell + tolerance) - cells_and_points.emplace_back(cell, p_unit); - } - catch (typename Mapping::ExcTransformationFailed &) - {} - } + const double original_distance_to_unit_cell = + GeometryInfo::distance_to_unit_cell(unit_point); + for (const auto &cell : cells_to_add) + { + if (cell != my_cell) + try + { + const Point p_unit = + mapping.transform_real_to_unit_cell(cell, p); + if (GeometryInfo::distance_to_unit_cell(p_unit) < + original_distance_to_unit_cell + tolerance) + cells_and_points.emplace_back(cell, p_unit); + } + catch (typename Mapping::ExcTransformationFailed &) + {} } + std::sort( cells_and_points.begin(), cells_and_points.end(), diff --git a/source/grid/grid_tools_dof_handlers.inst.in b/source/grid/grid_tools_dof_handlers.inst.in index 400f0c2984..30ca321b91 100644 --- a/source/grid/grid_tools_dof_handlers.inst.in +++ b/source/grid/grid_tools_dof_handlers.inst.in @@ -56,6 +56,19 @@ for (X : TRIANGULATION_AND_DOFHANDLERS; deal_II_dimension : DIMENSIONS; const Point &, const std::vector &); + template std::vector< + std::pair::type, + Point>> + find_all_active_cells_around_point( + const Mapping &, + const X &, + const Point &, + const double, + const std::pair> &); + template std::vector< std::pair +#include + +#include + +#include +#include + +#include "../tests.h" + + + +template +void +print_result(const Mapping & mapping, + const Triangulation &tria, + const Point p, + const double tolerance = 1.e-10) +{ + deallog << "Testing " << dim << "D with point " << p << " tolerance " + << tolerance << std::endl; + auto first_cell = GridTools::find_active_cell_around_point(mapping, tria, p); + auto c_p = GridTools::find_all_active_cells_around_point( + mapping, tria, p, tolerance, first_cell); + for (auto i : c_p) + deallog << "Cell: " << i.first->id() << " unit point " << i.second + << std::endl; + deallog << std::endl; +} + + + +template +void +test(unsigned int n_ref) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(n_ref); + MappingQGeneric mapping(3); + + Point p; + { + for (unsigned int d = 0; d < dim; ++d) + p[d] = 0.5; + print_result(mapping, tria, p); + } +} + + + +int +main() +{ + initlog(); + deallog << std::setprecision(10); + + test<3, 3>(1); +} diff --git a/tests/grid/find_all_active_cells_around_point_01_b.output b/tests/grid/find_all_active_cells_around_point_01_b.output new file mode 100644 index 0000000000..70cfe1449f --- /dev/null +++ b/tests/grid/find_all_active_cells_around_point_01_b.output @@ -0,0 +1,11 @@ + +DEAL::Testing 3D with point 0.5000000000 0.5000000000 0.5000000000 tolerance 1.000000000e-10 +DEAL::Cell: 0_1:0 unit point 1.000000000 1.000000000 1.000000000 +DEAL::Cell: 0_1:1 unit point 0.000000000 1.000000000 1.000000000 +DEAL::Cell: 0_1:2 unit point 1.000000000 0.000000000 1.000000000 +DEAL::Cell: 0_1:3 unit point 0.000000000 0.000000000 1.000000000 +DEAL::Cell: 0_1:4 unit point 1.000000000 1.000000000 0.000000000 +DEAL::Cell: 0_1:5 unit point 0.000000000 1.000000000 0.000000000 +DEAL::Cell: 0_1:6 unit point 1.000000000 0.000000000 0.000000000 +DEAL::Cell: 0_1:7 unit point 0.000000000 0.000000000 0.000000000 +DEAL::