From: Luca Heltai <luca.heltai@sissa.it>
Date: Thu, 11 Apr 2019 13:53:58 +0000 (+0200)
Subject: Make sure rtree works with mappings that do not preserve vertices.
X-Git-Tag: v9.1.0-rc1~199^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F7915%2Fhead;p=dealii.git

Make sure rtree works with mappings that do not preserve vertices.
---

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. 
+<br>
+(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<std::pair<
           BoundingBox<spacedim>,
           typename Triangulation<dim, spacedim>::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<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);
+              }
+          }
         cell_bounding_boxes_rtree = pack_rtree(boxes);
       }
     return cell_bounding_boxes_rtree;
   }
 
 
-
 #ifdef DEAL_II_WITH_NANOFLANN
   template <int dim, int spacedim>
   const KDTree<spacedim> &
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 <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/grid_tools_cache.h>
+#include <deal.II/grid/tria.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 = 2.0;
+
+  MappingQEulerian<dim> euler(2, dof_handler, displacements);
+
+  GridTools::Cache<dim> cache0(triangulation);
+  GridTools::Cache<dim> 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