From: Luca Heltai Date: Thu, 21 Nov 2019 19:04:36 +0000 (+0100) Subject: Extract index sets from particle handlers X-Git-Tag: v9.2.0-rc1~783^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ccc5b49dfdb0fd3a8ffda52a486d2f0591ca3e4d;p=dealii.git 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 --- diff --git a/doc/news/changes/minor/20191118LucaHeltaiBrunoBlais b/doc/news/changes/minor/20191118LucaHeltaiBrunoBlais new file mode 100644 index 0000000000..9cc6b854af --- /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. +
+(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 1021fb5c99..4ebef151b6 100644 --- 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 8c6e6ab2cc..b1d725df35 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -706,6 +706,19 @@ namespace Particles + template + IndexSet + ParticleHandler::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 PropertyPool & ParticleHandler::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 index 0000000000..29eab28a0e --- /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 +#include + +#include + +#include + +#include +#include + +#include + +#include "../tests.h" + +template +void +test() +{ + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + + MappingQ mapping(1); + + Particles::ParticleHandler 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> points(n_points); + for (auto &p : points) + p = random_point(); + + 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 index 0000000000..9358953f6b --- /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 index 0000000000..f595308b3c --- /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:: +