Keep both owned/ghost particles in a single container

This is needed to keep _GLIBCXX_DEBUG happy when comparing iterators from cache.
This commit also addresses the review comments

diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h
index 4b3ff0910cd0bc580749a38116f44fe7dd56e771..d5a3b3b3aedce219c04629c70e5f88cd1c1c7348 100644 (file)
--- a/include/deal.II/particles/particle_accessor.h
+++ b/include/deal.II/particles/particle_accessor.h
@@ -46,8 +46,37 @@ namespace Particles
   {
   public:
     /**
-     * Data structure to describes the particles in a given cell. This is used
-     * inside an std::list in `particle_container`.
+     * Data structure to describe the particles in a given cell. This is used
+     * inside an std::list in `particle_container`. The storage of this field
+     * is typically handled by ParticleHandler, using an std::list of this
+     * structure.
+     *
+     * There are four main reasons for the present design:
+     * <ul>
+     * <li>The list represents a compact storage for all particles on a cell,
+     * handling e.g. the case where only a few cells hold particles without
+     * storage overhead.</li>
+     * <li>The particles on a cell are identified by their handle object,
+     * i.e., a single integer, again ensuring only 4 byte storage cost for the
+     * handling of the particles for the case of sufficiently many particles
+     * per cell.</li>
+     * <li>Combined with a cache object that holds an iterator to the
+     * `std::list` of particles (8 bytes per cell), or the `std::list::end` in
+     * case no particles are present on a cell, there is a fast access from a
+     * ParticleAccessor -> surrounding cell (by access to
+     * ParticlesInCell::cell_iterator), as well as fast access from cell ->
+     * all particles (through the cache). It also allows for a fast iteration
+     * through all particles, by incrementing either the index of particles
+     * within a cell, or, if at the end of the cell, to the next element in
+     * the outer list. The cache is simple to keep consistent because the
+     * iterators into `std::list` remain valid upon insertion or deletion of
+     * entries in the list, as specified by `std::list`'s API.
+     * <li>In order to detect the start and end of all cells with particles,
+     * the std::list contains some dummy entry with a `cell_iterator` past the
+     * end of valid cells, which is used as a criterion to terminate the loops
+     * of ParticleAccessor, again minimizing the computational cost of handling
+     * the loop over particles.</li>
+     * </ul>
      */
     struct ParticlesInCell
     {
@@ -62,10 +91,9 @@ namespace Particles
       ParticlesInCell(
         const std::vector<typename PropertyPool<dim, spacedim>::Handle>
           &particles,
-        const typename Triangulation<dim, spacedim>::active_cell_iterator
-          &cell_iterator)
+        const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
         : particles(particles)
-        , cell_iterator(cell_iterator)
+        , cell(cell)
       {}
 
       /**
@@ -76,7 +104,7 @@ namespace Particles
       /**
        * The underlying cell.
        */
-      typename Triangulation<dim, spacedim>::active_cell_iterator cell_iterator;
+      typename Triangulation<dim, spacedim>::active_cell_iterator cell;
     };
 
     /**
@@ -391,7 +419,6 @@ namespace Particles
      * friend classes.
      */
     ParticleAccessor(
-      const particle_container &                  particles,
       const typename particle_container::iterator particles_in_cell,
       const PropertyPool<dim, spacedim> &         property_pool,
       const unsigned int                          particle_index_within_cell);
@@ -411,13 +438,8 @@ namespace Particles
     get_handle();
 
     /**
-     * A pointer to the container that stores the particles.
-     */
-    const particle_container *particles;
-
-    /**
-     * An iterator to the currently active particles within the particles
-     * object.
+     * An iterator to the particles in the current cell within the
+     * particle_container object.
      */
     typename particle_container::iterator particles_in_cell;
 
@@ -501,8 +523,7 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor()
-    : particles(nullptr)
-    , particles_in_cell(typename particle_container::iterator())
+    : particles_in_cell(typename particle_container::iterator())
     , property_pool(nullptr)
     , particle_index_within_cell(numbers::invalid_unsigned_int)
   {}
@@ -511,12 +532,10 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor(
-    const particle_container &                  particles,
     const typename particle_container::iterator particles_in_cell,
     const PropertyPool<dim, spacedim> &         property_pool,
     const unsigned int                          particle_index_within_cell)
-    : particles(&particles)
-    , particles_in_cell(particles_in_cell)
+    : particles_in_cell(particles_in_cell)
     , property_pool(const_cast<PropertyPool<dim, spacedim> *>(&property_pool))
     , particle_index_within_cell(particle_index_within_cell)
   {}
@@ -776,10 +795,10 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::get_surrounding_cell() const
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
-    Assert(particles_in_cell->cell_iterator.state() == IteratorState::valid,
+    Assert(particles_in_cell->cell.state() == IteratorState::valid,
            ExcInternalError());
 
-    return particles_in_cell->cell_iterator;
+    return particles_in_cell->cell;
   }
 
 
@@ -790,10 +809,10 @@ namespace Particles
     const Triangulation<dim, spacedim> & /*triangulation*/) const
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
-    Assert(particles_in_cell->cell_iterator.state() == IteratorState::valid,
+    Assert(particles_in_cell->cell.state() == IteratorState::valid,
            ExcInternalError());
 
-    return particles_in_cell->cell_iterator;
+    return particles_in_cell->cell;
   }
 
 
@@ -854,17 +873,10 @@ namespace Particles
       --particle_index_within_cell;
     else
       {
-        if (particles_in_cell == particles->begin())
-          {
-            particles_in_cell =
-              const_cast<particle_container *>(particles)->end();
-          }
-        else
-          {
-            --particles_in_cell;
-            particle_index_within_cell =
-              particles_in_cell->particles.size() - 1;
-          }
+        --particles_in_cell;
+        particle_index_within_cell = particles_in_cell->particles.empty() ?
+                                       0 :
+                                       particles_in_cell->particles.size() - 1;
       }
   }
 
@@ -896,11 +908,12 @@ namespace Particles
   inline IteratorState::IteratorStates
   ParticleAccessor<dim, spacedim>::state() const
   {
-    if (particles != nullptr && property_pool != nullptr &&
-        particles_in_cell != particles->end() &&
+    if (property_pool != nullptr &&
+        particles_in_cell->cell.state() == IteratorState::valid &&
         particle_index_within_cell < particles_in_cell->particles.size())
       return IteratorState::valid;
-    else if (particles != nullptr && particles_in_cell == particles->end() &&
+    else if (property_pool != nullptr &&
+             particles_in_cell->cell.state() == IteratorState::past_the_end &&
              particle_index_within_cell == 0)
       return IteratorState::past_the_end;
     else

diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index ddc28596211f4dfff53963a766f145c7eb3e29c2..1b6e5b5f45e79640f17f9cc0376fd26c5800d834 100644 (file)
--- a/include/deal.II/particles/particle_handler.h
+++ b/include/deal.II/particles/particle_handler.h
@@ -880,13 +880,25 @@ namespace Particles
 
     /**
      * Perform the local insertion operation into the particle container. This
-     * function is used during the higher-level functions inserting particles.
+     * function is used in the higher-level functions inserting particles.
      */
-    void
+    particle_iterator
     insert_particle(
       const typename PropertyPool<dim, spacedim>::Handle          handle,
-      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-      particle_container &                                        particles);
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell);
+
+    /**
+     * Delete all entries in the particles container, and then set the three
+     * anchor entries used to distinguish between owned and ghost cells: We
+     * add one item to the front, one between and one as the last element in
+     * the particle container to be able to iterate across particles without
+     * `if` statements, solely relying on ParticleAccessor::operator== to
+     * terminate operations, and using the `cell_iterator` inside the particle
+     * container to check for valid states.
+     */
+    void
+    reset_particle_container(const Triangulation<dim, spacedim> *trianulation,
+                             particle_container &                particles);
 
     /**
      * Address of the triangulation to work on.
@@ -912,14 +924,15 @@ namespace Particles
     std::unique_ptr<PropertyPool<dim, spacedim>> property_pool;
 
     /**
-     * Set of particles currently living in the locally owned cells.
+     * Set of particles currently living in the locally owned or ghost cells.
      */
-    particle_container owned_particles;
+    particle_container particles;
 
     /**
-     * Set of particles currently living in the ghost cells.
+     * Iterator to the start of the end of the list elements of
+     * particle_container which belong to locally owned elements.
      */
-    particle_container ghost_particles;
+    typename particle_container::iterator owned_particles_end;
 
     /**
      * List from the active cells on the present MPI process to positions in
@@ -1018,15 +1031,11 @@ namespace Particles
      * Transfer particles that have crossed subdomain boundaries to other
      * processors.
      * All received particles and their new cells will be appended to the
-     * @p received_particles vector.
+     * class variable `particles` at the right slot.
      *
      * @param [in] particles_to_send All particles that should be sent and
      * their new subdomain_ids are in this map.
      *
-     * @param [in,out] received_particles Particle container that stores all received
-     * particles. Note that it is not required nor checked that the container
-     * is empty, received particles are simply inserted into the container.
-     *
      * @param [in] new_cells_for_particles Optional vector of cell
      * iterators with the same structure as @p particles_to_send. If this
      * parameter is given it should contain the cell iterator for every
@@ -1044,8 +1053,7 @@ namespace Particles
     void
     send_recv_particles(
       const std::map<types::subdomain_id, std::vector<particle_iterator>>
-        &                 particles_to_send,
-      particle_container &received_particles,
+        &particles_to_send,
       const std::map<
         types::subdomain_id,
         std::vector<
@@ -1057,27 +1065,20 @@ namespace Particles
       const bool enable_cache = false);
 
     /**
-     * Transfer ghost particles' position and properties assuming that
-     * the particles have not changed cells. This routine uses the
+     * Transfer ghost particles' position and properties assuming that the
+     * particles have not changed cells. This routine uses the
      * GhostParticlePartitioner as a caching structure to know which particles
-     * are ghost to other processes, and where they need to be sent.
-     * It inherently assumes that particles cannot have changed cell.
+     * are ghost to other processes, and where they need to be sent.  It
+     * inherently assumes that particles cannot have changed cell, and writes
+     * the result back to the `particles` member variable.
      *
      * @param [in] particles_to_send All particles for which information
      * should be sent and their new subdomain_ids are in this map.
-     *
-     * @param [in,out] received_particles A map with all received
-     * particles. Note that it is not required nor checked that the container
-     * is empty, received particles are simply inserted into
-     * the container.
-     *
      */
     void
     send_recv_particles_properties_and_location(
       const std::map<types::subdomain_id, std::vector<particle_iterator>>
-        &                 particles_to_send,
-      particle_container &received_particles);
-
+        &particles_to_send);
 
 #endif
 
@@ -1171,13 +1172,10 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::begin()
   {
-    if (owned_particles.empty())
-      return end();
-    else
-      return particle_iterator(owned_particles,
-                               owned_particles.begin(),
-                               *property_pool,
-                               0);
+    Assert(!particles.empty(), ExcInternalError());
+
+    // jump over the first entry which is used for book-keeping
+    return particle_iterator(++particles.begin(), *property_pool, 0);
   }
 
 
@@ -1195,10 +1193,7 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end()
   {
-    return particle_iterator(owned_particles,
-                             owned_particles.end(),
-                             *property_pool,
-                             0);
+    return particle_iterator(owned_particles_end, *property_pool, 0);
   }
 
 
@@ -1216,13 +1211,13 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::begin_ghost()
   {
-    if (ghost_particles.empty())
-      return end_ghost();
-    else
-      return particle_iterator(ghost_particles,
-                               ghost_particles.begin(),
-                               *property_pool,
-                               0);
+    Assert(!particles.empty(), ExcInternalError());
+    // find the start of ghost particles as one past the end of owned indices;
+    // the index in between is used for book-keeping
+    typename particle_container::iterator ghost_particles_begin =
+      owned_particles_end;
+    ++ghost_particles_begin;
+    return particle_iterator(ghost_particles_begin, *property_pool, 0);
   }
 
 
@@ -1240,10 +1235,8 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end_ghost()
   {
-    return particle_iterator(ghost_particles,
-                             ghost_particles.end(),
-                             *property_pool,
-                             0);
+    // skip the last entry which is used for book-keeping
+    return particle_iterator(--particles.end(), *property_pool, 0);
   }
 
 

diff --git a/include/deal.II/particles/particle_iterator.h b/include/deal.II/particles/particle_iterator.h
index ea89a4368ab919d1b13d8844a2355158795d792f..75fe59acb38d99c17a3f646c1eff10da319cc974 100644 (file)
--- a/include/deal.II/particles/particle_iterator.h
+++ b/include/deal.II/particles/particle_iterator.h
@@ -56,7 +56,6 @@ namespace Particles
      * cell.
      */
     ParticleIterator(
-      const particle_container &                  particles,
       const typename particle_container::iterator particles_in_cell,
       const PropertyPool<dim, spacedim> &         property_pool,
       const unsigned int                          particle_index_within_cell);
@@ -163,14 +162,10 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleIterator<dim, spacedim>::ParticleIterator(
-    const particle_container &                  particles,
     const typename particle_container::iterator particles_in_cell,
     const PropertyPool<dim, spacedim> &         property_pool,
     const unsigned int                          particle_index_within_cell)
-    : accessor(particles,
-               particles_in_cell,
-               property_pool,
-               particle_index_within_cell)
+    : accessor(particles_in_cell, property_pool, particle_index_within_cell)
   {}
 
 

diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index a75d95a98130a488e4eea5617f6ffb3c7090a81f..c00f560faf3f9d9051c595a3776bacfe899610b1 100644 (file)
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -51,6 +51,8 @@ namespace Particles
     }
   } // namespace
 
+
+
   template <int dim, int spacedim>
   ParticleHandler<dim, spacedim>::ParticleHandler()
     : triangulation()
@@ -65,7 +67,9 @@ namespace Particles
     , load_callback()
     , handle(numbers::invalid_unsigned_int)
     , tria_listeners()
-  {}
+  {
+    reset_particle_container(nullptr, particles);
+  }
 
 
 
@@ -77,7 +81,7 @@ namespace Particles
     : triangulation(&triangulation, typeid(*this).name())
     , mapping(&mapping, typeid(*this).name())
     , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
-    , cells_to_particle_cache(triangulation.n_active_cells())
+    , cells_to_particle_cache(triangulation.n_active_cells(), particles.end())
     , global_number_of_particles(0)
     , number_of_locally_owned_particles(0)
     , global_max_particles_per_cell(0)
@@ -91,6 +95,7 @@ namespace Particles
                                                           mapping))
     , tria_listeners()
   {
+    reset_particle_container(&triangulation, particles);
     connect_to_triangulation_signals();
   }
 
@@ -119,6 +124,8 @@ namespace Particles
     triangulation = &new_triangulation;
     mapping       = &new_mapping;
 
+    reset_particle_container(&new_triangulation, particles);
+
     // Create the memory pool that will store all particle properties
     property_pool = std::make_unique<PropertyPool<dim, spacedim>>(n_properties);
 
@@ -128,7 +135,8 @@ namespace Particles
       std::make_unique<GridTools::Cache<dim, spacedim>>(new_triangulation,
                                                         new_mapping);
 
-    cells_to_particle_cache.resize(triangulation->n_active_cells());
+    cells_to_particle_cache.resize(triangulation->n_active_cells(),
+                                   particles.end());
 
     connect_to_triangulation_signals();
   }
@@ -140,9 +148,6 @@ namespace Particles
   ParticleHandler<dim, spacedim>::copy_from(
     const ParticleHandler<dim, spacedim> &particle_handler)
   {
-    // clear this object before copying particles
-    clear();
-
     const unsigned int n_properties =
       particle_handler.property_pool->n_properties_per_slot();
     initialize(*particle_handler.triangulation,
@@ -160,9 +165,16 @@ namespace Particles
     global_max_particles_per_cell =
       particle_handler.global_max_particles_per_cell;
     next_free_particle_index = particle_handler.next_free_particle_index;
-    cells_to_particle_cache  = particle_handler.cells_to_particle_cache;
-    owned_particles          = particle_handler.owned_particles;
-    ghost_particles          = particle_handler.ghost_particles;
+    particles                = particle_handler.particles;
+    owned_particles_end      = particles.begin();
+    for (auto it = particle_handler.particles.begin();
+         it != particle_handler.owned_particles_end;
+         ++it)
+      ++owned_particles_end;
+
+    for (auto it = particles.begin(); it != particles.end(); ++it)
+      if (!it->particles.empty())
+        cells_to_particle_cache[it->cell->active_cell_index()] = it;
 
     ghost_particles_cache.ghost_particles_by_domain =
       particle_handler.ghost_particles_cache.ghost_particles_by_domain;
@@ -188,21 +200,20 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::clear_particles()
   {
-    for (auto &particles_in_cell : owned_particles)
-      for (auto &particle : particles_in_cell.particles)
-        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
-          property_pool->deregister_particle(particle);
-
-    for (auto &particles_in_cell : ghost_particles)
+    for (auto &particles_in_cell : particles)
       for (auto &particle : particles_in_cell.particles)
         if (particle != PropertyPool<dim, spacedim>::invalid_handle)
           property_pool->deregister_particle(particle);
 
     cells_to_particle_cache.clear();
-    owned_particles.clear();
-    ghost_particles.clear();
     if (triangulation != nullptr)
-      cells_to_particle_cache.resize(triangulation->n_active_cells());
+      {
+        cells_to_particle_cache.resize(triangulation->n_active_cells(),
+                                       particles.end());
+        reset_particle_container(&*triangulation, particles);
+      }
+    else
+      reset_particle_container(nullptr, particles);
 
     // the particle properties have already been deleted by their destructor,
     // but the memory is still allocated. Return the memory as well.
@@ -211,69 +222,144 @@ namespace Particles
 
 
 
-  namespace
+  template <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::reset_particle_container(
+    const Triangulation<dim, spacedim> *triangulation,
+    particle_container &                given_particles)
   {
-    template <int dim, int spacedim, typename ParticleContainer>
-    std::array<types::particle_index, 3>
-    my_update_cached_numbers(const PropertyPool<dim, spacedim> &property_pool,
-                             ParticleContainer &                particles)
-    {
-      // sort the particles by the active cell index
-      particles.sort([](const typename ParticleContainer::value_type &a,
-                        const typename ParticleContainer::value_type &b) {
-        return a.cell_iterator < b.cell_iterator;
-      });
-
-      std::array<types::particle_index, 3> result = {};
-      for (const auto &particles_in_cell : particles)
-        {
-          const types::particle_index n_particles_in_cell =
-            particles_in_cell.particles.size();
-
-          // local_max_particles_per_cell
-          result[0] = std::max(result[0], n_particles_in_cell);
+    TriaActiveIterator<CellAccessor<dim, spacedim>> past_the_end_iterator(
+      triangulation, -1, -1);
+    given_particles.clear();
+    for (unsigned int i = 0; i < 3; ++i)
+      given_particles.emplace_back(
+        std::vector<typename PropertyPool<dim, spacedim>::Handle>(),
+        past_the_end_iterator);
+
+    // Set the end of owned particles to the middle of the three elements
+    owned_particles_end = ++given_particles.begin();
+  }
 
-          // number of locally owned particles
-          result[2] += n_particles_in_cell;
 
-          // local_max_particle_index
-          for (const auto &particle : particles_in_cell.particles)
+  template <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::update_cached_numbers()
+  {
+    // first sort the owned particles by the active cell index
+    bool sort_is_necessary = false;
+    {
+      auto previous = ++particles.begin();
+      for (auto next = previous; next != owned_particles_end; ++next)
+        {
+          if (previous->cell.state() == IteratorState::valid &&
+              next->cell.state() == IteratorState::valid &&
+              previous->cell > next->cell)
             {
-              result[1] = std::max(result[1], property_pool.get_id(particle));
+              sort_is_necessary = true;
+              break;
             }
+          previous = next;
         }
-
-      return result;
     }
-  } // namespace
-
+    if (sort_is_necessary)
+      {
+        // we could have tried to call std::list::sort with a custom
+        // comparator to get things sorted, but things are complicated by the
+        // three anchor entries that we do not want to move, and we would
+        // hence pay an O(N log(N)) algorithm with a large constant in front
+        // of it (on the order of 20+ instructions with many
+        // difficult-to-predict branches). Therefore, we simply copy the list
+        // into a new one (keeping alive the possibly large vectors with
+        // particles on cells) into a new container.
+        particle_container sorted_particles;
+        reset_particle_container(&*triangulation, sorted_particles);
+
+        // iterate over cells and insert the entries in the new order
+        for (const auto &cell : triangulation->active_cell_iterators())
+          if (!cell->is_artificial())
+            if (cells_to_particle_cache[cell->active_cell_index()] !=
+                particles.end())
+              {
+                // observe that owned_particles_end was already set to point
+                // into the new sorted_particles array
+                typename particle_container::iterator insert_position =
+                  cell->is_locally_owned() ? owned_particles_end :
+                                             --sorted_particles.end();
+                typename particle_container::iterator new_entry =
+                  sorted_particles.insert(
+                    insert_position, typename particle_container::value_type());
+                new_entry->cell = cell;
+                new_entry->particles =
+                  std::move(cells_to_particle_cache[cell->active_cell_index()]
+                              ->particles);
+              }
+        particles = std::move(sorted_particles);
+
+        // refresh cells_to_particle_cache
+        cells_to_particle_cache.clear();
+        cells_to_particle_cache.resize(triangulation->n_active_cells(),
+                                       particles.end());
+        for (auto it = particles.begin(); it != particles.end(); ++it)
+          if (!it->particles.empty())
+            cells_to_particle_cache[it->cell->active_cell_index()] = it;
+      }
 
-  template <int dim, int spacedim>
-  void
-  ParticleHandler<dim, spacedim>::update_cached_numbers()
-  {
-    // first compute local result with the function above and then compute the
-    // collective results
-    const std::array<types::particle_index, 3> result_ghosted =
-      my_update_cached_numbers(*property_pool, ghost_particles);
-    std::array<types::particle_index, 3> result_owned =
-      my_update_cached_numbers(*property_pool, owned_particles);
+    // Ensure that we did not accidentally modify the anchor entries with
+    // special purpose.
+    Assert(particles.front().cell.state() == IteratorState::past_the_end &&
+             particles.front().particles.empty() &&
+             particles.back().cell.state() == IteratorState::past_the_end &&
+             particles.back().particles.empty() &&
+             owned_particles_end->cell.state() == IteratorState::past_the_end &&
+             owned_particles_end->particles.empty(),
+           ExcInternalError());
 
 #ifdef DEBUG
+    // check that no cache element hits the three anchor states in the list of
+    // particles
+    for (const auto &it : cells_to_particle_cache)
+      Assert(it != particles.begin() && it != owned_particles_end &&
+               it != --(particles.end()),
+             ExcInternalError());
+
+    // check that we only have locally owned particles in the first region of
+    // cells; note that we skip the very first anchor element
+    for (auto it = ++particles.begin(); it != owned_particles_end; ++it)
+      Assert(it->cell->is_locally_owned(), ExcInternalError());
+
+    // check that the cache is consistent with the iterators
     std::vector<typename particle_container::iterator> verify_cache(
-      triangulation->n_active_cells());
-    for (auto it = owned_particles.begin(); it != owned_particles.end(); ++it)
-      verify_cache[it->cell_iterator->active_cell_index()] = it;
-    for (auto it = ghost_particles.begin(); it != ghost_particles.end(); ++it)
-      verify_cache[it->cell_iterator->active_cell_index()] = it;
+      triangulation->n_active_cells(), particles.end());
+    for (auto it = particles.begin(); it != particles.end(); ++it)
+      if (!it->particles.empty())
+        verify_cache[it->cell->active_cell_index()] = it;
 
     for (unsigned int i = 0; i < verify_cache.size(); ++i)
       Assert(verify_cache[i] == cells_to_particle_cache[i], ExcInternalError());
 #endif
 
-    number_of_locally_owned_particles = result_owned[2];
-    result_owned[0] = std::max(result_owned[0], result_ghosted[0]);
-    result_owned[1] = std::max(result_owned[1], result_ghosted[1]);
+    // now compute local result with the function above and then compute the
+    // collective results
+    number_of_locally_owned_particles = 0;
+
+    types::particle_index result[2] = {};
+    for (const auto &particles_in_cell : particles)
+      {
+        const types::particle_index n_particles_in_cell =
+          particles_in_cell.particles.size();
+
+        // local_max_particles_per_cell
+        result[0] = std::max(result[0], n_particles_in_cell);
+
+        // number of locally owned particles
+        if (n_particles_in_cell > 0 &&
+            particles_in_cell.cell->is_locally_owned())
+          number_of_locally_owned_particles += n_particles_in_cell;
+
+        // local_max_particle_index
+        for (const auto &particle : particles_in_cell.particles)
+          result[1] = std::max(result[1], property_pool->get_id(particle));
+      }
 
     global_number_of_particles =
       dealii::Utilities::MPI::sum(number_of_locally_owned_particles,
@@ -286,13 +372,10 @@ namespace Particles
       }
     else
       {
-        Utilities::MPI::max(
-          ArrayView<const types::particle_index>(result_owned.data(), 2),
-          triangulation->get_communicator(),
-          make_array_view(result_owned.data(), result_owned.data() + 2));
+        Utilities::MPI::max(result, triangulation->get_communicator(), result);
 
-        next_free_particle_index      = result_owned[1] + 1;
-        global_max_particles_per_cell = result_owned[0];
+        next_free_particle_index      = result[1] + 1;
+        global_max_particles_per_cell = result[0];
       }
   }
 
@@ -310,7 +393,7 @@ namespace Particles
     if (cell->is_artificial() == false)
       {
         return cells_to_particle_cache[cell->active_cell_index()] !=
-                   typename particle_container::iterator() ?
+                   particles.end() ?
                  cells_to_particle_cache[cell->active_cell_index()]
                    ->particles.size() :
                  0;
@@ -347,35 +430,23 @@ namespace Particles
 
     if (cell->is_artificial() == false)
       {
-        if (cells_to_particle_cache[active_cell_index] ==
-            typename particle_container::iterator())
+        if (cells_to_particle_cache[active_cell_index] == particles.end())
           {
             return boost::make_iterator_range(
-              particle_iterator(
-                owned_particles, owned_particles.end(), *property_pool, 0),
-              particle_iterator(
-                owned_particles, owned_particles.end(), *property_pool, 0));
+              particle_iterator(particles.begin(), *property_pool, 0),
+              particle_iterator(particles.begin(), *property_pool, 0));
           }
         else
           {
-            const typename particle_container::iterator particles =
-              cells_to_particle_cache[active_cell_index];
-            typename particle_container::iterator next = particles;
-            ++next;
-            if (cell->is_locally_owned())
-              {
-                return boost::make_iterator_range(
-                  particle_iterator(
-                    owned_particles, particles, *property_pool, 0),
-                  particle_iterator(owned_particles, next, *property_pool, 0));
-              }
-            else
-              {
-                return boost::make_iterator_range(
-                  particle_iterator(
-                    ghost_particles, particles, *property_pool, 0),
-                  particle_iterator(ghost_particles, next, *property_pool, 0));
-              }
+            const typename particle_container::iterator
+              particles_in_current_cell =
+                cells_to_particle_cache[active_cell_index];
+            typename particle_container::iterator particles_in_next_cell =
+              particles_in_current_cell;
+            ++particles_in_next_cell;
+            return boost::make_iterator_range(
+              particle_iterator(particles_in_current_cell, *property_pool, 0),
+              particle_iterator(particles_in_next_cell, *property_pool, 0));
           }
       }
     else
@@ -419,12 +490,8 @@ namespace Particles
       }
     else
       {
-        if (owned_cell)
-          owned_particles.erase(particle->particles_in_cell);
-        else
-          ghost_particles.erase(particle->particles_in_cell);
-        cells_to_particle_cache[cell->active_cell_index()] =
-          typename particle_container::iterator();
+        particles.erase(particle->particles_in_cell);
+        cells_to_particle_cache[cell->active_cell_index()] = particles.end();
       }
   }
 
@@ -436,24 +503,47 @@ namespace Particles
     const std::vector<ParticleHandler<dim, spacedim>::particle_iterator>
       &particles_to_remove)
   {
-    // sort particles in reverse order on each cell to ensure the removal with
-    // the other function is valid
-    std::vector<ParticleHandler<dim, spacedim>::particle_iterator> copy(
-      particles_to_remove);
-    std::sort(copy.begin(),
-              copy.end(),
-              [&](const ParticleHandler<dim, spacedim>::particle_iterator &a,
-                  const ParticleHandler<dim, spacedim>::particle_iterator &b) {
-                return a->particles_in_cell->cell_iterator <
-                         b->particles_in_cell->cell_iterator ||
-                       (a->particles_in_cell->cell_iterator ==
-                          b->particles_in_cell->cell_iterator &&
-                        a->particle_index_within_cell >
-                          b->particle_index_within_cell);
-              });
-
-    for (const auto &particle : copy)
-      remove_particle(particle);
+    // We need to remove particles backwards on each cell to keep the particle
+    // iterators alive as we keep removing particles on the same cell. To
+    // ensure that this is safe, we either check if we already have sorted
+    // iterators or if we need to manually sort
+    const auto check_greater = [](const particle_iterator &a,
+                                  const particle_iterator &b) {
+      return a->particles_in_cell->cell > b->particles_in_cell->cell ||
+             (a->particles_in_cell->cell == b->particles_in_cell->cell &&
+              a->particle_index_within_cell > b->particle_index_within_cell);
+    };
+
+    bool particles_are_sorted = true;
+    auto previous             = particles_to_remove.begin();
+    for (auto next = previous; next != particles_to_remove.end(); ++next)
+      {
+        if (check_greater(*previous, *next))
+          {
+            particles_are_sorted = false;
+            break;
+          }
+        previous = next;
+      }
+    if (particles_are_sorted)
+      {
+        // pass along backwards in array
+        for (auto it = particles_to_remove.rbegin();
+             it != particles_to_remove.rend();
+             ++it)
+          remove_particle(*it);
+      }
+    else
+      {
+        std::vector<ParticleHandler<dim, spacedim>::particle_iterator>
+          sorted_particles(particles_to_remove);
+        std::sort(sorted_particles.begin(),
+                  sorted_particles.end(),
+                  check_greater);
+
+        for (const auto &particle : sorted_particles)
+          remove_particle(particle);
+      }
 
     update_cached_numbers();
   }
@@ -476,25 +566,31 @@ namespace Particles
 
 
   template <int dim, int spacedim>
-  void
+  typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::insert_particle(
     const typename PropertyPool<dim, spacedim>::Handle          handle,
-    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
-    particle_container &                                        particles)
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell)
   {
     const unsigned int active_cell_index = cell->active_cell_index();
     typename particle_container::iterator &cache =
       cells_to_particle_cache[active_cell_index];
-    if (cache == typename particle_container::iterator())
-      cache = particles.emplace(
-        particles.end(),
-        std::vector<typename PropertyPool<dim, spacedim>::Handle>{handle},
-        cell);
+    if (cache == particles.end())
+      {
+        const typename particle_container::iterator insert_position =
+          cell->is_locally_owned() ? owned_particles_end : --particles.end();
+        cache = particles.emplace(
+          insert_position,
+          std::vector<typename PropertyPool<dim, spacedim>::Handle>{handle},
+          cell);
+      }
     else
       {
         cache->particles.push_back(handle);
-        Assert(cache->cell_iterator == cell, ExcInternalError());
+        Assert(cache->cell == cell, ExcInternalError());
       }
+    return particle_iterator(cache,
+                             *property_pool,
+                             cache->particles.size() - 1);
   }
 
 
@@ -512,13 +608,8 @@ namespace Particles
            ExcMessage("You tried to insert particles into a cell that is not "
                       "locally owned. This is not supported."));
 
-    insert_particle(property_pool->register_particle(), cell, owned_particles);
-    const typename particle_container::iterator &cache =
-      cells_to_particle_cache[cell->active_cell_index()];
-    particle_iterator particle_it(owned_particles,
-                                  cache,
-                                  *property_pool,
-                                  cache->particles.size() - 1);
+    particle_iterator particle_it =
+      insert_particle(property_pool->register_particle(), cell);
 
     data = particle_it->read_particle_data_from_memory(data);
 
@@ -546,13 +637,8 @@ namespace Particles
            ExcMessage("You tried to insert particles into a cell that is not "
                       "locally owned. This is not supported."));
 
-    insert_particle(property_pool->register_particle(), cell, owned_particles);
-    const typename particle_container::iterator &cache =
-      cells_to_particle_cache[cell->active_cell_index()];
-    particle_iterator particle_it(owned_particles,
-                                  cache,
-                                  *property_pool,
-                                  cache->particles.size() - 1);
+    particle_iterator particle_it =
+      insert_particle(property_pool->register_particle(), cell);
 
     particle_it->set_location(position);
     particle_it->set_reference_location(reference_position);
@@ -1310,7 +1396,7 @@ namespace Particles
               auto &old =
                 out_particle->particles_in_cell
                   ->particles[out_particle->particle_index_within_cell];
-              insert_particle(old, current_cell, owned_particles);
+              insert_particle(old, current_cell);
 
               // Avoid deallocating the memory of this particle
               old = PropertyPool<dim, spacedim>::invalid_handle;
@@ -1332,7 +1418,7 @@ namespace Particles
       {
         if (dealii::Utilities::MPI::n_mpi_processes(
               parallel_triangulation->get_communicator()) > 1)
-          send_recv_particles(moved_particles, owned_particles, moved_cells);
+          send_recv_particles(moved_particles, moved_cells);
       }
 #endif
 
@@ -1344,13 +1430,7 @@ namespace Particles
     unsorted_handles.reserve(property_pool->n_registered_slots());
 
     typename PropertyPool<dim, spacedim>::Handle sorted_handle = 0;
-    for (auto &particles_in_cell : owned_particles)
-      for (auto &particle : particles_in_cell.particles)
-        {
-          unsorted_handles.push_back(particle);
-          particle = sorted_handle++;
-        }
-    for (auto &particles_in_cell : ghost_particles)
+    for (auto &particles_in_cell : particles)
       for (auto &particle : particles_in_cell.particles)
         {
           unsorted_handles.push_back(particle);
@@ -1387,11 +1467,10 @@ namespace Particles
     // Clear ghost particles and their properties
     for (const auto &cell : triangulation->active_cell_iterators())
       if (cell->is_ghost() &&
-          cells_to_particle_cache[cell->active_cell_index()] !=
-            typename particle_container::iterator())
+          cells_to_particle_cache[cell->active_cell_index()] != particles.end())
         {
-          Assert(cells_to_particle_cache[cell->active_cell_index()]
-                     ->cell_iterator == cell,
+          Assert(cells_to_particle_cache[cell->active_cell_index()]->cell ==
+                   cell,
                  ExcInternalError());
           // Clear particle properties
           for (auto &ghost_particle :
@@ -1399,10 +1478,8 @@ namespace Particles
             property_pool->deregister_particle(ghost_particle);
 
           // Clear particles themselves
-          ghost_particles.erase(
-            cells_to_particle_cache[cell->active_cell_index()]);
-          cells_to_particle_cache[cell->active_cell_index()] =
-            typename particle_container::iterator();
+          particles.erase(cells_to_particle_cache[cell->active_cell_index()]);
+          cells_to_particle_cache[cell->active_cell_index()] = particles.end();
         }
 
     // Clear ghost particles cache and invalidate it
@@ -1450,7 +1527,6 @@ namespace Particles
 
     send_recv_particles(
       ghost_particles_cache.ghost_particles_by_domain,
-      ghost_particles,
       std::map<
         types::subdomain_id,
         std::vector<
@@ -1487,7 +1563,7 @@ namespace Particles
 
 
     send_recv_particles_properties_and_location(
-      ghost_particles_cache.ghost_particles_by_domain, ghost_particles);
+      ghost_particles_cache.ghost_particles_by_domain);
 #endif
   }
 
@@ -1498,8 +1574,7 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::send_recv_particles(
     const std::map<types::subdomain_id, std::vector<particle_iterator>>
-      &                 particles_to_send,
-    particle_container &received_particles,
+      &particles_to_send,
     const std::map<
       types::subdomain_id,
       std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
@@ -1761,15 +1836,12 @@ namespace Particles
         const typename Triangulation<dim, spacedim>::active_cell_iterator cell =
           triangulation->create_cell_iterator(id);
 
-        insert_particle(property_pool->register_particle(),
-                        cell,
-                        received_particles);
+        insert_particle(property_pool->register_particle(), cell);
         const typename particle_container::iterator &cache =
           cells_to_particle_cache[cell->active_cell_index()];
-        Assert(cache->cell_iterator == cell, ExcInternalError());
+        Assert(cache->cell == cell, ExcInternalError());
 
-        particle_iterator particle_it(received_particles,
-                                      cache,
+        particle_iterator particle_it(cache,
                                       *property_pool,
                                       cache->particles.size() - 1);
 
@@ -1797,8 +1869,7 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::send_recv_particles_properties_and_location(
     const std::map<types::subdomain_id, std::vector<particle_iterator>>
-      &                 particles_to_send,
-    particle_container &updated_particles)
+      &particles_to_send)
   {
     const auto parallel_triangulation =
       dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
@@ -1890,13 +1961,12 @@ namespace Particles
         recv_data_it =
           recv_particle->read_particle_data_from_memory(recv_data_it);
 
-        Assert(recv_particle->particles_in_cell->cell_iterator->is_ghost(),
+        Assert(recv_particle->particles_in_cell->cell->is_ghost(),
                ExcInternalError());
 
         if (load_callback)
           recv_data_it = load_callback(
-            particle_iterator(updated_particles,
-                              recv_particle->particles_in_cell,
+            particle_iterator(recv_particle->particles_in_cell,
                               *property_pool,
                               recv_particle->particle_index_within_cell),
             recv_data_it);
@@ -1987,7 +2057,8 @@ namespace Particles
     // Resize the container if it is possible without
     // transferring particles
     if (number_of_locally_owned_particles == 0)
-      cells_to_particle_cache.resize(triangulation->n_active_cells());
+      cells_to_particle_cache.resize(triangulation->n_active_cells(),
+                                     particles.end());
   }
 
 
@@ -2159,7 +2230,6 @@ namespace Particles
 
             for (unsigned int i = 0; i < n_particles; ++i)
               stored_particles_on_cell.push_back(particle_iterator(
-                owned_particles,
                 cells_to_particle_cache[cell->active_cell_index()],
                 *property_pool,
                 i));
@@ -2180,8 +2250,8 @@ namespace Particles
                 const typename particle_container::iterator &cache =
                   cells_to_particle_cache[child->active_cell_index()];
                 for (unsigned int i = 0; i < n_particles; ++i)
-                  stored_particles_on_cell.push_back(particle_iterator(
-                    owned_particles, cache, *property_pool, i));
+                  stored_particles_on_cell.push_back(
+                    particle_iterator(cache, *property_pool, i));
               }
           }
           break;
@@ -2203,25 +2273,30 @@ namespace Particles
     const typename Triangulation<dim, spacedim>::CellStatus         status,
     const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
   {
+    if (data_range.begin() == data_range.end())
+      return;
+
     const auto cell_to_store_particles =
       (status != parallel::TriangulationBase<dim, spacedim>::CELL_REFINE) ?
         cell :
         cell->child(0);
 
     // deserialize particles and insert into local storage
-    {
-      const void *data = static_cast<const void *>(&(*data_range.begin()));
-
-      while (data < &(*data_range.end()))
-        insert_particle(data, cell_to_store_particles);
-
-      Assert(
-        data == &(*data_range.end()),
-        ExcMessage(
-          "The particle data could not be deserialized successfully. "
-          "Check that when deserializing the particles you expect the same "
-          "number of properties that were serialized."));
-    }
+    if (data_range.begin() != data_range.end())
+      {
+        const void *data = static_cast<const void *>(&(*data_range.begin()));
+        const void *end  = static_cast<const void *>(
+          &(*data_range.begin()) + (data_range.end() - data_range.begin()));
+
+        while (data < end)
+          insert_particle(data, cell_to_store_particles);
+
+        Assert(data == end,
+               ExcMessage(
+                 "The particle data could not be deserialized successfully. "
+                 "Check that when deserializing the particles you expect "
+                 "the same number of properties that were serialized."));
+      }
 
     auto loaded_particles_on_cell = particles_in_cell(cell_to_store_particles);
 
@@ -2252,15 +2327,15 @@ namespace Particles
           {
             // we need to find the correct child to store the particles and
             // their reference location has changed
-            const typename particle_container::iterator &cache =
+            typename particle_container::iterator &cache =
               cells_to_particle_cache[cell_to_store_particles
                                         ->active_cell_index()];
 
+            // make sure that the call above has inserted an entry
+            Assert(cache != particles.end(), ExcInternalError());
+
             // Cannot use range-based loop, because number of particles in cell
             // is going to change
-            if (cache == typename particle_container::iterator())
-              break;
-
             auto particle = loaded_particles_on_cell.begin();
             for (unsigned int i = 0; i < cache->particles.size();)
               {
@@ -2286,9 +2361,7 @@ namespace Particles
                             // redo the loop; otherwise move on to next particle
                             if (child_index != 0)
                               {
-                                insert_particle(cache->particles[i],
-                                                child,
-                                                owned_particles);
+                                insert_particle(cache->particles[i], child);
 
                                 cache->particles[i] = cache->particles.back();
                                 cache->particles.resize(
@@ -2306,6 +2379,12 @@ namespace Particles
                       {}
                   }
               }
+            // clean up in case child 0 has no particle left
+            if (cache->particles.empty())
+              {
+                particles.erase(cache);
+                cache = particles.end();
+              }
           }
           break;
 

diff --git a/tests/particles/particle_iterator_02.cc b/tests/particles/particle_iterator_02.cc
index 17cc1c082f82e48bb362e1330d40a2ec89c446e3..0fce6bd3e6bb6c9c1b212b2a9c855132e31aeb98 100644 (file)
--- a/tests/particles/particle_iterator_02.cc
+++ b/tests/particles/particle_iterator_02.cc
@@ -68,6 +68,10 @@ test()
     particle_container.emplace_back(
       std::vector<typename Particles::PropertyPool<dim>::Handle>(1, particle),
       tr.begin_active());
+    particle_container.emplace_front(
+      std::vector<typename Particles::PropertyPool<dim>::Handle>(), tr.end());
+    particle_container.emplace_back(
+      std::vector<typename Particles::PropertyPool<dim>::Handle>(), tr.end());
 
     pool.set_location(particle, position);
     pool.set_reference_location(particle, reference_position);
@@ -78,18 +82,15 @@ test()
               properties.end(),
               particle_properties.begin());
 
-    Particles::ParticleIterator<dim> particle_begin(particle_container,
-                                                    particle_container.begin(),
-                                                    pool,
-                                                    0);
-    Particles::ParticleIterator<dim> particle_end(particle_container,
-                                                  particle_container.end(),
+    Particles::ParticleIterator<dim> particle_begin(
+      ++particle_container.begin(), pool, 0);
+    Particles::ParticleIterator<dim> particle_end(--particle_container.end(),
                                                   pool,
                                                   0);
     Particles::ParticleIterator<dim> particle_nonexistent1(
-      particle_container, particle_container.begin(), pool, 1);
+      ++particle_container.begin(), pool, 1);
     Particles::ParticleIterator<dim> particle_nonexistent2(
-      particle_container, particle_container.end(), pool, 1);
+      --particle_container.end(), pool, 1);
     Particles::ParticleIterator<dim> particle_invalid;
 
     Assert(particle_begin->state() == IteratorState::valid, ExcInternalError());