From: Rene Gassmoeller <rene.gassmoeller@mailbox.org>
Date: Mon, 7 Jun 2021 23:03:53 +0000 (-0400)
Subject: simplify particle data structure
X-Git-Tag: v9.4.0-rc1~1256^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F12409%2Fhead;p=dealii.git

simplify particle data structure
---

diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h
index 3057232ecf..d2c7455e1f 100644
--- a/include/deal.II/particles/particle_accessor.h
+++ b/include/deal.II/particles/particle_accessor.h
@@ -24,6 +24,7 @@
 #include <deal.II/grid/tria_iterator_base.h>
 
 #include <deal.II/particles/particle.h>
+#include <deal.II/particles/property_pool.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -48,7 +49,7 @@ namespace Particles
      * A type for the storage container for particles.
      */
     using particle_container =
-      std::vector<std::vector<Particle<dim, spacedim>>>;
+      std::vector<std::vector<typename PropertyPool<dim, spacedim>::Handle>>;
 
     /**
      * @copydoc Particle::write_particle_data_to_memory
@@ -118,20 +119,16 @@ namespace Particles
     get_reference_location() const;
 
     /**
-     * Return the ID number of this particle.
+     * Set the ID number of this particle.
      */
-    types::particle_index
-    get_id() const;
+    void
+    set_id(const types::particle_index &new_id);
 
     /**
-     * Tell the particle where to store its properties (even if it does not
-     * own properties). Usually this is only done once per particle, but
-     * since the particle generator does not know about the properties
-     * we want to do it not at construction time. Another use for this
-     * function is after particle transfer to a new process.
+     * Return the ID number of this particle.
      */
-    void
-    set_property_pool(PropertyPool<dim, spacedim> &property_pool);
+    types::particle_index
+    get_id() const;
 
     /**
      * Return whether this particle has a valid property pool and a valid
@@ -223,12 +220,41 @@ namespace Particles
       const Triangulation<dim, spacedim> &triangulation) const;
 
     /**
-     * Serialize the contents of this class using the [BOOST serialization
+     * Write the data of this object to a stream for the purpose of
+     * serialization using the [BOOST serialization
+     * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html).
+     */
+    template <class Archive>
+    void
+    save(Archive &ar, const unsigned int version) const;
+
+    /**
+     * Read the data of this object from a stream for the purpose of
+     * serialization using the [BOOST serialization
      * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html).
+     * Note that in order to store the properties correctly, the property pool
+     * of this particle has to be known at the time of reading, i.e.
+     * set_property_pool() has to have been called, before this function is
+     * called.
      */
     template <class Archive>
     void
-    serialize(Archive &ar, const unsigned int version);
+    load(Archive &ar, const unsigned int version);
+
+#ifdef DOXYGEN
+    /**
+     * Write and read the data of this object from a stream for the purpose
+     * of serialization using the [BOOST serialization
+     * library](https://www.boost.org/doc/libs/1_74_0/libs/serialization/doc/index.html).
+     */
+    template <class Archive>
+    void
+    serialize(Archive &archive, const unsigned int version);
+#else
+    // This macro defines the serialize() method that is compatible with
+    // the templated save() and load() method that have been implemented.
+    BOOST_SERIALIZATION_SPLIT_MEMBER()
+#endif
 
     /**
      * Advance the ParticleAccessor to the next particle.
@@ -273,7 +299,8 @@ namespace Particles
      * friend classes.
      */
     ParticleAccessor(
-      const particle_container &particles,
+      const particle_container &         particles,
+      const PropertyPool<dim, spacedim> &property_pool,
       const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
       const unsigned int particle_index_within_cell);
 
@@ -282,14 +309,14 @@ namespace Particles
      * internal structure may change this is not intended for public use and
      * only a convenience function for internal purposes.
      */
-    const Particle<dim, spacedim> &
-    get_particle() const;
+    const typename PropertyPool<dim, spacedim>::Handle &
+    get_handle() const;
 
     /**
      * Non-@p const version of the function with the same name above.
      */
-    Particle<dim, spacedim> &
-    get_particle();
+    typename PropertyPool<dim, spacedim>::Handle &
+    get_handle();
 
     /**
      * A pointer to the container that stores the particles. Obviously,
@@ -297,6 +324,11 @@ namespace Particles
      */
     particle_container *particles;
 
+    /**
+     * A pointer to the property pool that stores the actual particle data.
+     */
+    PropertyPool<dim, spacedim> *property_pool;
+
     /**
      * Cell iterator to the cell of the current particle.
      */
@@ -326,11 +358,57 @@ namespace Particles
 
   template <int dim, int spacedim>
   template <class Archive>
-  void
-  ParticleAccessor<dim, spacedim>::serialize(Archive &          ar,
-                                             const unsigned int version)
+  inline void
+  ParticleAccessor<dim, spacedim>::load(Archive &ar, const unsigned int)
+  {
+    unsigned int n_properties = 0;
+
+    Point<spacedim>       location;
+    Point<dim>            reference_location;
+    types::particle_index id;
+    ar &location &reference_location &id &n_properties;
+
+    set_location(location);
+    set_reference_location(reference_location);
+    set_id(id);
+
+    if (n_properties > 0)
+      {
+        ArrayView<double> properties(get_properties());
+        Assert(
+          properties.size() == n_properties,
+          ExcMessage(
+            "This particle was serialized with " +
+            std::to_string(n_properties) +
+            " properties, but the new property handler provides space for " +
+            std::to_string(properties.size()) +
+            " properties. Deserializing a particle only works for matching property sizes."));
+
+        ar &boost::serialization::make_array(properties.data(), n_properties);
+      }
+  }
+
+
+
+  template <int dim, int spacedim>
+  template <class Archive>
+  inline void
+  ParticleAccessor<dim, spacedim>::save(Archive &ar, const unsigned int) const
   {
-    return get_particle().serialize(ar, version);
+    unsigned int n_properties = 0;
+    if ((property_pool != nullptr) &&
+        (get_handle() != PropertyPool<dim, spacedim>::invalid_handle))
+      n_properties = get_properties().size();
+
+    Point<spacedim>       location           = get_location();
+    Point<dim>            reference_location = get_reference_location();
+    types::particle_index id                 = get_id();
+
+    ar &location &reference_location &id &n_properties;
+
+    if (n_properties > 0)
+      ar &boost::serialization::make_array(get_properties().data(),
+                                           n_properties);
   }
 
 
@@ -339,6 +417,7 @@ namespace Particles
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor()
     : particles(nullptr)
+    , property_pool(nullptr)
     , cell()
     , active_cell_index(numbers::invalid_unsigned_int)
     , particle_index_within_cell(numbers::invalid_unsigned_int)
@@ -348,10 +427,12 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor(
-    const particle_container &particles,
+    const particle_container &         particles,
+    const PropertyPool<dim, spacedim> &property_pool,
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
     const unsigned int particle_index_within_cell)
     : particles(const_cast<particle_container *>(&particles))
+    , property_pool(const_cast<PropertyPool<dim, spacedim> *>(&property_pool))
     , cell(cell)
     , particle_index_within_cell(particle_index_within_cell)
   {
@@ -372,7 +453,32 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().read_particle_data_from_memory(data);
+    const types::particle_index *id_data =
+      static_cast<const types::particle_index *>(data);
+    set_id(*id_data++);
+    const double *pdata = reinterpret_cast<const double *>(id_data);
+
+    Point<spacedim> location;
+    for (unsigned int i = 0; i < spacedim; ++i)
+      location(i) = *pdata++;
+    set_location(location);
+
+    Point<dim> reference_location;
+    for (unsigned int i = 0; i < dim; ++i)
+      reference_location(i) = *pdata++;
+    set_reference_location(reference_location);
+
+    // See if there are properties to load
+    if (has_properties())
+      {
+        const ArrayView<double> particle_properties =
+          property_pool->get_properties(get_handle());
+        const unsigned int size = particle_properties.size();
+        for (unsigned int i = 0; i < size; ++i)
+          particle_properties[i] = *pdata++;
+      }
+
+    return static_cast<const void *>(pdata);
   }
 
 
@@ -384,7 +490,29 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().write_particle_data_to_memory(data);
+    types::particle_index *id_data = static_cast<types::particle_index *>(data);
+    *id_data                       = get_id();
+    ++id_data;
+    double *pdata = reinterpret_cast<double *>(id_data);
+
+    // Write location
+    for (unsigned int i = 0; i < spacedim; ++i, ++pdata)
+      *pdata = get_location()[i];
+
+    // Write reference location
+    for (unsigned int i = 0; i < dim; ++i, ++pdata)
+      *pdata = get_reference_location()[i];
+
+    // Write properties
+    if (has_properties())
+      {
+        const ArrayView<double> particle_properties =
+          property_pool->get_properties(get_handle());
+        for (unsigned int i = 0; i < particle_properties.size(); ++i, ++pdata)
+          *pdata = particle_properties[i];
+      }
+
+    return static_cast<void *>(pdata);
   }
 
 
@@ -395,7 +523,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    get_particle().set_location(new_loc);
+    property_pool->set_location(get_handle(), new_loc);
   }
 
 
@@ -406,7 +534,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().get_location();
+    return property_pool->get_location(get_handle());
   }
 
 
@@ -418,7 +546,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    get_particle().set_reference_location(new_loc);
+    property_pool->set_reference_location(get_handle(), new_loc);
   }
 
 
@@ -429,7 +557,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().get_reference_location();
+    return property_pool->get_reference_location(get_handle());
   }
 
 
@@ -440,19 +568,18 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().get_id();
+    return property_pool->get_id(get_handle());
   }
 
 
 
   template <int dim, int spacedim>
   inline void
-  ParticleAccessor<dim, spacedim>::set_property_pool(
-    PropertyPool<dim, spacedim> &new_property_pool)
+  ParticleAccessor<dim, spacedim>::set_id(const types::particle_index &new_id)
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    get_particle().set_property_pool(new_property_pool);
+    property_pool->set_id(get_handle(), new_id);
   }
 
 
@@ -463,7 +590,14 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().has_properties();
+    // Particles always have a property pool associated with them,
+    // but we can access properties only if there is a valid handle.
+    // The only way a particle can have no valid handle if it has
+    // been moved-from -- but that leaves an object in an invalid
+    // state, and so we can just assert that that can't be the case.
+    Assert((get_handle() != PropertyPool<dim, spacedim>::invalid_handle),
+           ExcInternalError());
+    return (property_pool->n_properties_per_slot() > 0);
   }
 
 
@@ -475,7 +609,8 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    get_particle().set_properties(new_properties);
+    set_properties(
+      ArrayView<const double>(new_properties.data(), new_properties.size()));
   }
 
 
@@ -487,7 +622,22 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    get_particle().set_properties(new_properties);
+    const ArrayView<double> property_values =
+      property_pool->get_properties(get_handle());
+
+    Assert(new_properties.size() == property_values.size(),
+           ExcMessage(
+             "You are trying to assign properties with an incompatible length. "
+             "The particle has space to store " +
+             std::to_string(property_values.size()) +
+             " properties, but you are trying to assign " +
+             std::to_string(new_properties.size()) +
+             " properties. This is not allowed."));
+
+    if (property_values.size() > 0)
+      std::copy(new_properties.begin(),
+                new_properties.end(),
+                property_values.begin());
   }
 
 
@@ -498,7 +648,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().get_properties();
+    return property_pool->get_properties(get_handle());
   }
 
 
@@ -532,7 +682,7 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().get_properties();
+    return property_pool->get_properties(get_handle());
   }
 
 
@@ -543,7 +693,14 @@ namespace Particles
   {
     Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return get_particle().serialized_size_in_bytes();
+    std::size_t size = sizeof(get_id()) + sizeof(get_location()) +
+                       sizeof(get_reference_location());
+
+    if (has_properties())
+      {
+        size += sizeof(double) * get_properties().size();
+      }
+    return size;
   }
 
 
@@ -621,7 +778,8 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::
   operator==(const ParticleAccessor<dim, spacedim> &other) const
   {
-    return (particles == other.particles) && (cell == other.cell) &&
+    return (particles == other.particles) &&
+           (property_pool == other.property_pool) && (cell == other.cell) &&
            (particle_index_within_cell == other.particle_index_within_cell);
   }
 
@@ -631,7 +789,8 @@ namespace Particles
   inline IteratorState::IteratorStates
   ParticleAccessor<dim, spacedim>::state() const
   {
-    if (particles != nullptr && cell.state() == IteratorState::valid &&
+    if (particles != nullptr && property_pool != nullptr &&
+        cell.state() == IteratorState::valid &&
         particle_index_within_cell < (*particles)[active_cell_index].size())
       return IteratorState::valid;
     else if (particles != nullptr &&
@@ -647,8 +806,8 @@ namespace Particles
 
 
   template <int dim, int spacedim>
-  inline Particle<dim, spacedim> &
-  ParticleAccessor<dim, spacedim>::get_particle()
+  inline typename PropertyPool<dim, spacedim>::Handle &
+  ParticleAccessor<dim, spacedim>::get_handle()
   {
     return (*particles)[active_cell_index][particle_index_within_cell];
   }
@@ -656,8 +815,8 @@ namespace Particles
 
 
   template <int dim, int spacedim>
-  inline const Particle<dim, spacedim> &
-  ParticleAccessor<dim, spacedim>::get_particle() const
+  inline const typename PropertyPool<dim, spacedim>::Handle &
+  ParticleAccessor<dim, spacedim>::get_handle() const
   {
     return (*particles)[active_cell_index][particle_index_within_cell];
   }
diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h
index 7a131922c3..4cc8088af2 100644
--- a/include/deal.II/particles/particle_handler.h
+++ b/include/deal.II/particles/particle_handler.h
@@ -75,7 +75,7 @@ namespace Particles
      * A type for the storage container for particles.
      */
     using particle_container =
-      std::vector<std::vector<Particle<dim, spacedim>>>;
+      std::vector<std::vector<typename PropertyPool<dim, spacedim>::Handle>>;
 
     /**
      * Default constructor.
@@ -97,7 +97,7 @@ namespace Particles
     /**
      * Destructor.
      */
-    virtual ~ParticleHandler() override = default;
+    virtual ~ParticleHandler();
 
     /**
      * Initialize the particle handler. This function does not clear the
@@ -975,8 +975,8 @@ namespace Particles
     void
     send_recv_particles(
       const std::map<types::subdomain_id, std::vector<particle_iterator>>
-        &                                                particles_to_send,
-      std::vector<std::vector<Particle<dim, spacedim>>> &received_particles,
+        &                 particles_to_send,
+      particle_container &received_particles,
       const std::map<
         types::subdomain_id,
         std::vector<
@@ -1007,8 +1007,8 @@ namespace Particles
     void
     send_recv_particles_properties_and_location(
       const std::map<types::subdomain_id, std::vector<particle_iterator>>
-        &                                                particles_to_send,
-      std::vector<std::vector<Particle<dim, spacedim>>> &received_particles);
+        &                 particles_to_send,
+      particle_container &received_particles);
 
 
 #endif
@@ -1067,7 +1067,7 @@ namespace Particles
     for (const auto &cell : triangulation->active_cell_iterators())
       if (cell->is_locally_owned() &&
           particles[cell->active_cell_index()].size() != 0)
-        return particle_iterator(particles, cell, 0);
+        return particle_iterator(particles, *property_pool, cell, 0);
 
     return end();
   }
@@ -1087,7 +1087,10 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end()
   {
-    return particle_iterator(particles, triangulation->end(), 0);
+    return particle_iterator(particles,
+                             *property_pool,
+                             triangulation->end(),
+                             0);
   }
 
 
@@ -1111,7 +1114,7 @@ namespace Particles
     for (const auto &cell : triangulation->active_cell_iterators())
       if (cell->is_locally_owned() == false &&
           ghost_particles[cell->active_cell_index()].size() != 0)
-        return particle_iterator(ghost_particles, cell, 0);
+        return particle_iterator(ghost_particles, *property_pool, cell, 0);
 
     return end_ghost();
   }
@@ -1131,7 +1134,10 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end_ghost()
   {
-    return particle_iterator(ghost_particles, triangulation->end(), 0);
+    return particle_iterator(ghost_particles,
+                             *property_pool,
+                             triangulation->end(),
+                             0);
   }
 
 
diff --git a/include/deal.II/particles/particle_iterator.h b/include/deal.II/particles/particle_iterator.h
index b817880b25..0e81a5e4de 100644
--- a/include/deal.II/particles/particle_iterator.h
+++ b/include/deal.II/particles/particle_iterator.h
@@ -39,6 +39,12 @@ namespace Particles
   class ParticleIterator
   {
   public:
+    /**
+     * A type for the storage container for particles.
+     */
+    using particle_container =
+      std::vector<std::vector<typename PropertyPool<dim, spacedim>::Handle>>;
+
     /**
      * Empty constructor. Such an object is not usable!
      */
@@ -50,7 +56,8 @@ namespace Particles
      * cell.
      */
     ParticleIterator(
-      const std::vector<std::vector<Particle<dim, spacedim>>> &   particles,
+      const particle_container &                                  particles,
+      const PropertyPool<dim, spacedim> &                         property_pool,
       const typename Triangulation<dim, spacedim>::cell_iterator &cell,
       const unsigned int particle_index_within_cell);
 
@@ -152,10 +159,11 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleIterator<dim, spacedim>::ParticleIterator(
-    const std::vector<std::vector<Particle<dim, spacedim>>> &   particles,
+    const particle_container &                                  particles,
+    const PropertyPool<dim, spacedim> &                         property_pool,
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const unsigned int particle_index_within_cell)
-    : accessor(particles, cell, particle_index_within_cell)
+    : accessor(particles, property_pool, cell, particle_index_within_cell)
   {}
 
 
diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc
index 5f5efeda64..ebab1fa080 100644
--- a/source/particles/particle_handler.cc
+++ b/source/particles/particle_handler.cc
@@ -31,7 +31,7 @@ namespace Particles
   {
     template <int dim, int spacedim>
     std::vector<char>
-    pack_particles(std::vector<Particle<dim, spacedim>> &particles)
+    pack_particles(std::vector<ParticleIterator<dim, spacedim>> &particles)
     {
       std::vector<char> buffer;
 
@@ -39,12 +39,12 @@ namespace Particles
         return buffer;
 
       buffer.resize(particles.size() *
-                    particles.front().serialized_size_in_bytes());
+                    particles.front()->serialized_size_in_bytes());
       void *current_data = buffer.data();
 
       for (const auto &particle : particles)
         {
-          current_data = particle.write_particle_data_to_memory(current_data);
+          current_data = particle->write_particle_data_to_memory(current_data);
         }
 
       return buffer;
@@ -132,6 +132,14 @@ namespace Particles
 
 
 
+  template <int dim, int spacedim>
+  ParticleHandler<dim, spacedim>::~ParticleHandler()
+  {
+    clear_particles();
+  }
+
+
+
   template <int dim, int spacedim>
   void
   ParticleHandler<dim, spacedim>::initialize(
@@ -171,7 +179,8 @@ namespace Particles
                *particle_handler.mapping,
                n_properties);
 
-    property_pool->reserve(particle_handler.n_locally_owned_particles());
+    property_pool = std::make_unique<PropertyPool<dim, spacedim>>(
+      *(particle_handler.property_pool));
 
     // copy static members
     global_number_of_particles = particle_handler.global_number_of_particles;
@@ -185,13 +194,6 @@ namespace Particles
     ghost_particles_cache.ghost_particles_by_domain =
       particle_handler.ghost_particles_cache.ghost_particles_by_domain;
     handle = particle_handler.handle;
-
-    // copy dynamic properties
-    for (auto &particle : *this)
-      particle.set_property_pool(*property_pool);
-
-    for (auto ghost = begin_ghost(); ghost != end_ghost(); ++ghost)
-      ghost->set_property_pool(*property_pool);
   }
 
 
@@ -213,10 +215,20 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::clear_particles()
   {
+    for (auto &particles_in_cell : particles)
+      for (auto &particle : particles_in_cell)
+        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
+          property_pool->deregister_particle(particle);
+
     particles.clear();
     if (triangulation != nullptr)
       particles.resize(triangulation->n_active_cells());
 
+    for (auto &particles_in_cell : ghost_particles)
+      for (auto &particle : particles_in_cell)
+        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
+          property_pool->deregister_particle(particle);
+
     ghost_particles.clear();
     if (triangulation != nullptr)
       ghost_particles.resize(triangulation->n_active_cells());
@@ -248,7 +260,8 @@ namespace Particles
         for (const auto &particle : particles_in_cell)
           {
             local_max_particle_index =
-              std::max(local_max_particle_index, particle.get_id());
+              std::max(local_max_particle_index,
+                       property_pool->get_id(particle));
           }
       }
 
@@ -344,13 +357,14 @@ namespace Particles
         if (container[active_cell_index].size() == 0)
           {
             return boost::make_iterator_range(
-              particle_iterator(container, cell, 0),
-              particle_iterator(container, cell, 0));
+              particle_iterator(container, *property_pool, cell, 0),
+              particle_iterator(container, *property_pool, cell, 0));
           }
         else
           {
-            particle_iterator begin(container, cell, 0);
+            particle_iterator begin(container, *property_pool, cell, 0);
             particle_iterator end(container,
+                                  *property_pool,
                                   cell,
                                   container[active_cell_index].size() - 1);
             // end needs to point to the particle after the last one in the
@@ -380,6 +394,15 @@ namespace Particles
 
     particle_container &container = *(particle->particles);
 
+    // if the particle has an invalid handle (e.g. because it has
+    // been duplicated before calling this function) do not try
+    // to deallocate its memory again
+    auto handle = particle->get_handle();
+    if (handle != PropertyPool<dim, spacedim>::invalid_handle)
+      {
+        property_pool->deregister_particle(handle);
+      }
+
     if (container[active_cell_index].size() > 1)
       {
         container[active_cell_index][particle->particle_index_within_cell] =
@@ -428,6 +451,14 @@ namespace Particles
                 particles_to_remove[n_particles_removed]
                     ->particle_index_within_cell == move_from)
               {
+                auto handle =
+                  particles_to_remove[n_particles_removed]->get_handle();
+
+                // if some particles have invalid handles (e.g. because they are
+                // multiple times in the vector, or have been moved from
+                // before), do not try to deallocate their memory again
+                if (handle != PropertyPool<dim, spacedim>::invalid_handle)
+                  property_pool->deregister_particle(handle);
                 ++n_particles_removed;
                 continue;
               }
@@ -459,11 +490,19 @@ namespace Particles
       particles.resize(triangulation->n_active_cells());
 
     const unsigned int active_cell_index = cell->active_cell_index();
-    particles[active_cell_index].push_back(particle);
+    particles[active_cell_index].push_back(property_pool->register_particle());
     particle_iterator particle_it(particles,
+                                  *property_pool,
                                   cell,
                                   particles[active_cell_index].size() - 1);
-    particle_it->set_property_pool(*property_pool);
+
+    particle_it->set_reference_location(particle.get_reference_location());
+    particle_it->set_location(particle.get_location());
+    particle_it->set_id(particle.get_id());
+
+    if (particle.has_properties())
+      particle_it->set_properties(particle.get_properties());
+
 
     ++local_number_of_particles;
 
@@ -483,13 +522,8 @@ namespace Particles
     if (particles.size() == 0)
       particles.resize(triangulation->n_active_cells());
 
-    for (const auto &particle : new_particles)
-      {
-        const unsigned int active_cell_index =
-          particle.first->active_cell_index();
-        particles[active_cell_index].push_back(particle.second);
-        particles[active_cell_index].back().set_property_pool(*property_pool);
-      }
+    for (const auto &cell_and_particle : new_particles)
+      insert_particle(cell_and_particle.second, cell_and_particle.first);
 
     update_cached_numbers();
   }
@@ -554,17 +588,14 @@ namespace Particles
 
     for (unsigned int i = 0; i < cells.size(); ++i)
       {
-        const unsigned int active_cell_index = cells[i]->active_cell_index();
-
         for (unsigned int p = 0; p < local_positions[i].size(); ++p)
           {
-            particles[active_cell_index].push_back(
-              Particle<dim, spacedim>(positions[index_map[i][p]],
-                                      local_positions[i][p],
-                                      local_start_index + index_map[i][p]));
+            const Particle<dim, spacedim> particle(positions[index_map[i][p]],
+                                                   local_positions[i][p],
+                                                   local_start_index +
+                                                     index_map[i][p]);
 
-            particles[active_cell_index].back().set_property_pool(
-              *property_pool);
+            insert_particle(particle, cells[i]);
           }
       }
 
@@ -1059,7 +1090,7 @@ namespace Particles
             for (unsigned int i = 0; i < n_pic; ++i)
               {
                 particles_out_of_cell.push_back(
-                  particle_iterator(particles, cell, i));
+                  particle_iterator(particles, *property_pool, cell, i));
               }
             continue;
           }
@@ -1239,9 +1270,15 @@ namespace Particles
             {
               const unsigned int active_cell_index =
                 current_cell->active_cell_index();
+
               particles[active_cell_index].push_back(
-                std::move(particles[out_particle->active_cell_index]
-                                   [out_particle->particle_index_within_cell]));
+                particles[out_particle->active_cell_index]
+                         [out_particle->particle_index_within_cell]);
+
+              // Avoid deallocating the memory of this particle
+              particles[out_particle->active_cell_index]
+                       [out_particle->particle_index_within_cell] =
+                         PropertyPool<dim, spacedim>::invalid_handle;
             }
           else
             {
@@ -1292,6 +1329,10 @@ namespace Particles
     (void)enable_cache;
 #else
     // First clear the current ghost_particle information
+    for (auto &particles_in_cell : ghost_particles)
+      for (auto &particle : particles_in_cell)
+        if (particle != PropertyPool<dim, spacedim>::invalid_handle)
+          property_pool->deregister_particle(particle);
     ghost_particles.clear();
     ghost_particles.resize(triangulation->n_active_cells());
 
@@ -1653,14 +1694,16 @@ namespace Particles
           triangulation->create_cell_iterator(id);
 
         const unsigned int active_cell_index = cell->active_cell_index();
-        received_particles[active_cell_index].push_back(
-          Particle<dim, spacedim>(recv_data_it, property_pool.get()));
-
+        received_particles[active_cell_index].emplace_back(
+          property_pool->register_particle());
         particle_iterator particle_it(
           received_particles,
+          *property_pool,
           cell,
           received_particles[active_cell_index].size() - 1);
-        particle_it->set_property_pool(*property_pool);
+
+        recv_data_it =
+          particle_it->read_particle_data_from_memory(recv_data_it);
 
         if (load_callback)
           recv_data_it = load_callback(particle_it, recv_data_it);
@@ -1779,6 +1822,7 @@ namespace Particles
         if (load_callback)
           recv_data_it = load_callback(
             particle_iterator(updated_particles,
+                              *property_pool,
                               recv_particle->cell,
                               recv_particle->particle_index_within_cell),
             recv_data_it);
@@ -1914,7 +1958,7 @@ namespace Particles
     const typename Triangulation<dim, spacedim>::cell_iterator &cell,
     const typename Triangulation<dim, spacedim>::CellStatus     status) const
   {
-    std::vector<Particle<dim, spacedim>> stored_particles_on_cell;
+    std::vector<particle_iterator> stored_particles_on_cell;
 
     switch (status)
       {
@@ -1923,7 +1967,13 @@ namespace Particles
           // If the cell persist or is refined store all particles of the
           // current cell.
           {
-            stored_particles_on_cell = particles[cell->active_cell_index()];
+            const unsigned int n_particles =
+              particles[cell->active_cell_index()].size();
+            stored_particles_on_cell.reserve(n_particles);
+
+            for (unsigned int i = 0; i < n_particles; ++i)
+              stored_particles_on_cell.push_back(
+                particle_iterator(particles, *property_pool, cell, i));
           }
           break;
 
@@ -1933,8 +1983,6 @@ namespace Particles
           {
             for (const auto &child : cell->child_iterators())
               {
-                const unsigned int active_cell_index =
-                  child->active_cell_index();
                 const unsigned int n_particles = n_particles_in_cell(child);
 
                 stored_particles_on_cell.reserve(
@@ -1942,7 +1990,7 @@ namespace Particles
 
                 for (unsigned int i = 0; i < n_particles; ++i)
                   stored_particles_on_cell.push_back(
-                    particles[active_cell_index][i]);
+                    particle_iterator(particles, *property_pool, child, i));
               }
           }
           break;
@@ -1969,34 +2017,24 @@ namespace Particles
     std::vector<Particle<dim, spacedim>> loaded_particles_on_cell =
       unpack_particles<dim, spacedim>(data_range, *property_pool);
 
-    // Update the reference to the current property pool for all particles.
-    // This was not stored, because they might be transported across process
-    // domains.
-    for (auto &particle : loaded_particles_on_cell)
-      particle.set_property_pool(*property_pool);
-
     switch (status)
       {
         case parallel::TriangulationBase<dim, spacedim>::CELL_PERSIST:
           {
-            particles[cell->active_cell_index()].insert(
-              particles[cell->active_cell_index()].end(),
-              loaded_particles_on_cell.begin(),
-              loaded_particles_on_cell.end());
+            for (const auto &particle : loaded_particles_on_cell)
+              insert_particle(particle, cell);
           }
           break;
 
         case parallel::TriangulationBase<dim, spacedim>::CELL_COARSEN:
           {
-            const unsigned int active_cell_index = cell->active_cell_index();
-
             for (auto &particle : loaded_particles_on_cell)
               {
                 const Point<dim> p_unit =
                   mapping->transform_real_to_unit_cell(cell,
                                                        particle.get_location());
                 particle.set_reference_location(p_unit);
-                particles[active_cell_index].emplace_back(std::move(particle));
+                insert_particle(particle, cell);
               }
           }
           break;
@@ -2010,9 +2048,7 @@ namespace Particles
                      ++child_index)
                   {
                     const typename Triangulation<dim, spacedim>::cell_iterator
-                                       child = cell->child(child_index);
-                    const unsigned int active_cell_index =
-                      child->active_cell_index();
+                      child = cell->child(child_index);
 
                     try
                       {
@@ -2024,8 +2060,7 @@ namespace Particles
                             particle.set_reference_location(p_unit);
                             // Use std::multimap::emplace_hint to speed up
                             // insertion of particles.
-                            particles[active_cell_index].emplace_back(
-                              std::move(particle));
+                            insert_particle(particle, child);
                             break;
                           }
                       }
diff --git a/source/particles/property_pool.cc b/source/particles/property_pool.cc
index dce1aa3ec2..cd0bffc73a 100644
--- a/source/particles/property_pool.cc
+++ b/source/particles/property_pool.cc
@@ -123,6 +123,13 @@ namespace Particles
       ExcMessage(
         "This handle is invalid and cannot be deallocated. This can happen if the "
         "handle was deallocated already before calling this function."));
+
+    Assert(
+      currently_available_handles.size() < locations.size(),
+      ExcMessage(
+        "Trying to deallocate a particle when none are allocated. This can happen if all "
+        "handles were deallocated already before calling this function."));
+
     currently_available_handles.push_back(handle);
     handle = invalid_handle;
 
diff --git a/tests/particles/particle_handler_08.cc b/tests/particles/particle_handler_08.cc
index e098b28bf1..1be745e584 100644
--- a/tests/particles/particle_handler_08.cc
+++ b/tests/particles/particle_handler_08.cc
@@ -125,23 +125,19 @@ test()
 
     create_regular_particle_distribution(particle_handler, tr);
 
-    for (auto particle = particle_handler.begin();
-         particle != particle_handler.end();
-         ++particle)
-      deallog << "Before refinement particle id " << particle->get_id()
-              << " has first property " << particle->get_properties()[0]
-              << " and second property " << particle->get_properties()[1]
+    for (const auto &particle : particle_handler)
+      deallog << "Before refinement particle id " << particle.get_id()
+              << " has first property " << particle.get_properties()[0]
+              << " and second property " << particle.get_properties()[1]
               << std::endl;
 
     // Check that all particles are moved to children
     tr.refine_global(1);
 
-    for (auto particle = particle_handler.begin();
-         particle != particle_handler.end();
-         ++particle)
-      deallog << "After refinement particle id " << particle->get_id()
-              << " has first property " << particle->get_properties()[0]
-              << " and second property " << particle->get_properties()[1]
+    for (const auto &particle : particle_handler)
+      deallog << "After refinement particle id " << particle.get_id()
+              << " has first property " << particle.get_properties()[0]
+              << " and second property " << particle.get_properties()[1]
               << std::endl;
 
     // Reverse the refinement and check again
@@ -150,12 +146,10 @@ test()
 
     tr.execute_coarsening_and_refinement();
 
-    for (auto particle = particle_handler.begin();
-         particle != particle_handler.end();
-         ++particle)
-      deallog << "After coarsening particle id " << particle->get_id()
-              << " has first property " << particle->get_properties()[0]
-              << " and second property " << particle->get_properties()[1]
+    for (const auto &particle : particle_handler)
+      deallog << "After coarsening particle id " << particle.get_id()
+              << " has first property " << particle.get_properties()[0]
+              << " and second property " << particle.get_properties()[1]
               << std::endl;
   }
 
diff --git a/tests/particles/particle_iterator_01.cc b/tests/particles/particle_iterator_01.cc
index d47d99c639..bb133fc800 100644
--- a/tests/particles/particle_iterator_01.cc
+++ b/tests/particles/particle_iterator_01.cc
@@ -26,6 +26,7 @@
 #include <deal.II/grid/tria.h>
 
 #include <deal.II/particles/particle.h>
+#include <deal.II/particles/particle_handler.h>
 #include <deal.II/particles/particle_iterator.h>
 
 #include "../tests.h"
@@ -39,10 +40,12 @@ test()
     Triangulation<dim> tr;
 
     GridGenerator::hyper_cube(tr);
-    MappingQ<dim> mapping(1);
+    MappingQ<dim>      mapping(1);
+    const unsigned int n_properties_per_particle = 3;
 
-    const unsigned int           n_properties_per_particle = 3;
-    Particles::PropertyPool<dim> pool(n_properties_per_particle);
+    Particles::ParticleHandler<dim> particle_handler(tr,
+                                                     mapping,
+                                                     n_properties_per_particle);
 
     Point<dim> position;
     position(0) = 0.3;
@@ -61,32 +64,26 @@ test()
     std::vector<double> properties = {0.15, 0.45, 0.75};
 
     Particles::Particle<dim> particle(position, reference_position, index);
-    particle.set_property_pool(pool);
-    particle.set_properties(
-      ArrayView<double>(&properties[0], properties.size()));
-
-    std::vector<std::vector<Particles::Particle<dim>>> particle_container;
-    particle_container.resize(1);
 
-    particle_container[0].push_back(particle);
+    // Insert two identical particles
+    Particles::ParticleIterator<dim> particle_it =
+      particle_handler.insert_particle(particle, tr.begin());
+    particle_it->set_properties(
+      ArrayView<double>(&properties[0], properties.size()));
 
-    particle.get_properties()[0] = 0.05;
-    particle_container[0].push_back(particle);
+    particle_it = particle_handler.insert_particle(particle, tr.begin());
+    particle_it->set_properties(
+      ArrayView<double>(&properties[0], properties.size()));
 
-    Particles::ParticleIterator<dim> particle_it(particle_container,
-                                                 tr.begin(),
-                                                 0);
-    Particles::ParticleIterator<dim> particle_end(particle_container,
-                                                  tr.end(),
-                                                  0);
+    // Modify the properties of the second particle
+    particle_it->get_properties()[0] = 0.05;
 
-    for (; particle_it != particle_end; ++particle_it)
+    for (const auto &particle : particle_handler)
       {
-        deallog << "Particle position: " << (*particle_it).get_location()
-                << std::endl
+        deallog << "Particle position: " << particle.get_location() << std::endl
                 << "Particle properties: "
-                << std::vector<double>(particle_it->get_properties().begin(),
-                                       particle_it->get_properties().end())
+                << std::vector<double>(particle.get_properties().begin(),
+                                       particle.get_properties().end())
                 << std::endl;
       }
   }
diff --git a/tests/particles/particle_iterator_02.cc b/tests/particles/particle_iterator_02.cc
index fa4a1132d9..2907dce62c 100644
--- a/tests/particles/particle_iterator_02.cc
+++ b/tests/particles/particle_iterator_02.cc
@@ -25,7 +25,9 @@
 #include <deal.II/grid/tria.h>
 
 #include <deal.II/particles/particle.h>
+#include <deal.II/particles/particle_handler.h>
 #include <deal.II/particles/particle_iterator.h>
+#include <deal.II/particles/property_pool.h>
 
 #include "../tests.h"
 
@@ -59,26 +61,36 @@ test()
 
     std::vector<double> properties = {0.15, 0.45, 0.75};
 
-    Particles::Particle<dim> particle(position, reference_position, index);
-    particle.set_property_pool(pool);
-    particle.set_properties(
-      ArrayView<double>(&properties[0], properties.size()));
-
-    std::vector<std::vector<Particles::Particle<dim>>> particle_container;
+    typename Particles::ParticleHandler<dim>::particle_container
+      particle_container;
     particle_container.resize(1);
 
+    auto particle = pool.register_particle();
     particle_container[0].push_back(particle);
 
+    pool.set_location(particle, position);
+    pool.set_reference_location(particle, reference_position);
+    pool.set_id(particle, index);
+    auto particle_properties = pool.get_properties(particle);
+
+    std::copy(properties.begin(),
+              properties.end(),
+              particle_properties.begin());
+
     Particles::ParticleIterator<dim> particle_begin(particle_container,
+                                                    pool,
                                                     tr.begin(),
                                                     0);
     Particles::ParticleIterator<dim> particle_end(particle_container,
+                                                  pool,
                                                   tr.end(),
                                                   0);
     Particles::ParticleIterator<dim> particle_nonexistent1(particle_container,
+                                                           pool,
                                                            tr.begin(),
                                                            1);
     Particles::ParticleIterator<dim> particle_nonexistent2(particle_container,
+                                                           pool,
                                                            tr.end(),
                                                            1);
     Particles::ParticleIterator<dim> particle_invalid;
@@ -116,6 +128,8 @@ test()
 
     Assert(particle_iterator->state() == IteratorState::past_the_end,
            ExcInternalError());
+
+    pool.deregister_particle(particle);
   }
 
   deallog << "OK" << std::endl;