From: Luca Heltai Date: Sat, 16 Jan 2021 21:19:16 +0000 (+0100) Subject: Return the box of support points instead of vertices. X-Git-Tag: v9.3.0-rc1~603^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5b874b72accd61340b2753bdd89a09668754c7c6;p=dealii.git Return the box of support points instead of vertices. --- diff --git a/doc/news/changes/minor/20210116LucaHeltai b/doc/news/changes/minor/20210116LucaHeltai new file mode 100644 index 0000000000..5f54d90c37 --- /dev/null +++ b/doc/news/changes/minor/20210116LucaHeltai @@ -0,0 +1,4 @@ +New: Mapping::get_bounding_box() now returns the bounding box of the support points, when the mapping +is of type MappingQ, MappingQGeneric, and MappingQEulerian. +
+(Luca Heltai, 2021/01/16) diff --git a/include/deal.II/fe/mapping.h b/include/deal.II/fe/mapping.h index f04ee0890a..b6d587aa04 100644 --- a/include/deal.II/fe/mapping.h +++ b/include/deal.II/fe/mapping.h @@ -375,11 +375,11 @@ public: * 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. Beware of the fact that - * for higher order mappings this bounding box is only an approximation of the - * true bounding box, since it does not take into account curved faces, and it - * may be smaller than the true bounding box. + * For linear mappings, this function returns the bounding box containing all + * the vertices of the cell, as returned by the get_vertices() method. For + * higher order mappings defined through support points, the bounding box is + * only guaranteed to contain all the support points, and it is, in general, + * only an approximation of the true bounding box, which may be larger. * * @param[in] cell The cell for which you want to compute the bounding box */ diff --git a/include/deal.II/fe/mapping_fe.h b/include/deal.II/fe/mapping_fe.h index 0865491ee0..24f4ad439c 100644 --- a/include/deal.II/fe/mapping_fe.h +++ b/include/deal.II/fe/mapping_fe.h @@ -79,6 +79,12 @@ public: unsigned int get_degree() const; + // for documentation, see the Mapping base class + virtual BoundingBox + get_bounding_box(const typename Triangulation::cell_iterator + &cell) const override; + + /** * Always returns @p true because the default implementation of functions in * this class preserves vertex locations. diff --git a/include/deal.II/fe/mapping_q.h b/include/deal.II/fe/mapping_q.h index 47ad2bf5d8..168191bce4 100644 --- a/include/deal.II/fe/mapping_q.h +++ b/include/deal.II/fe/mapping_q.h @@ -119,6 +119,11 @@ public: virtual bool preserves_vertex_locations() const override; + // for documentation, see the Mapping base class + virtual BoundingBox + get_bounding_box(const typename Triangulation::cell_iterator + &cell) const override; + /** * Transform the point @p p on the unit cell to the point @p p_real on the * real cell @p cell and returns @p p_real. diff --git a/include/deal.II/fe/mapping_q_eulerian.h b/include/deal.II/fe/mapping_q_eulerian.h index 10bd13b13f..bf4d090417 100644 --- a/include/deal.II/fe/mapping_q_eulerian.h +++ b/include/deal.II/fe/mapping_q_eulerian.h @@ -140,6 +140,11 @@ public: virtual bool preserves_vertex_locations() const override; + // for documentation, see the Mapping base class + virtual BoundingBox + get_bounding_box(const typename Triangulation::cell_iterator + &cell) const override; + /** * Exception which is thrown when the mapping is being evaluated at * non-active cell. diff --git a/include/deal.II/fe/mapping_q_generic.h b/include/deal.II/fe/mapping_q_generic.h index a3000b187b..8a2d6ddb0e 100644 --- a/include/deal.II/fe/mapping_q_generic.h +++ b/include/deal.II/fe/mapping_q_generic.h @@ -165,6 +165,11 @@ public: virtual bool preserves_vertex_locations() const override; + // for documentation, see the Mapping base class + virtual BoundingBox + get_bounding_box(const typename Triangulation::cell_iterator + &cell) const override; + /** * @name Mapping points between reference and real cells * @{ diff --git a/source/fe/mapping_fe.cc b/source/fe/mapping_fe.cc index 3f7148dbca..1aea54c729 100644 --- a/source/fe/mapping_fe.cc +++ b/source/fe/mapping_fe.cc @@ -2157,6 +2157,16 @@ MappingFE::compute_mapping_support_points( +template +BoundingBox +MappingFE::get_bounding_box( + const typename Triangulation::cell_iterator &cell) const +{ + return BoundingBox(this->compute_mapping_support_points(cell)); +} + + + //--------------------------- Explicit instantiations ----------------------- #include "mapping_fe.inst" diff --git a/source/fe/mapping_q.cc b/source/fe/mapping_q.cc index 69276944a3..38f412e293 100644 --- a/source/fe/mapping_q.cc +++ b/source/fe/mapping_q.cc @@ -544,6 +544,21 @@ MappingQ::clone() const +template +BoundingBox +MappingQ::get_bounding_box( + const typename Triangulation::cell_iterator &cell) const +{ + if (cell->has_boundary_lines() || use_mapping_q_on_all_cells || + (dim != spacedim)) + return BoundingBox( + qp_mapping->compute_mapping_support_points(cell)); + else + return BoundingBox(q1_mapping->get_vertices(cell)); +} + + + // explicit instantiations #include "mapping_q.inst" diff --git a/source/fe/mapping_q_eulerian.cc b/source/fe/mapping_q_eulerian.cc index d2990e8cb9..dcd0d57c4f 100644 --- a/source/fe/mapping_q_eulerian.cc +++ b/source/fe/mapping_q_eulerian.cc @@ -271,6 +271,18 @@ MappingQEulerian::fill_fe_values( +template +BoundingBox +MappingQEulerian::get_bounding_box( + const typename Triangulation::cell_iterator &cell) const +{ + return BoundingBox( + dynamic_cast(*this->qp_mapping) + .compute_mapping_support_points(cell)); +} + + + // explicit instantiations #include "mapping_q_eulerian.inst" diff --git a/source/fe/mapping_q_generic.cc b/source/fe/mapping_q_generic.cc index 4559fb6c9c..73c1678654 100644 --- a/source/fe/mapping_q_generic.cc +++ b/source/fe/mapping_q_generic.cc @@ -1643,6 +1643,16 @@ MappingQGeneric::compute_mapping_support_points( +template +BoundingBox +MappingQGeneric::get_bounding_box( + const typename Triangulation::cell_iterator &cell) const +{ + return BoundingBox(this->compute_mapping_support_points(cell)); +} + + + //--------------------------- Explicit instantiations ----------------------- #include "mapping_q_generic.inst" diff --git a/tests/mappings/mapping_get_bounding_box_01.cc b/tests/mappings/mapping_get_bounding_box_01.cc new file mode 100644 index 0000000000..6c495f32ea --- /dev/null +++ b/tests/mappings/mapping_get_bounding_box_01.cc @@ -0,0 +1,80 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2006 - 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 that bounding boxes are created correctly with high order mappings + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +#include + +#include "../tests.h" + + +template +void +test_bounding_boxes(const Mapping &mapping, + const unsigned int degree) +{ + deallog << "Testing " << boost::core::demangle(typeid(mapping).name()) << "(" + << degree << ")" << std::endl; + + Triangulation triangulation; + GridGenerator::hyper_ball(triangulation); + GridTools::Cache cache(triangulation, mapping); + + const auto boxes = cache.get_cell_bounding_boxes_rtree(); + + { + std::string fname = "boxes_" + std::to_string(dim) + "_" + + std::to_string(spacedim) + "_" + + std::to_string(degree) + ".vtk"; + std::ofstream ofile(fname); + BoundingBoxDataOut data_out; + DataOutBase::VtkFlags flags; + flags.print_date_and_time = false; + data_out.set_flags(flags); + data_out.build_patches(boxes); + data_out.write_vtk(ofile); + cat_file(fname.c_str()); + } +} + + +int +main() +{ + initlog(); + { + unsigned int degree = 2; + MappingQ<2> mapping(degree); + test_bounding_boxes(mapping, degree); + } + { + unsigned int degree = 2; + MappingQGeneric<2> mapping(degree); + test_bounding_boxes(mapping, degree); + } +} diff --git a/tests/mappings/mapping_get_bounding_box_01.output b/tests/mappings/mapping_get_bounding_box_01.output new file mode 100644 index 0000000000..6ceb716d34 --- /dev/null +++ b/tests/mappings/mapping_get_bounding_box_01.output @@ -0,0 +1,79 @@ + +DEAL::Testing dealii::MappingQ<2, 2>(2) +# vtk DataFile Version 3.0 +#This file was generated by the deal.II library. +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 20 double +-0.707107 -1 0 +0.707107 -1 0 +-0.707107 -0.292893 0 +0.707107 -0.292893 0 +-1 -0.707107 0 +-0.292893 -0.707107 0 +-1 0.707107 0 +-0.292893 0.707107 0 +-0.292893 -0.292893 0 +0.292893 -0.292893 0 +-0.292893 0.292893 0 +0.292893 0.292893 0 +0.292893 -0.707107 0 +1 -0.707107 0 +0.292893 0.707107 0 +1 0.707107 0 +-0.707107 0.292893 0 +0.707107 0.292893 0 +-0.707107 1 0 +0.707107 1 0 + +CELLS 5 25 +4 0 1 3 2 +4 4 5 7 6 +4 8 9 11 10 +4 12 13 15 14 +4 16 17 19 18 + +CELL_TYPES 5 + 9 9 9 9 9 +POINT_DATA 20 + +DEAL::Testing dealii::MappingQGeneric<2, 2>(2) +# vtk DataFile Version 3.0 +#This file was generated by the deal.II library. +ASCII +DATASET UNSTRUCTURED_GRID + +POINTS 20 double +-0.707107 -1 0 +0.707107 -1 0 +-0.707107 -0.292893 0 +0.707107 -0.292893 0 +-1 -0.707107 0 +-0.292893 -0.707107 0 +-1 0.707107 0 +-0.292893 0.707107 0 +-0.292893 -0.292893 0 +0.292893 -0.292893 0 +-0.292893 0.292893 0 +0.292893 0.292893 0 +0.292893 -0.707107 0 +1 -0.707107 0 +0.292893 0.707107 0 +1 0.707107 0 +-0.707107 0.292893 0 +0.707107 0.292893 0 +-0.707107 1 0 +0.707107 1 0 + +CELLS 5 25 +4 0 1 3 2 +4 4 5 7 6 +4 8 9 11 10 +4 12 13 15 14 +4 16 17 19 18 + +CELL_TYPES 5 + 9 9 9 9 9 +POINT_DATA 20 +