From: Luca Heltai Date: Fri, 10 Jan 2020 14:07:31 +0000 (+0100) Subject: Get and set particle positions. X-Git-Tag: v9.2.0-rc1~635^2~3 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=50eb6ea83398c8bc174513790ba82c85ac0c7008;p=dealii.git Get and set particle positions. --- diff --git a/doc/news/changes/minor/20200110LucaHeltai b/doc/news/changes/minor/20200110LucaHeltai new file mode 100644 index 0000000000..5841c01667 --- /dev/null +++ b/doc/news/changes/minor/20200110LucaHeltai @@ -0,0 +1,4 @@ +New: ParticleHandler::get_particle_positions() and ParticleHandler::set_particle_positions() allow +to set and get particle positions from various types of sources. +
+(Bruno Blais, Luca Heltai, 2020/01/10) diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index 4ebef151b6..ef402cbd38 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -20,6 +20,7 @@ #include #include +#include #include #include #include @@ -285,6 +286,136 @@ namespace Particles & global_bounding_boxes, const std::vector> &properties = {}); + /** + * Set the position of the particles by using the values contained in the + * vector @p input_vector. + * + * The vector @p input_vector should have read access to the indices + * created by extracting the locally relevant ids with + * locally_relevant_ids(), and taking its tensor + * product with the index set representing the range `[0, spacedim)`, i.e.: + * @code + * IndexSet ids = particle_handler.locally_relevant_ids(). + * tensor_product(complete_index_set(spacedim)); + * @endcode + * + * The position of the particle with global index `id` is read from + * spacedim consecutive entries starting from + * `input_vector[id*spacedim]`. + * + * If the argument @p displace_particles is set to false, then the new + * position is computed by setting it to the values contained in + * @p input_vector. By default, the particles are displaced of the + * amount contained in the @p input_vector. + * + * After setting the new position, this function calls internally the method + * sort_particles_into_subdomains_and_cells(). You should + * make sure you satisfy the requirements of that function. + * + * @param[in] input_vector A parallel distributed vector containing + * the displacement to apply to each particle, or their new absolute + * position. + * + * @param[in] displace_particles Control if the @p input_vector should + * be interpreted as a displacement vector, or a vector of absolute + * positions. + * + * @authors Luca Heltai, Bruno Blais, 2019. + */ + template + void + set_particle_positions(const VectorType &input_vector, + const bool displace_particles = true); + + /** + * Set the position of the particles within the particle handler using a + * vector of points. The new set of point defined by the + * vector has to be sufficiently close to the original one to ensure that + * the sort_particles_into_subdomains_and_cells() function manages to find + * the new cells in which the particles belong + * + * @param [in] new_positions A vector of points of dimension + * particle_handler.n_locally_owned_particles() + * + * @param [in] displace_particles When true, this add the value of the + * vector of points to the + * current position of the particle, thus displacing them by the + * amount given by the function. When false, the position of the + * particle is replaced by the value in the vector + * + * @authors Bruno Blais, Luca Heltai (2019) + */ + void + set_particle_positions(const std::vector> &new_positions, + const bool displace_particles = true); + + + /** + * Set the position of the particles within the particle handler using a + * function with spacedim components. The new set of point defined by the + * fuction has to be sufficiently close to the original one to ensure that + * the sort_particles_into_subdomains_and_cells algorithm manages to find + * the new cells in which the particles belong. + * + * @param [in] function A function that has n_components==spacedim that + * describes either the displacement or the new position of the particles + * + * @param [in] displace_particles When true, this add the results of the + * function to the current position of the particle, thus displacing them by + * the amount given by the function. When false, the position of the + * particle is is replaced by the value of the function + * + * @authors Bruno Blais, Luca Heltai (2019) + */ + void + set_particle_positions(const Function &function, + const bool displace_particles = true); + + /** + * Read the position of the particles and store them into the distributed + * vector @p output_vector. By default the + * @p output_vector is overwritten by this operation, but you can add to + * its entries by setting @p add_to_output_vector to `true`. + * + * This is the reverse operation of the set_particle_positions() function. + * The position of the particle with global index `id` is written to + * spacedim consecutive entries starting from + * `output_vector[id*spacedim]`. + * + * @param[out] output_vector A parallel distributed vector containing + * the positions of the particles, or updated with the positions of the + * particles. + * + * @param[in] add_to_output_vector Control if the function should set the + * entries of the @p output_vector or if should add to them. + * + * @author Luca Heltai, Bruno Blais, 2019. + */ + template + void + get_particle_positions(VectorType &output_vector, + const bool add_to_output_vector = false); + + /** + * Gather the position of the particles within the particle handler in + * a vector of points. + * + * @param [in,out] positions A vector preallocated at size + * (particle_handler.n_locally_owned_articles) and whose points will become + * the positions of the locally owned particles + * + * @param [in] add_to_output_vector When true, the value of the point of + * the particles is added to the positions vector. When false, + * the value of the points in the positions vector are replaced by the + * position of the particles. + * + * @authors Bruno Blais, Luca Heltai (2019) + * + */ + void + get_particle_positions(std::vector> &positions, + const bool add_to_output_vector = false); + /** * This function allows to register three additional functions that are * called every time a particle is transferred to another process @@ -619,6 +750,58 @@ namespace Particles &global_number_of_particles &global_max_particles_per_cell & next_free_particle_index; } + + + + template + template + void + ParticleHandler::set_particle_positions( + const VectorType &input_vector, + const bool displace_particles) + { + AssertDimension(input_vector.size(), + get_next_free_particle_index() * spacedim); + for (auto &p : *this) + { + const auto &point = p.get_location(); + auto new_point = point * (displace_particles ? 1.0 : 0.0); + const auto id = p.get_id(); + for (unsigned int i = 0; i < spacedim; ++i) + new_point[i] += input_vector[id * spacedim + i]; + p.set_location(new_point); + } + sort_particles_into_subdomains_and_cells(); + } + + + + template + template + void + ParticleHandler::get_particle_positions( + VectorType &output_vector, + const bool add_to_output_vector) + { + AssertDimension(output_vector.size(), + get_next_free_particle_index() * spacedim); + for (const auto &p : *this) + { + auto & point = p.get_location(); + const auto id = p.get_id(); + if (add_to_output_vector) + for (unsigned int i = 0; i < spacedim; ++i) + output_vector[id * spacedim + i] += point[i]; + else + for (unsigned int i = 0; i < spacedim; ++i) + output_vector[id * spacedim + i] = point[i]; + } + if (add_to_output_vector) + output_vector.compress(VectorOperation::add); + else + output_vector.compress(VectorOperation::insert); + } + } // namespace Particles DEAL_II_NAMESPACE_CLOSE diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index b1d725df35..06ee06ba5d 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -719,6 +719,77 @@ namespace Particles + template + void + ParticleHandler::get_particle_positions( + std::vector> &positions, + const bool add_to_output_vector) + { + // There should be one point per particle to gather + AssertDimension(positions.size(), n_locally_owned_particles()); + + unsigned int i = 0; + for (auto it = begin(); it != end(); ++it, ++i) + { + if (add_to_output_vector) + positions[i] = positions[i] + it->get_location(); + else + positions[i] = it->get_location(); + } + } + + + + template + void + ParticleHandler::set_particle_positions( + const std::vector> &new_positions, + const bool displace_particles) + { + // There should be one point per particle to fix the new position + AssertDimension(new_positions.size(), n_locally_owned_particles()); + + unsigned int i = 0; + for (auto it = begin(); it != end(); ++it, ++i) + { + if (displace_particles) + it->set_location(it->get_location() + new_positions[i]); + else + it->set_location(new_positions[i]); + } + sort_particles_into_subdomains_and_cells(); + } + + + + template + void + ParticleHandler::set_particle_positions( + const Function &function, + const bool displace_particles) + { + // The function should have sufficient components to displace the + // particles + AssertDimension(function.n_components, spacedim); + + Vector new_position(spacedim); + for (auto &particle : *this) + { + Point particle_location = particle.get_location(); + function.vector_value(particle_location, new_position); + if (displace_particles) + for (unsigned int d = 0; d < spacedim; ++d) + particle_location[d] = particle_location[d] + new_position[d]; + else + for (unsigned int d = 0; d < spacedim; ++d) + particle_location[d] = new_position[d]; + particle.set_location(particle_location); + } + sort_particles_into_subdomains_and_cells(); + } + + + template PropertyPool & ParticleHandler::get_property_pool() const diff --git a/tests/particles/set_particle_positions_01.cc b/tests/particles/set_particle_positions_01.cc new file mode 100644 index 0000000000..920a26bc9c --- /dev/null +++ b/tests/particles/set_particle_positions_01.cc @@ -0,0 +1,108 @@ +// --------------------------------------------------------------------- +// +// 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 get and set particle positions + +#include +#include + +#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().tensor_product( + complete_index_set(spacedim)); + + TrilinosWrappers::MPI::Vector vector(set, MPI_COMM_WORLD); + particle_handler.get_particle_positions(vector); + deallog << "Local set: "; + set.print(deallog); + deallog << std::endl + << "Local vector norm: " << vector.linfty_norm() << std::endl; + + // With this we know we don't get out of the domain + vector *= 1e-3; + particle_handler.set_particle_positions(vector); + + // Make sure we have a new index set + + const auto new_set = particle_handler.locally_relevant_ids().tensor_product( + complete_index_set(spacedim)); + + TrilinosWrappers::MPI::Vector new_vector(new_set, MPI_COMM_WORLD); + + particle_handler.get_particle_positions(new_vector); + deallog << "New position norm: " << vector.linfty_norm() << 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/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=1.output b/tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=1.output new file mode 100644 index 0000000000..59d6994e96 --- /dev/null +++ b/tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=1.output @@ -0,0 +1,9 @@ + +DEAL:0:2d/2d::Local set: {[0,19]} +DEAL:0:2d/2d:: +DEAL:0:2d/2d::Local vector norm: 0.952230 +DEAL:0:2d/2d::New position norm: 0.000952230 +DEAL:0:3d/3d::Local set: {[0,29]} +DEAL:0:3d/3d:: +DEAL:0:3d/3d::Local vector norm: 0.998925 +DEAL:0:3d/3d::New position norm: 0.000998925 diff --git a/tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=2.output b/tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=2.output new file mode 100644 index 0000000000..e425cd53f7 --- /dev/null +++ b/tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=2.output @@ -0,0 +1,19 @@ + +DEAL:0:2d/2d::Local set: {[0,1], [4,5], [22,27], [30,35], [38,39]} +DEAL:0:2d/2d:: +DEAL:0:2d/2d::Local vector norm: 0.990603 +DEAL:0:2d/2d::New position norm: 0.000990603 +DEAL:0:3d/3d::Local set: {[3,8], [18,20], [24,35], [39,41], [51,53], [57,59]} +DEAL:0:3d/3d:: +DEAL:0:3d/3d::Local vector norm: 0.998925 +DEAL:0:3d/3d::New position norm: 0.000998925 + +DEAL:1:2d/2d::Local set: {[2,3], [6,21], [28,29], [36,37]} +DEAL:1:2d/2d:: +DEAL:1:2d/2d::Local vector norm: 0.990603 +DEAL:1:2d/2d::New position norm: 0.000990603 +DEAL:1:3d/3d::Local set: {[0,2], [9,17], [21,23], [36,38], [42,50], [54,56]} +DEAL:1:3d/3d:: +DEAL:1:3d/3d::Local vector norm: 0.998925 +DEAL:1:3d/3d::New position norm: 0.000998925 + diff --git a/tests/particles/set_particle_positions_02.cc b/tests/particles/set_particle_positions_02.cc new file mode 100644 index 0000000000..fcbbb3a209 --- /dev/null +++ b/tests/particles/set_particle_positions_02.cc @@ -0,0 +1,110 @@ +// --------------------------------------------------------------------- +// +// 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 insert_global_particles. Make sure we don't lose particles +// along the way + +#include +#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 = 2; + 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); + + // Create the bounding boxes for the global insertion + 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(0.1, 0.2); + + auto cpu_to_index = + particle_handler.insert_global_particles(points, global_bounding_boxes); + + for (const auto &particle : particle_handler) + { + deallog << "Particle id : " << particle.get_id(); + deallog << " Particle location: " << particle.get_location() << std::endl; + } + + Tensor<1, spacedim> displacement; + displacement[0] = 0.1; + std::vector> new_particle_positions( + particle_handler.n_locally_owned_particles()); + + std::vector> old_particle_positions( + particle_handler.n_locally_owned_particles()); + particle_handler.get_particle_positions(old_particle_positions); + + for (unsigned int i = 0; i < old_particle_positions.size(); ++i) + new_particle_positions[i] = old_particle_positions[i] + displacement; + + particle_handler.set_particle_positions(new_particle_positions, false); + + deallog << "Displacement" << std::endl; + for (const auto &particle : particle_handler) + { + deallog << "Particle id : " << particle.get_id(); + deallog << " Particle location: " << particle.get_location() << 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 << "---" << std::endl; + deallog.push("3d/3d"); + test<3, 3>(); + deallog.pop(); +} diff --git a/tests/particles/set_particle_positions_02.with_p4est=true.mpirun=1.output b/tests/particles/set_particle_positions_02.with_p4est=true.mpirun=1.output new file mode 100644 index 0000000000..54287987a4 --- /dev/null +++ b/tests/particles/set_particle_positions_02.with_p4est=true.mpirun=1.output @@ -0,0 +1,12 @@ + +DEAL:0:2d/2d::Particle id : 0 Particle location: 0.184019 0.139438 +DEAL:0:2d/2d::Particle id : 1 Particle location: 0.178310 0.179844 +DEAL:0:2d/2d::Displacement +DEAL:0:2d/2d::Particle id : 0 Particle location: 0.284019 0.139438 +DEAL:0:2d/2d::Particle id : 1 Particle location: 0.278310 0.179844 +DEAL:0::--- +DEAL:0:3d/3d::Particle id : 0 Particle location: 0.184019 0.139438 0.178310 +DEAL:0:3d/3d::Particle id : 1 Particle location: 0.179844 0.191165 0.119755 +DEAL:0:3d/3d::Displacement +DEAL:0:3d/3d::Particle id : 0 Particle location: 0.284019 0.139438 0.178310 +DEAL:0:3d/3d::Particle id : 1 Particle location: 0.279844 0.191165 0.119755