From 77fb2394f658eac51d0be277f583ede16b6c2993 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Thu, 11 Apr 2019 15:53:58 +0200 Subject: [PATCH] Make sure rtree works with mappings that do not preserve vertices. --- doc/news/changes/minor/20190411LucaHeltai-c | 3 + source/grid/grid_tools_cache.cc | 29 +++++-- tests/grid/grid_tools_cache_05.cc | 83 +++++++++++++++++++++ tests/grid/grid_tools_cache_05.output | 10 +++ 4 files changed, 120 insertions(+), 5 deletions(-) create mode 100644 doc/news/changes/minor/20190411LucaHeltai-c create mode 100644 tests/grid/grid_tools_cache_05.cc create mode 100644 tests/grid/grid_tools_cache_05.output diff --git a/doc/news/changes/minor/20190411LucaHeltai-c b/doc/news/changes/minor/20190411LucaHeltai-c new file mode 100644 index 0000000000..e8a86eb0ff --- /dev/null +++ b/doc/news/changes/minor/20190411LucaHeltai-c @@ -0,0 +1,3 @@ +Fixed: Make sure that GridTools::Cache::get_cell_bounding_boxes_rtree() honors mappings. +
+(Luca Heltai, 2019/04/11) diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc index d1098c6986..b0e1be61cf 100644 --- a/source/grid/grid_tools_cache.cc +++ b/source/grid/grid_tools_cache.cc @@ -131,17 +131,36 @@ namespace GridTools std::vector, typename Triangulation::active_cell_iterator>> - boxes(tria->n_active_cells()); - unsigned int i = 0; - for (const auto &cell : tria->active_cell_iterators()) - boxes[i++] = std::make_pair(cell->bounding_box(), cell); + boxes(tria->n_active_cells()); + if (mapping->preserves_vertex_locations()) + { + unsigned int i = 0; + for (const auto &cell : tria->active_cell_iterators()) + boxes[i++] = std::make_pair(cell->bounding_box(), cell); + } + else + { + unsigned int i = 0; + for (const auto &cell : tria->active_cell_iterators()) + { + const auto vertices = mapping->get_vertices(cell); + Point p0 = vertices[0], p1 = vertices[0]; + for (unsigned int j = 1; j < vertices.size(); ++j) + for (unsigned int d = 0; d < spacedim; ++d) + { + p0[d] = std::min(p0[d], vertices[j][d]); + p1[d] = std::max(p1[d], vertices[j][d]); + } + boxes[i++] = + std::make_pair(BoundingBox({p0, p1}), cell); + } + } cell_bounding_boxes_rtree = pack_rtree(boxes); } return cell_bounding_boxes_rtree; } - #ifdef DEAL_II_WITH_NANOFLANN template const KDTree & diff --git a/tests/grid/grid_tools_cache_05.cc b/tests/grid/grid_tools_cache_05.cc new file mode 100644 index 0000000000..c4aef67106 --- /dev/null +++ b/tests/grid/grid_tools_cache_05.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 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 that GridTools::Cache::get_cell_bounding_boxes_rtree() +// works in conjunction with MappingQEulerian + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +using namespace dealii; + + +template +void +test() +{ + deallog << "dim=" << dim << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -1, 1); + + FESystem fe(FE_Q(2), dim); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + Vector displacements(dof_handler.n_dofs()); + + displacements = 2.0; + + MappingQEulerian euler(2, dof_handler, displacements); + + GridTools::Cache cache0(triangulation); + GridTools::Cache cache1(triangulation, euler); + + auto &tree0 = cache0.get_cell_bounding_boxes_rtree(); + auto &tree1 = cache1.get_cell_bounding_boxes_rtree(); + + const auto box0 = tree0.begin()->first; + const auto box1 = tree1.begin()->first; + + deallog << "No Mapping : " << box0.get_boundary_points().first << ", " + << box0.get_boundary_points().second << std::endl + << "MappingQEulerian: " << box1.get_boundary_points().first << ", " + << box1.get_boundary_points().second << std::endl; +} + + + +int +main() +{ + initlog(); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/grid/grid_tools_cache_05.output b/tests/grid/grid_tools_cache_05.output new file mode 100644 index 0000000000..a12320a83c --- /dev/null +++ b/tests/grid/grid_tools_cache_05.output @@ -0,0 +1,10 @@ + +DEAL::dim=1 +DEAL::No Mapping : -1.00000, 1.00000 +DEAL::MappingQEulerian: 1.00000, 3.00000 +DEAL::dim=2 +DEAL::No Mapping : -1.00000 -1.00000, 1.00000 1.00000 +DEAL::MappingQEulerian: 1.00000 1.00000, 3.00000 3.00000 +DEAL::dim=3 +DEAL::No Mapping : -1.00000 -1.00000 -1.00000, 1.00000 1.00000 1.00000 +DEAL::MappingQEulerian: 1.00000 1.00000 1.00000, 3.00000 3.00000 3.00000 -- 2.39.5