]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Get and set particle positions.
authorLuca Heltai <luca.heltai@sissa.it>
Fri, 10 Jan 2020 14:07:31 +0000 (15:07 +0100)
committerLuca Heltai <luca.heltai@sissa.it>
Sun, 19 Jan 2020 18:21:12 +0000 (19:21 +0100)
doc/news/changes/minor/20200110LucaHeltai [new file with mode: 0644]
include/deal.II/particles/particle_handler.h
source/particles/particle_handler.cc
tests/particles/set_particle_positions_01.cc [new file with mode: 0644]
tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=1.output [new file with mode: 0644]
tests/particles/set_particle_positions_01.with_p4est=true.with_trilinos=true.mpirun=2.output [new file with mode: 0644]
tests/particles/set_particle_positions_02.cc [new file with mode: 0644]
tests/particles/set_particle_positions_02.with_p4est=true.mpirun=1.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20200110LucaHeltai b/doc/news/changes/minor/20200110LucaHeltai
new file mode 100644 (file)
index 0000000..5841c01
--- /dev/null
@@ -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.
+<br>
+(Bruno Blais, Luca Heltai, 2020/01/10)
index 4ebef151b63d8f6c3abc2bd39bca247384f07e89..ef402cbd38bb7a13cde3e27bd40b90a825a83650 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <deal.II/base/array_view.h>
 #include <deal.II/base/bounding_box.h>
+#include <deal.II/base/function.h>
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/smartpointer.h>
 #include <deal.II/base/subscriptor.h>
@@ -285,6 +286,136 @@ namespace Particles
         &                                     global_bounding_boxes,
       const std::vector<std::vector<double>> &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 <class VectorType>
+    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<Point<spacedim>> &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<spacedim> &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 <class VectorType>
+    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<Point<spacedim>> &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 <int dim, int spacedim>
+  template <class VectorType>
+  void
+  ParticleHandler<dim, spacedim>::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 <int dim, int spacedim>
+  template <class VectorType>
+  void
+  ParticleHandler<dim, spacedim>::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
index b1d725df358ceeb2753ea66ad497e6416654f0aa..06ee06ba5d9bbd5b2890cbf24c8880f782ae514a 100644 (file)
@@ -719,6 +719,77 @@ namespace Particles
 
 
 
+  template <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::get_particle_positions(
+    std::vector<Point<spacedim>> &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 <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::set_particle_positions(
+    const std::vector<Point<spacedim>> &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 <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::set_particle_positions(
+    const Function<spacedim> &function,
+    const bool                displace_particles)
+  {
+    // The function should have sufficient components to displace the
+    // particles
+    AssertDimension(function.n_components, spacedim);
+
+    Vector<double> new_position(spacedim);
+    for (auto &particle : *this)
+      {
+        Point<spacedim> 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 <int dim, int spacedim>
   PropertyPool &
   ParticleHandler<dim, spacedim>::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 (file)
index 0000000..920a26b
--- /dev/null
@@ -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 <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/filtered_iterator.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/lac/trilinos_vector.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().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 (file)
index 0000000..59d6994
--- /dev/null
@@ -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 (file)
index 0000000..e425cd5
--- /dev/null
@@ -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 (file)
index 0000000..fcbbb3a
--- /dev/null
@@ -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 <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/filtered_iterator.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 = 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<Point<spacedim>> points(n_points);
+  for (auto &p : points)
+    p = random_point<spacedim>(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<Point<spacedim>> new_particle_positions(
+    particle_handler.n_locally_owned_particles());
+
+  std::vector<Point<spacedim>> 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 (file)
index 0000000..5428798
--- /dev/null
@@ -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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.