Fix ParticleHandler::sort_particles_into_subdomains_and_cells() for lost particles

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 629c63804ac6eec74de0e3f8f884a1e3b2d6c50c..d7dd259f8d83a4b68aeb0fbfd21612e8b9a4395c 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -1067,6 +1067,11 @@ namespace Particles
                                                           cell_centers);
                     });
 
+          // make a copy of the current cell, since we will modify the variable
+          // current_cell in the following but we need a backup in the case
+          // the particle is not found
+          const auto previous_cell_of_particle = current_cell;
+
           // Search all of the cells adjacent to the closest vertex of the
           // previous cell Most likely we will find the particle in them.
           for (unsigned int i = 0; i < n_neighbor_cells; ++i)
@@ -1112,7 +1117,8 @@ namespace Particles
                   // We can find no cell for this particle. It has left the
                   // domain due to an integration error or an open boundary.
                   // Signal the loss and move on.
-                  signals.particle_lost(out_particle, current_cell);
+                  signals.particle_lost(out_particle,
+                                        previous_cell_of_particle);
                   continue;
                 }
             }

diff --git a/tests/particles/lost_01.cc b/tests/particles/lost_01.cc
new file mode 100644 (file)
index 0000000..09d17b0
--- /dev/null
+++ b/tests/particles/lost_01.cc
@@ -0,0 +1,54 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 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 the signal that is called if a particle is lost.
+
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/particles/particle_handler.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+int
+main()
+{
+  initlog();
+
+  const int dim = 2;
+
+  Triangulation<dim> tria;
+  GridGenerator::hyper_cube(tria);
+
+  MappingQ1<dim> mapping;
+
+  Particles::ParticleHandler<dim> particle_handler(tria, mapping);
+
+  std::vector<Point<dim>> particle_positions;
+  particle_positions.emplace_back(0.5, 0.5);
+  particle_handler.insert_particles(particle_positions);
+
+  particle_handler.signals.particle_lost.connect(
+    [](const auto &particle, const auto &cell) {
+      deallog << particle->get_id() << " " << cell->id() << std::endl;
+    });
+
+  particle_handler.begin()->set_location(Point<dim>{2.0, 2.0});
+
+  particle_handler.sort_particles_into_subdomains_and_cells();
+}

diff --git a/tests/particles/lost_01.output b/tests/particles/lost_01.output
new file mode 100644 (file)
index 0000000..6cc2f57
--- /dev/null
+++ b/tests/particles/lost_01.output
@@ -0,0 +1,2 @@
+
+DEAL::0 0_0: