Fix bug in remove particle

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 81fcefaf007416690a6a4386489046ed46b838f1..2269b73c74bba32315b5b7f2d3e173a96ef986f1 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -357,7 +357,7 @@ namespace Particles
   ParticleHandler<dim, spacedim>::remove_particle(
     const ParticleHandler<dim, spacedim>::particle_iterator &particle)
   {
-    auto &container = *(particle->particles_on_cell);
+    auto &particles_on_cell = *(particle->particles_on_cell);
 
     // if the particle has an invalid handle (e.g. because it has
     // been duplicated before calling this function) do not try
@@ -368,15 +368,15 @@ namespace Particles
         property_pool->deregister_particle(handle);
       }
 
-    if (container.size() > 1)
+    if (particles_on_cell.size() > 1)
       {
-        container[particle->particle_index_within_cell] =
-          std::move(container.back());
-        container.resize(particles.size() - 1);
+        particles_on_cell[particle->particle_index_within_cell] =
+          std::move(particles_on_cell.back());
+        particles_on_cell.resize(particles_on_cell.size() - 1);
       }
     else
       {
-        container.clear();
+        particles_on_cell.clear();
       }
 
     --local_number_of_particles;

diff --git a/tests/particles/remove_particle_01.cc b/tests/particles/remove_particle_01.cc
new file mode 100644 (file)
index 0000000..edf8334
--- /dev/null
+++ b/tests/particles/remove_particle_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.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the removal of a single particle works as expected
+
+#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/generators.h>
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+template <int dim, int spacedim>
+void
+test()
+{
+  {
+    Triangulation<dim, spacedim> tr;
+
+    GridGenerator::hyper_cube(tr);
+    tr.refine_global(1);
+    MappingQ<dim, spacedim> mapping(1);
+
+    Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
+
+    std::vector<Point<dim>> particle_reference_locations(3, Point<dim>());
+
+    for (unsigned int i = 0; i < dim; ++i)
+      {
+        particle_reference_locations[0](i) = 0.25;
+        particle_reference_locations[1](i) = 0.5;
+        particle_reference_locations[2](i) = 0.75;
+      }
+
+    Particles::Generators::regular_reference_locations(
+      tr, particle_reference_locations, particle_handler);
+
+    deallog << "Particle number: " << particle_handler.n_global_particles()
+            << std::endl;
+
+    for (const auto &cell : tr.active_cell_iterators())
+      {
+        const unsigned int n_particles_to_remove =
+          cell->active_cell_index() % 4;
+
+        for (unsigned int i = 0; i < n_particles_to_remove; ++i)
+          {
+            const unsigned int particle_index_to_remove =
+              Testing::rand() % particle_handler.n_particles_in_cell(cell);
+            auto particle_to_remove =
+              particle_handler.particles_in_cell(cell).begin();
+            std::advance(particle_to_remove, particle_index_to_remove);
+
+            deallog << "Removing particle index: "
+                    << particle_to_remove->get_id()
+                    << ". Advanced by: " << particle_index_to_remove
+                    << std::endl;
+
+            particle_handler.remove_particle(particle_to_remove);
+          }
+
+        deallog << "Cell index: " << cell->active_cell_index()
+                << ". N Particles: "
+                << particle_handler.n_particles_in_cell(cell) << std::endl;
+      }
+
+    for (const auto &particle : particle_handler)
+      {
+        deallog << "Particle id: " << particle.get_id() << std::endl;
+      }
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main(int argc, char *argv[])
+{
+  initlog();
+
+  deallog.push("2d/2d");
+  test<2, 2>();
+  deallog.pop();
+  deallog.push("3d/3d");
+  test<3, 3>();
+  deallog.pop();
+}

diff --git a/tests/particles/remove_particle_01.output b/tests/particles/remove_particle_01.output
new file mode 100644 (file)
index 0000000..af546c9
--- /dev/null
+++ b/tests/particles/remove_particle_01.output
@@ -0,0 +1,53 @@
+
+DEAL:2d/2d::Particle number: 12
+DEAL:2d/2d::Cell index: 0. N Particles: 3
+DEAL:2d/2d::Removing particle index: 4. Advanced by: 1
+DEAL:2d/2d::Cell index: 1. N Particles: 2
+DEAL:2d/2d::Removing particle index: 7. Advanced by: 1
+DEAL:2d/2d::Removing particle index: 8. Advanced by: 1
+DEAL:2d/2d::Cell index: 2. N Particles: 1
+DEAL:2d/2d::Removing particle index: 10. Advanced by: 1
+DEAL:2d/2d::Removing particle index: 11. Advanced by: 1
+DEAL:2d/2d::Removing particle index: 9. Advanced by: 0
+DEAL:2d/2d::Cell index: 3. N Particles: 0
+DEAL:2d/2d::Particle id: 0
+DEAL:2d/2d::Particle id: 1
+DEAL:2d/2d::Particle id: 2
+DEAL:2d/2d::Particle id: 3
+DEAL:2d/2d::Particle id: 5
+DEAL:2d/2d::Particle id: 6
+DEAL:2d/2d::OK
+DEAL:3d/3d::Particle number: 24
+DEAL:3d/3d::Cell index: 0. N Particles: 3
+DEAL:3d/3d::Removing particle index: 4. Advanced by: 1
+DEAL:3d/3d::Cell index: 1. N Particles: 2
+DEAL:3d/3d::Removing particle index: 6. Advanced by: 0
+DEAL:3d/3d::Removing particle index: 7. Advanced by: 1
+DEAL:3d/3d::Cell index: 2. N Particles: 1
+DEAL:3d/3d::Removing particle index: 10. Advanced by: 1
+DEAL:3d/3d::Removing particle index: 9. Advanced by: 0
+DEAL:3d/3d::Removing particle index: 11. Advanced by: 0
+DEAL:3d/3d::Cell index: 3. N Particles: 0
+DEAL:3d/3d::Cell index: 4. N Particles: 3
+DEAL:3d/3d::Removing particle index: 17. Advanced by: 2
+DEAL:3d/3d::Cell index: 5. N Particles: 2
+DEAL:3d/3d::Removing particle index: 19. Advanced by: 1
+DEAL:3d/3d::Removing particle index: 20. Advanced by: 1
+DEAL:3d/3d::Cell index: 6. N Particles: 1
+DEAL:3d/3d::Removing particle index: 22. Advanced by: 1
+DEAL:3d/3d::Removing particle index: 21. Advanced by: 0
+DEAL:3d/3d::Removing particle index: 23. Advanced by: 0
+DEAL:3d/3d::Cell index: 7. N Particles: 0
+DEAL:3d/3d::Particle id: 0
+DEAL:3d/3d::Particle id: 1
+DEAL:3d/3d::Particle id: 2
+DEAL:3d/3d::Particle id: 3
+DEAL:3d/3d::Particle id: 5
+DEAL:3d/3d::Particle id: 8
+DEAL:3d/3d::Particle id: 12
+DEAL:3d/3d::Particle id: 13
+DEAL:3d/3d::Particle id: 14
+DEAL:3d/3d::Particle id: 15
+DEAL:3d/3d::Particle id: 16
+DEAL:3d/3d::Particle id: 18
+DEAL:3d/3d::OK