From: Luca Heltai <luca.heltai@sissa.it>
Date: Fri, 3 May 2019 17:53:06 +0000 (+0200)
Subject: Mapping::get_bounding_box()
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fheads%2Fmapping-bbox;p=dealii.git

Mapping::get_bounding_box()
---

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.
+<br>
+(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<dim, spacedim>::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<spacedim>
+  get_bounding_box(
+    const typename Triangulation<dim, spacedim>::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<dim, spacedim>::get_center(
 
 
 
+template <int dim, int spacedim>
+BoundingBox<spacedim>
+Mapping<dim, spacedim>::get_bounding_box(
+  const typename Triangulation<dim, spacedim>::cell_iterator &cell) const
+{
+  if (preserves_vertex_locations())
+    return cell->bounding_box();
+  else
+    {
+      const auto      vertices = get_vertices(cell);
+      Point<spacedim> 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<spacedim>({p0, p1});
+    }
+}
+
+
+
 template <int dim, int spacedim>
 Point<dim - 1>
 Mapping<dim, spacedim>::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<std::pair<
           BoundingBox<spacedim>,
           typename Triangulation<dim, spacedim>::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<spacedim> 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<spacedim>({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 <deal.II/base/function.h>
+
+#include <deal.II/dofs/dof_accessor.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_system.h>
+#include <deal.II/fe/mapping_q_eulerian.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/numerics/vector_tools.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test()
+{
+  deallog << "dim=" << dim << std::endl;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube(triangulation, -1, 1);
+
+  FESystem<dim>   fe(FE_Q<dim>(2), dim);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(fe);
+
+  Vector<double> displacements(dof_handler.n_dofs());
+  displacements = 1.0;
+
+  MappingQEulerian<dim> 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]