--- /dev/null
+New: Mapping::get_bounding_box() returns the correct bounding box for mappings that do not preserve vertex
+locations.
+<br>
+(Luca Heltai, 2019/05/03)
* 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
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
+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(
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;
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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]
projects / dealii.git / commitdiff