From 5f2ba3ddf5c46a23a617751fbb697165d3d5c68d Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Sun, 10 May 2020 02:21:05 +0200 Subject: [PATCH] Add mapping argument to find_closest_vertex_of_cell(). Closes #7922. --- doc/news/changes/minor/20200510LucaHeltai | 4 + include/deal.II/grid/grid_tools.h | 11 ++- source/grid/grid_tools.cc | 16 ++-- source/grid/grid_tools.inst.in | 3 +- .../grid/grid_tools_find_closest_vertex_02.cc | 78 +++++++++++++++++++ .../grid_tools_find_closest_vertex_02.output | 21 +++++ 6 files changed, 123 insertions(+), 10 deletions(-) create mode 100644 doc/news/changes/minor/20200510LucaHeltai create mode 100644 tests/grid/grid_tools_find_closest_vertex_02.cc create mode 100644 tests/grid/grid_tools_find_closest_vertex_02.output diff --git a/doc/news/changes/minor/20200510LucaHeltai b/doc/news/changes/minor/20200510LucaHeltai new file mode 100644 index 0000000000..d4b9fc7bc2 --- /dev/null +++ b/doc/news/changes/minor/20200510LucaHeltai @@ -0,0 +1,4 @@ +Fixed: GridTools::find_closest_vertex_of_cell() now takes an optional mapping argument, to take into +account of all those cases where a mapping may modify the location of the cell vertices. +
+(Luca Heltai, 2020/05/10) diff --git a/include/deal.II/grid/grid_tools.h b/include/deal.II/grid/grid_tools.h index b3935cd15a..9a8016766c 100644 --- a/include/deal.II/grid/grid_tools.h +++ b/include/deal.II/grid/grid_tools.h @@ -1824,15 +1824,20 @@ namespace GridTools /** * Return the local vertex index of cell @p cell that is closest to - * the given location @p position. + * the given location @p position. The location of the vertices is extracted + * from the (optional) @p mapping argument, to guarantee that the correct + * answer is returned when the underlying mapping modifies the position of the + * vertices. * - * @author Rene Gassmoeller, Luca Heltai, 2017. + * @author Rene Gassmoeller, Luca Heltai, 2017, 2020. */ template unsigned int find_closest_vertex_of_cell( const typename Triangulation::active_cell_iterator &cell, - const Point &position); + const Point & position, + const Mapping & mapping = + StaticMappingQ1::mapping); /** * Compute a globally unique index for each vertex and hanging node diff --git a/source/grid/grid_tools.cc b/source/grid/grid_tools.cc index 685c975582..6a78bf705e 100644 --- a/source/grid/grid_tools.cc +++ b/source/grid/grid_tools.cc @@ -1604,9 +1604,12 @@ namespace GridTools // invalid cell, then query for the closest global vertex if (current_cell.state() == IteratorState::valid) { + const auto cell_vertices = mapping.get_vertices(current_cell); const unsigned int closest_vertex = - find_closest_vertex_of_cell(current_cell, p); - vertex_to_point = p - current_cell->vertex(closest_vertex); + find_closest_vertex_of_cell(current_cell, + p, + mapping); + vertex_to_point = p - cell_vertices[closest_vertex]; closest_vertex_index = current_cell->vertex_index(closest_vertex); } else @@ -1736,15 +1739,16 @@ namespace GridTools unsigned int find_closest_vertex_of_cell( const typename Triangulation::active_cell_iterator &cell, - const Point & position) + const Point & position, + const Mapping & mapping) { - double minimum_distance = position.distance_square(cell->vertex(0)); + const auto vertices = mapping.get_vertices(cell); + double minimum_distance = position.distance_square(vertices[0]); unsigned int closest_vertex = 0; for (unsigned int v = 1; v < GeometryInfo::vertices_per_cell; ++v) { - const double vertex_distance = - position.distance_square(cell->vertex(v)); + const double vertex_distance = position.distance_square(vertices[v]); if (vertex_distance < minimum_distance) { closest_vertex = v; diff --git a/source/grid/grid_tools.inst.in b/source/grid/grid_tools.inst.in index c080375436..77b083198d 100644 --- a/source/grid/grid_tools.inst.in +++ b/source/grid/grid_tools.inst.in @@ -68,7 +68,8 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const typename Triangulation< deal_II_dimension, deal_II_space_dimension>::active_cell_iterator &, - const Point &); + const Point &, + const Mapping &); template std::map compute_local_to_global_vertex_index_map( diff --git a/tests/grid/grid_tools_find_closest_vertex_02.cc b/tests/grid/grid_tools_find_closest_vertex_02.cc new file mode 100644 index 0000000000..28ed4793c6 --- /dev/null +++ b/tests/grid/grid_tools_find_closest_vertex_02.cc @@ -0,0 +1,78 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 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. +// +// --------------------------------------------------------------------- + +// Check find_closest_vertex_of_cell when a mapping is present + +#include + +#include +#include +#include + +#include +#include +#include +#include + +#include "../tests.h" + + +template +void +test(const Point &p, double displacement) +{ + deallog << "dim: " << dim << ", spacedim: " << spacedim + << " with displacement " << displacement << std::endl; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(); + + // Vector fe + FESystem fe{FE_Q(1), spacedim}; + DoFHandler dh(tria); + dh.distribute_dofs(fe); + + // Displacement vector. + Vector disp(dh.n_dofs()); + for (unsigned int i = 0; i < dh.n_dofs(); ++i) + disp[i] = displacement; + + // mapping + MappingQEulerian, spacedim> mapping(1, dh, disp); + + for (const auto cell : tria.active_cell_iterators()) + { + const auto i = + GridTools::find_closest_vertex_of_cell(cell, p, mapping); + const auto &v = mapping.get_vertices(cell); + + deallog << "Closest vertex of cell " << cell << " to " << p + << " is cell->vertex(" << i << ") :" << v[i] << std::endl; + } +}; + + +int +main() +{ + initlog(); + test<2, 2>(Point<2>(.45, .45), 0); + test<2, 2>(Point<2>(.45, .45), 0.5); + + test<2, 3>(Point<3>(.9, .9, .1), 0.0); + test<2, 3>(Point<3>(.9, .9, .1), 1.0); + return 0; +} diff --git a/tests/grid/grid_tools_find_closest_vertex_02.output b/tests/grid/grid_tools_find_closest_vertex_02.output new file mode 100644 index 0000000000..3839578fab --- /dev/null +++ b/tests/grid/grid_tools_find_closest_vertex_02.output @@ -0,0 +1,21 @@ + +DEAL::dim: 2, spacedim: 2 with displacement 0.00000 +DEAL::Closest vertex of cell 1.0 to 0.450000 0.450000 is cell->vertex(3) :0.500000 0.500000 +DEAL::Closest vertex of cell 1.1 to 0.450000 0.450000 is cell->vertex(2) :0.500000 0.500000 +DEAL::Closest vertex of cell 1.2 to 0.450000 0.450000 is cell->vertex(1) :0.500000 0.500000 +DEAL::Closest vertex of cell 1.3 to 0.450000 0.450000 is cell->vertex(0) :0.500000 0.500000 +DEAL::dim: 2, spacedim: 2 with displacement 0.500000 +DEAL::Closest vertex of cell 1.0 to 0.450000 0.450000 is cell->vertex(0) :0.500000 0.500000 +DEAL::Closest vertex of cell 1.1 to 0.450000 0.450000 is cell->vertex(0) :1.00000 0.500000 +DEAL::Closest vertex of cell 1.2 to 0.450000 0.450000 is cell->vertex(0) :0.500000 1.00000 +DEAL::Closest vertex of cell 1.3 to 0.450000 0.450000 is cell->vertex(0) :1.00000 1.00000 +DEAL::dim: 2, spacedim: 3 with displacement 0.00000 +DEAL::Closest vertex of cell 1.0 to 0.900000 0.900000 0.100000 is cell->vertex(3) :0.500000 0.500000 0.00000 +DEAL::Closest vertex of cell 1.1 to 0.900000 0.900000 0.100000 is cell->vertex(3) :1.00000 0.500000 0.00000 +DEAL::Closest vertex of cell 1.2 to 0.900000 0.900000 0.100000 is cell->vertex(3) :0.500000 1.00000 0.00000 +DEAL::Closest vertex of cell 1.3 to 0.900000 0.900000 0.100000 is cell->vertex(3) :1.00000 1.00000 0.00000 +DEAL::dim: 2, spacedim: 3 with displacement 1.00000 +DEAL::Closest vertex of cell 1.0 to 0.900000 0.900000 0.100000 is cell->vertex(0) :1.00000 1.00000 1.00000 +DEAL::Closest vertex of cell 1.1 to 0.900000 0.900000 0.100000 is cell->vertex(0) :1.50000 1.00000 1.00000 +DEAL::Closest vertex of cell 1.2 to 0.900000 0.900000 0.100000 is cell->vertex(0) :1.00000 1.50000 1.00000 +DEAL::Closest vertex of cell 1.3 to 0.900000 0.900000 0.100000 is cell->vertex(0) :1.50000 1.50000 1.00000 -- 2.39.5