Extract index sets from particle handlers

Extract an IndexSet with global dimensions equal to
`n_comps*particles.get_next_free_particle_index()`, containing for each
particle index, a set of `n_comps*` consecutive indices associated to
the particles that are locally owned.

The indices associated to a particle with global index `n` will be the
half open range `[n_comps*n`, n_comps*(n+1))`.

Co-authored-by: Bruno Blais <blais.bruno@gmail.com>

diff --git a/doc/news/changes/minor/20191118LucaHeltaiBrunoBlais b/doc/news/changes/minor/20191118LucaHeltaiBrunoBlais
new file mode 100644 (file)
index 0000000..9cc6b85
--- /dev/null
+++ b/doc/news/changes/minor/20191118LucaHeltaiBrunoBlais
@@ -0,0 +1,4 @@
+New: ParticleHandler::locally_relevant_ids() generates an IndexSet of the locally owned 
+particles. This can be used to construct linear algebra vectors to be used with particles.
+<br>
+(Luca Heltai, Bruno Blais, 2019/11/18)

diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index 1021fb5c99ca8691725ee904063bcd46191019d8..4ebef151b63d8f6c3abc2bd39bca247384f07e89 100644 (file)
--- a/include/deal.II/particles/particle_handler.h
+++ b/include/deal.II/particles/particle_handler.h
@@ -354,6 +354,27 @@ namespace Particles
     types::particle_index
     get_next_free_particle_index() const;
 
+    /**
+     * Extract an IndexSet with global dimensions equal to
+     * get_next_free_particle_index(), containing the locally owned
+     * particle indices.
+     *
+     * This function can be used to construct distributed vectors and matrices
+     * to manipulate particles using linear algebra operations.
+     *
+     * Notice that it is the user's responsability to guarantee that particle
+     * indices are unique, and no check is performed to verify that this is the
+     * case, nor that the union of all IndexSet objects on each mpi process is
+     * complete.
+     *
+     * @return An IndexSet of size get_next_free_particle_index(), containing
+     * n_locally_owned_particle() indices.
+     *
+     * @author Luca Heltai, Bruno Blais, 2019.
+     */
+    IndexSet
+    locally_relevant_ids() const;
+
     /**
      * Return the number of properties each particle has.
      */

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 8c6e6ab2cc4dce20f7c5f5fb919526e7870622e7..b1d725df358ceeb2753ea66ad497e6416654f0aa 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -706,6 +706,19 @@ namespace Particles
 
 
 
+  template <int dim, int spacedim>
+  IndexSet
+  ParticleHandler<dim, spacedim>::locally_relevant_ids() const
+  {
+    IndexSet set(get_next_free_particle_index());
+    for (const auto &p : *this)
+      set.add_index(p.get_id());
+    set.compress();
+    return set;
+  }
+
+
+
   template <int dim, int spacedim>
   PropertyPool &
   ParticleHandler<dim, spacedim>::get_property_pool() const

diff --git a/tests/particles/extract_index_set_01.cc b/tests/particles/extract_index_set_01.cc
new file mode 100644 (file)
index 0000000..29eab28
--- /dev/null
+++ b/tests/particles/extract_index_set_01.cc
@@ -0,0 +1,86 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 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 extract_index_set.
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/distributed/tria.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+test()
+{
+  parallel::distributed::Triangulation<dim, spacedim> tr(MPI_COMM_WORLD);
+
+  GridGenerator::hyper_cube(tr);
+  tr.refine_global(2);
+
+  MappingQ<dim, spacedim> mapping(1);
+
+  Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+  const unsigned int n_points = 10;
+  const unsigned int my_cpu = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  const unsigned int n_cpus = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+  Testing::srand(my_cpu + 1);
+
+  // Distribute the local points to the processor that owns them
+  // on the triangulation
+  auto my_bounding_box = GridTools::compute_mesh_predicate_bounding_box(
+    tr, IteratorFilters::LocallyOwnedCell());
+
+  auto global_bounding_boxes =
+    Utilities::MPI::all_gather(MPI_COMM_WORLD, my_bounding_box);
+
+  std::vector<Point<spacedim>> points(n_points);
+  for (auto &p : points)
+    p = random_point<spacedim>();
+
+  auto cpu_to_index =
+    particle_handler.insert_global_particles(points, global_bounding_boxes);
+
+  const auto set = particle_handler.locally_relevant_ids();
+  deallog << "Local set: ";
+  set.print(deallog);
+  deallog << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  MPILogInitAll all;
+
+  deallog.push("2d/2d");
+  test<2, 2>();
+  deallog.pop();
+  deallog.push("3d/3d");
+  test<3, 3>();
+  deallog.pop();
+}

diff --git a/tests/particles/extract_index_set_01.with_p4est=true.mpirun=1.output b/tests/particles/extract_index_set_01.with_p4est=true.mpirun=1.output
new file mode 100644 (file)
index 0000000..9358953
--- /dev/null
+++ b/tests/particles/extract_index_set_01.with_p4est=true.mpirun=1.output
@@ -0,0 +1,5 @@
+
+DEAL:0:2d/2d::Local set: {[0,9]}
+DEAL:0:2d/2d::
+DEAL:0:3d/3d::Local set: {[0,9]}
+DEAL:0:3d/3d::

diff --git a/tests/particles/extract_index_set_01.with_p4est=true.mpirun=2.output b/tests/particles/extract_index_set_01.with_p4est=true.mpirun=2.output
new file mode 100644 (file)
index 0000000..f595308
--- /dev/null
+++ b/tests/particles/extract_index_set_01.with_p4est=true.mpirun=2.output
@@ -0,0 +1,11 @@
+
+DEAL:0:2d/2d::Local set: {0, 2, [11,13], [15,17], 19}
+DEAL:0:2d/2d::
+DEAL:0:3d/3d::Local set: {[1,2], 6, [8,11], 13, 17, 19}
+DEAL:0:3d/3d::
+
+DEAL:1:2d/2d::Local set: {1, [3,10], 14, 18}
+DEAL:1:2d/2d::
+DEAL:1:3d/3d::Local set: {0, [3,5], 7, 12, [14,16], 18}
+DEAL:1:3d/3d::
+