--- /dev/null
+New: GridTools::Cache::get_locally_owned_cell_bounding_boxes_rtree() extracts
+a tree of bounding boxes covering the locally owned cells of a triangulation.
+This can be used in conjunction with GridTools::Cache::get_covering_rtree() to
+make approximate geometrical queries on who owns what spacial region.
+<br>
+(Luca Heltai, 2020/05/26)
get_used_vertices_rtree() const;
/**
- * Return the cached RTree object of the cell bounging boxes, constructed
- * using the active cell iterators of the stored triangulation.
+ * Return the cached RTree object of the cell bounding boxes, constructed
+ * using the active cell iterators of the stored triangulation. For
+ * parallel::distributed::Triangulation objects, this function will return
+ * also the bounding boxes of ghost and artificial cells.
*/
const RTree<
std::pair<BoundingBox<spacedim>,
typename Triangulation<dim, spacedim>::active_cell_iterator>> &
get_cell_bounding_boxes_rtree() const;
+ /**
+ * Return the cached RTree object of bounding boxes containing locally owned
+ * active cells, constructed using the active cell iterators of the stored
+ * triangulation.
+ *
+ * In contrast to the previous function, this function builds the RTree
+ * using only the locally owned cells, i.e., not including ghost or
+ * artificial cells. The two functions return the same result in serial
+ * computations.
+ */
+ const RTree<
+ std::pair<BoundingBox<spacedim>,
+ typename Triangulation<dim, spacedim>::active_cell_iterator>> &
+ get_locally_owned_cell_bounding_boxes_rtree() const;
+
/**
* Return a reference to the stored triangulation.
*/
* object -- an Rtree -- are pairs of bounding boxes denoting
* areas that cover all or parts of the local portion of a
* parallel triangulation, and an unsigned int representing
- * the process or subdomain that owns these cells.
+ * the process or subdomain that owns the cells falling within the given
+ * bounding box.
+ *
* Given a point on a parallel::TriangulationBase, this tree
* allows to identify one, or few candidate processes, for
* which the point lies on a locally owned cell.
* Therefore this function must be called by all processes
* at the same time.
*
- * While each box may only cover part of a process's locally
- * owned part of the triangulation, the boxes associated with
- * each process jointly cover the entire local portion.
+ * The local bounding boxes are constructed by extracting the
+ * specified @p level from the rtree object returned by the
+ * get_locally_owned_cell_bounding_boxes_rtree(). Notice that @p level
+ * in this context does not refer to the level of the triangulation, but
+ * refer to the level of the RTree object (see, e.g.,
+ * https://en.wikipedia.org/wiki/R-tree).
*/
const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &
- get_covering_rtree() const;
+ get_covering_rtree(const unsigned int level = 0) const;
private:
/**
vertex_to_cell_centers;
/**
- * An rtree object covering the whole mesh.
+ * A collection of rtree objects covering the whole mesh.
+ *
+ * Each entry of the map is constructed from the function
+ * extract_rtree_level() applied to the rtree returned by the function
+ * get_locally_owned_cell_bounding_boxes_rtree(), with the specified level
+ * in input.
*/
- mutable RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
+ mutable std::map<unsigned int,
+ RTree<std::pair<BoundingBox<spacedim>, unsigned int>>>
covering_rtree;
/**
typename Triangulation<dim, spacedim>::active_cell_iterator>>
cell_bounding_boxes_rtree;
+ /**
+ * Store an RTree object, containing the bounding boxes of the locally owned
+ * cells of the triangulation.
+ */
+ mutable RTree<
+ std::pair<BoundingBox<spacedim>,
+ typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ locally_owned_cell_bounding_boxes_rtree;
+
/**
* Storage for the status of the triangulation signal.
*/
*/
update_covering_rtree = 0x040,
+ /**
+ * Update an RTree of locally owned cell bounding boxes.
+ */
+ update_locally_owned_cell_bounding_boxes_rtree = 0x080,
+
/**
* Update all objects.
*/
template <int dim, int spacedim>
- const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &
- Cache<dim, spacedim>::get_covering_rtree() const
+ const RTree<
+ std::pair<BoundingBox<spacedim>,
+ typename Triangulation<dim, spacedim>::active_cell_iterator>> &
+ Cache<dim, spacedim>::get_locally_owned_cell_bounding_boxes_rtree() const
{
- if (update_flags & update_covering_rtree)
+ if (update_flags & update_locally_owned_cell_bounding_boxes_rtree)
{
- // TO DO: the locally owned portion of the domain
- // might consists of more separate pieces: a single
- // bounding box might not always be the best choice.
+ std::vector<std::pair<
+ BoundingBox<spacedim>,
+ typename Triangulation<dim, spacedim>::active_cell_iterator>>
+ boxes;
+ boxes.reserve(tria->n_active_cells());
+ for (const auto &cell : tria->active_cell_iterators())
+ if (cell->is_locally_owned())
+ boxes.emplace_back(
+ std::make_pair(mapping->get_bounding_box(cell), cell));
- // Note: we create the local variable ptr here, because gcc 4.8.5
- // fails to compile if we pass the variable directly.
- IteratorFilters::LocallyOwnedCell locally_owned_cell_predicate;
- std::function<bool(
- const typename Triangulation<dim, spacedim>::active_cell_iterator &)>
- ptr(locally_owned_cell_predicate);
+ locally_owned_cell_bounding_boxes_rtree = pack_rtree(boxes);
+ update_flags =
+ update_flags & ~update_locally_owned_cell_bounding_boxes_rtree;
+ }
+ return locally_owned_cell_bounding_boxes_rtree;
+ }
- const BoundingBox<spacedim> bbox =
- GridTools::compute_bounding_box(get_triangulation(), ptr);
- std::vector<BoundingBox<spacedim>> bbox_v(1, bbox);
+ template <int dim, int spacedim>
+ const RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &
+ Cache<dim, spacedim>::get_covering_rtree(const unsigned int level) const
+ {
+ if (update_flags & update_covering_rtree ||
+ covering_rtree.find(level) == covering_rtree.end())
+ {
+ const auto boxes =
+ extract_rtree_level(get_locally_owned_cell_bounding_boxes_rtree(),
+ level);
+
if (const auto tria_mpi =
dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
&(*tria)))
{
- covering_rtree = GridTools::build_global_description_tree(
- bbox_v, tria_mpi->get_communicator());
+ covering_rtree[level] = GridTools::build_global_description_tree(
+ boxes, tria_mpi->get_communicator());
}
else
{
- covering_rtree =
- GridTools::build_global_description_tree(bbox_v, MPI_COMM_SELF);
+ covering_rtree[level] =
+ GridTools::build_global_description_tree(boxes, MPI_COMM_SELF);
}
update_flags = update_flags & ~update_covering_rtree;
}
- return covering_rtree;
+ return covering_rtree[level];
}
#include "grid_tools_cache.inst"
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test that Cache::get_locally_owned_cell_bounding_boxes_rtree() works
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/boost_adaptors/bounding_box.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include "../tests.h"
+
+namespace bgi = boost::geometry::index;
+
+template <int dim, int spacedim = dim>
+void
+test(unsigned int ref)
+{
+ const MPI_Comm &mpi_communicator = MPI_COMM_WORLD;
+
+ parallel::distributed::Triangulation<dim, spacedim> tria(mpi_communicator);
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(ref);
+
+ GridTools::Cache<dim, spacedim> cache(tria);
+
+ const auto &tree = cache.get_locally_owned_cell_bounding_boxes_rtree();
+ const auto box = random_box<dim>();
+
+ deallog << "Query box: " << box.get_boundary_points().first << "; "
+ << box.get_boundary_points().second << std::endl;
+
+ for (const auto p : tree | bgi::adaptors::queried(bgi::intersects(box)))
+ deallog << "Cell " << p.second << " intersects box with "
+ << p.first.get_boundary_points().first << "; "
+ << p.first.get_boundary_points().second << std::endl;
+}
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll log;
+
+ test<2>(3);
+
+ return 0;
+}
--- /dev/null
+
+DEAL:0::Query box: 0.783099 0.394383; 0.840188 0.798440
+DEAL:0::Cell 3.212 intersects box with 0.753761 0.349513; 0.923880 0.555570
+DEAL:0::Cell 3.213 intersects box with 0.655330 0.507759; 0.831470 0.707107
+DEAL:0::Cell 3.214 intersects box with 0.676052 0.316342; 0.836215 0.507759
--- /dev/null
+
+DEAL:0::Query box: 0.783099 0.394383; 0.840188 0.798440
+
+DEAL:1::Query box: 0.783099 0.394383; 0.840188 0.798440
+DEAL:1::Cell 3.114 intersects box with 0.676052 0.316342; 0.836215 0.507759
+DEAL:1::Cell 3.112 intersects box with 0.753761 0.349513; 0.923880 0.555570
+DEAL:1::Cell 3.113 intersects box with 0.655330 0.507759; 0.831470 0.707107
+
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 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.
+//
+// ---------------------------------------------------------------------
+
+
+// test that Cache::get_covering_rtree() works
+
+#include <deal.II/base/bounding_box.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/point.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/boost_adaptors/bounding_box.h>
+
+#include <deal.II/distributed/grid_refinement.h>
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+#include <deal.II/grid/manifold_lib.h>
+
+#include "../tests.h"
+
+namespace bgi = boost::geometry::index;
+
+template <int dim, int spacedim = dim>
+void
+test(unsigned int ref)
+{
+ const MPI_Comm &mpi_communicator = MPI_COMM_WORLD;
+
+ parallel::distributed::Triangulation<dim, spacedim> tria(mpi_communicator);
+ GridGenerator::hyper_ball(tria);
+ tria.refine_global(ref);
+
+ GridTools::Cache<dim, spacedim> cache(tria);
+
+ const auto &tree = cache.get_covering_rtree();
+ const auto point = random_point<dim>();
+
+ deallog << "Query point: " << point << std::endl;
+
+ for (const auto p : tree | bgi::adaptors::queried(bgi::intersects(point)))
+ deallog << "Processor " << p.second << " may own the point " << point
+ << ": it is within " << p.first.get_boundary_points().first << "; "
+ << p.first.get_boundary_points().second << std::endl;
+}
+
+int
+main(int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+ MPILogInitAll log;
+
+ test<2>(3);
+
+ return 0;
+}
--- /dev/null
+
+DEAL:0::Query point: 0.840188 0.394383
+DEAL:0::Processor 0 may own the point 0.840188 0.394383: it is within 0.316342 -0.923880; 1.00000 0.923880
--- /dev/null
+
+DEAL:0::Query point: 0.840188 0.394383
+DEAL:0::Processor 1 may own the point 0.840188 0.394383: it is within 0.316342 0.283171; 0.923880 0.923880
+DEAL:0::Processor 1 may own the point 0.840188 0.394383: it is within 0.372129 -1.11022e-16; 1.00000 0.417474
+
+DEAL:1::Query point: 0.840188 0.394383
+DEAL:1::Processor 1 may own the point 0.840188 0.394383: it is within 0.316342 0.283171; 0.923880 0.923880
+DEAL:1::Processor 1 may own the point 0.840188 0.394383: it is within 0.372129 -1.11022e-16; 1.00000 0.417474
+
#include <deal.II/base/config.h>
+#include <deal.II/base/bounding_box.h>
#include <deal.II/base/cuda.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/job_identifier.h>
+// Construct a uniformly distributed random box, with each coordinate
+// between min and max
+template <int dim>
+inline BoundingBox<dim>
+random_box(const double &min = 0.0, const double &max = 1.0)
+{
+ Assert(max >= min, ExcMessage("Make sure max>=min"));
+ std::vector<Point<dim>> p = {random_point<dim>(min, max),
+ random_point<dim>(min, max)};
+ return BoundingBox<dim>(p);
+}
+
+
+
// given the name of a file, copy it to deallog
// and then delete it
void
projects / dealii.git / commitdiff