From cd6d59cae1ef7849c2bb639ac5597d4bf5dff0d4 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 3 May 2019 19:53:06 +0200 Subject: [PATCH] Mapping::get_bounding_box() --- doc/news/changes/minor/20190503LucaHeltai | 4 ++ include/deal.II/fe/mapping.h | 22 +++++- source/fe/mapping.cc | 23 ++++++ source/grid/grid_tools_cache.cc | 29 ++------ tests/mappings/mapping_q_eulerian_13.cc | 77 +++++++++++++++++++++ tests/mappings/mapping_q_eulerian_13.output | 7 ++ 6 files changed, 137 insertions(+), 25 deletions(-) create mode 100644 doc/news/changes/minor/20190503LucaHeltai create mode 100644 tests/mappings/mapping_q_eulerian_13.cc create mode 100644 tests/mappings/mapping_q_eulerian_13.output diff --git a/doc/news/changes/minor/20190503LucaHeltai b/doc/news/changes/minor/20190503LucaHeltai new file mode 100644 index 0000000000..1590818e35 --- /dev/null +++ b/doc/news/changes/minor/20190503LucaHeltai @@ -0,0 +1,4 @@ +New: Mapping::get_bounding_box() returns the correct bounding box for mappings that do not preserve vertex +locations. +
+(Luca Heltai, 2019/05/03) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index 544ba862e2..0c91e615b0 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -353,7 +353,7 @@ public: * will be the average of the vertex locations, as returned by the * get_vertices() method. * - * @param[in] cell The cell for which you want to compute the center + * @param[in] cell The cell for which you want to compute the center; * @param[in] map_center_of_reference_cell A flag that switches the algorithm * for the computation of the cell center from * transform_unit_to_real_cell() applied to the center of the reference cell @@ -365,6 +365,26 @@ public: get_center(const typename Triangulation::cell_iterator &cell, const bool map_center_of_reference_cell = true) const; + /** + * Return the bounding box of a mapped cell. + * + * If you are using a (bi-,tri-)linear mapping that preserves vertex + * locations, this function simply returns the value also produced by + * `cell->bounding_box()`. However, there are also mappings that add + * displacements or choose completely different locations, e.g., + * MappingQEulerian, MappingQ1Eulerian, or MappingFEField. + * + * This function returns the bounding box containing all the vertices of the + * cell, as returned by the get_vertices() method. + * + * @param[in] cell The cell for which you want to compute the bounding box. + * + * @author Luca Heltai, 2019. + */ + virtual BoundingBox + get_bounding_box( + const typename Triangulation::cell_iterator &cell) const; + /** * Return whether the mapping preserves vertex locations. In other words, * this function returns whether the mapped location of the reference cell diff --git a/source/fe/mapping.cc b/source/fe/mapping.cc index 97a7a38f23..45a1371568 100644 --- a/source/fe/mapping.cc +++ b/source/fe/mapping.cc @@ -61,6 +61,29 @@ Mapping::get_center( +template +BoundingBox +Mapping::get_bounding_box( + const typename Triangulation::cell_iterator &cell) const +{ + if (preserves_vertex_locations()) + return cell->bounding_box(); + else + { + const auto vertices = 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]); + } + return BoundingBox({p0, p1}); + } +} + + + template Point Mapping::project_real_point_to_unit_point_on_face( diff --git a/source/grid/grid_tools_cache.cc b/source/grid/grid_tools_cache.cc index b0e1be61cf..4ed0b878c7 100644 --- a/source/grid/grid_tools_cache.cc +++ b/source/grid/grid_tools_cache.cc @@ -131,30 +131,11 @@ namespace GridTools std::vector, typename Triangulation::active_cell_iterator>> - 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); - } - } + boxes(tria->n_active_cells()); + unsigned int i = 0; + for (const auto &cell : tria->active_cell_iterators()) + boxes[i++] = std::make_pair(mapping->get_bounding_box(cell), cell); + cell_bounding_boxes_rtree = pack_rtree(boxes); } return cell_bounding_boxes_rtree; diff --git a/tests/mappings/mapping_q_eulerian_13.cc b/tests/mappings/mapping_q_eulerian_13.cc new file mode 100644 index 0000000000..ab005101cc --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_13.cc @@ -0,0 +1,77 @@ +// --------------------------------------------------------------------- +// +// 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 MappingQEulerian knows how to compute bounding boxes of cells + +#include + +#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 = 1.0; + + MappingQEulerian euler(2, dof_handler, displacements); + + // now the actual test + for (const auto &cell : dof_handler.active_cell_iterators()) + { + const auto p1 = cell->bounding_box().get_boundary_points(); + const auto p2 = euler.get_bounding_box(cell).get_boundary_points(); + deallog << "BBox: [" << p1.first << ", " << p1.second + << "], with mapping [" << p2.first << ", " << p2.second << "]" + << std::endl; + } +} + +int +main() +{ + initlog(); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/mappings/mapping_q_eulerian_13.output b/tests/mappings/mapping_q_eulerian_13.output new file mode 100644 index 0000000000..cbba8f0264 --- /dev/null +++ b/tests/mappings/mapping_q_eulerian_13.output @@ -0,0 +1,7 @@ + +DEAL::dim=1 +DEAL::BBox: [-1.00000, 1.00000], with mapping [0.00000, 2.00000] +DEAL::dim=2 +DEAL::BBox: [-1.00000 -1.00000, 1.00000 1.00000], with mapping [0.00000 0.00000, 2.00000 2.00000] +DEAL::dim=3 +DEAL::BBox: [-1.00000 -1.00000 -1.00000, 1.00000 1.00000 1.00000], with mapping [0.00000 0.00000 0.00000, 2.00000 2.00000 2.00000] -- 2.39.5