]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Store particles in a vector
authorRene Gassmoeller <rene.gassmoeller@mailbox.org>
Wed, 9 Dec 2020 23:59:29 +0000 (18:59 -0500)
committerRene Gassmoeller <rene.gassmoeller@mailbox.org>
Sun, 6 Jun 2021 16:00:43 +0000 (12:00 -0400)
18 files changed:
examples/step-68/step-68.cc
examples/step-70/step-70.cc
include/deal.II/particles/particle_accessor.h
include/deal.II/particles/particle_handler.h
include/deal.II/particles/particle_iterator.h
include/deal.II/particles/partitioner.h
include/deal.II/particles/utilities.h
source/particles/particle_handler.cc
source/particles/utilities.cc
tests/particles/data_out_05.cc
tests/particles/data_out_05.with_p4est=true.mpirun=2.output
tests/particles/particle_handler_06.cc
tests/particles/particle_handler_07.with_p4est=true.mpirun=1.output
tests/particles/particle_handler_serial_07.output
tests/particles/particle_iterator_01.cc
tests/particles/particle_iterator_02.cc [new file with mode: 0644]
tests/particles/particle_iterator_02.output [new file with mode: 0644]
tests/particles/properties_03.output

index 8b6ee0597b49b0e8d0223f25f9e4f52f1449288d..5c89f7b854f37752180e97efa92b2c153aaccc0f 100644 (file)
@@ -566,8 +566,7 @@ namespace Step68
     auto particle = particle_handler.begin();
     while (particle != particle_handler.end())
       {
-        const auto cell =
-          particle->get_surrounding_cell(background_triangulation);
+        const auto cell = particle->get_surrounding_cell();
         const auto dh_cell =
           typename DoFHandler<dim>::cell_iterator(*cell, &fluid_dh);
 
index 29390ba2685596a60a748b1edf2090e0743ba5a9..2ac8003840e5d7f20dc2299af9c75a92f3af818b 100644 (file)
@@ -1492,7 +1492,7 @@ namespace Step70
         // the particle itself. We can then assemble the additional
         // terms in the system matrix and the right hand side as we would
         // normally.
-        const auto &cell = particle->get_surrounding_cell(fluid_tria);
+        const auto &cell = particle->get_surrounding_cell();
         const auto &dh_cell =
           typename DoFHandler<spacedim>::cell_iterator(*cell, &fluid_dh);
         dh_cell->get_dof_indices(fluid_dof_indices);
index b3dd8db0f0091c6ec616aa976165e704160b7364..3057232ecff3bb7ef0d7bc29bf3a0da2d205df31 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/base/array_view.h>
 
 #include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator_base.h>
 
 #include <deal.II/particles/particle.h>
 
@@ -43,6 +44,12 @@ namespace Particles
   class ParticleAccessor
   {
   public:
+    /**
+     * A type for the storage container for particles.
+     */
+    using particle_container =
+      std::vector<std::vector<Particle<dim, spacedim>>>;
+
     /**
      * @copydoc Particle::write_particle_data_to_memory
      */
@@ -203,7 +210,15 @@ namespace Particles
      * but the triangulation itself is not stored in the particle this
      * operation requires a reference to the triangulation.
      */
-    typename Triangulation<dim, spacedim>::cell_iterator
+    const typename Triangulation<dim, spacedim>::cell_iterator &
+    get_surrounding_cell() const;
+
+    /**
+     * @deprecated: Deprecated version of the function with the same
+     * name above.
+     */
+    DEAL_II_DEPRECATED
+    const typename Triangulation<dim, spacedim>::cell_iterator &
     get_surrounding_cell(
       const Triangulation<dim, spacedim> &triangulation) const;
 
@@ -239,6 +254,12 @@ namespace Particles
     bool
     operator==(const ParticleAccessor<dim, spacedim> &other) const;
 
+    /**
+     * Return the state of the accessor.
+     */
+    IteratorState::IteratorStates
+    state() const;
+
   private:
     /**
      * Construct an invalid accessor. Such an object is not usable.
@@ -246,29 +267,52 @@ namespace Particles
     ParticleAccessor();
 
     /**
-     * Construct an accessor from a reference to a map and an iterator to the
-     * map. This constructor is `private` so that it can only be accessed by
+     * Construct an accessor from a reference to a container, an iterator to the
+     * current cell, and the particle index within that cell.
+     * This constructor is `private` so that it can only be accessed by
      * friend classes.
      */
     ParticleAccessor(
-      const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
-      const typename std::multimap<internal::LevelInd,
-                                   Particle<dim, spacedim>>::iterator
-        &particle);
+      const particle_container &particles,
+      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+      const unsigned int particle_index_within_cell);
+
+    /**
+     * Returns a reference to the current Particle. Because the
+     * 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;
+
+    /**
+     * Non-@p const version of the function with the same name above.
+     */
+    Particle<dim, spacedim> &
+    get_particle();
 
-  private:
     /**
      * A pointer to the container that stores the particles. Obviously,
      * this accessor is invalidated if the container changes.
      */
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>> *map;
+    particle_container *particles;
 
     /**
-     * An iterator into the container of particles. Obviously,
-     * this accessor is invalidated if the container changes.
+     * Cell iterator to the cell of the current particle.
+     */
+    typename Triangulation<dim, spacedim>::active_cell_iterator cell;
+
+    /**
+     * Index to the cell this particle is stored in at the moment. This could be
+     * read from the member variable cell, but is used in many performance
+     * critical functions and is therefore stored and updated separately.
+     */
+    unsigned int active_cell_index;
+
+    /**
+     * Local index of the particle within its current cell.
      */
-    typename std::multimap<internal::LevelInd,
-                           Particle<dim, spacedim>>::iterator particle;
+    unsigned int particle_index_within_cell;
 
     // Make ParticleIterator a friend to allow it constructing
     // ParticleAccessors.
@@ -286,7 +330,7 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::serialize(Archive &          ar,
                                              const unsigned int version)
   {
-    return particle->second.serialize(ar, version);
+    return get_particle().serialize(ar, version);
   }
 
 
@@ -294,21 +338,30 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor()
-    : map(nullptr)
-    , particle()
+    : particles(nullptr)
+    , cell()
+    , active_cell_index(numbers::invalid_unsigned_int)
+    , particle_index_within_cell(numbers::invalid_unsigned_int)
   {}
 
 
 
   template <int dim, int spacedim>
   inline ParticleAccessor<dim, spacedim>::ParticleAccessor(
-    const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
-    const typename std::multimap<internal::LevelInd,
-                                 Particle<dim, spacedim>>::iterator & particle)
-    : map(const_cast<
-          std::multimap<internal::LevelInd, Particle<dim, spacedim>> *>(&map))
-    , particle(particle)
-  {}
+    const particle_container &particles,
+    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+    const unsigned int particle_index_within_cell)
+    : particles(const_cast<particle_container *>(&particles))
+    , cell(cell)
+    , particle_index_within_cell(particle_index_within_cell)
+  {
+    if (cell.state() == IteratorState::valid)
+      active_cell_index = cell->active_cell_index();
+    else if (cell.state() == IteratorState::past_the_end)
+      active_cell_index = particles.size();
+    else
+      active_cell_index = numbers::invalid_unsigned_int;
+  }
 
 
 
@@ -317,9 +370,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::read_particle_data_from_memory(
     const void *data)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.read_particle_data_from_memory(data);
+    return get_particle().read_particle_data_from_memory(data);
   }
 
 
@@ -329,9 +382,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::write_particle_data_to_memory(
     void *data) const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.write_particle_data_to_memory(data);
+    return get_particle().write_particle_data_to_memory(data);
   }
 
 
@@ -340,9 +393,9 @@ namespace Particles
   inline void
   ParticleAccessor<dim, spacedim>::set_location(const Point<spacedim> &new_loc)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    particle->second.set_location(new_loc);
+    get_particle().set_location(new_loc);
   }
 
 
@@ -351,9 +404,9 @@ namespace Particles
   inline const Point<spacedim> &
   ParticleAccessor<dim, spacedim>::get_location() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.get_location();
+    return get_particle().get_location();
   }
 
 
@@ -363,9 +416,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::set_reference_location(
     const Point<dim> &new_loc)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    particle->second.set_reference_location(new_loc);
+    get_particle().set_reference_location(new_loc);
   }
 
 
@@ -374,9 +427,9 @@ namespace Particles
   inline const Point<dim> &
   ParticleAccessor<dim, spacedim>::get_reference_location() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.get_reference_location();
+    return get_particle().get_reference_location();
   }
 
 
@@ -385,9 +438,9 @@ namespace Particles
   inline types::particle_index
   ParticleAccessor<dim, spacedim>::get_id() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.get_id();
+    return get_particle().get_id();
   }
 
 
@@ -397,9 +450,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::set_property_pool(
     PropertyPool<dim, spacedim> &new_property_pool)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    particle->second.set_property_pool(new_property_pool);
+    get_particle().set_property_pool(new_property_pool);
   }
 
 
@@ -408,9 +461,9 @@ namespace Particles
   inline bool
   ParticleAccessor<dim, spacedim>::has_properties() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.has_properties();
+    return get_particle().has_properties();
   }
 
 
@@ -420,9 +473,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::set_properties(
     const std::vector<double> &new_properties)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    particle->second.set_properties(new_properties);
+    get_particle().set_properties(new_properties);
   }
 
 
@@ -432,9 +485,9 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::set_properties(
     const ArrayView<const double> &new_properties)
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    particle->second.set_properties(new_properties);
+    get_particle().set_properties(new_properties);
   }
 
 
@@ -443,22 +496,31 @@ namespace Particles
   inline const ArrayView<const double>
   ParticleAccessor<dim, spacedim>::get_properties() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
+
+    return get_particle().get_properties();
+  }
+
+
+
+  template <int dim, int spacedim>
+  inline const typename Triangulation<dim, spacedim>::cell_iterator &
+  ParticleAccessor<dim, spacedim>::get_surrounding_cell() const
+  {
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.get_properties();
+    return cell;
   }
 
 
 
   template <int dim, int spacedim>
-  inline typename Triangulation<dim, spacedim>::cell_iterator
+  inline const typename Triangulation<dim, spacedim>::cell_iterator &
   ParticleAccessor<dim, spacedim>::get_surrounding_cell(
-    const Triangulation<dim, spacedim> &triangulation) const
+    const Triangulation<dim, spacedim> & /*triangulation*/) const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    const typename Triangulation<dim, spacedim>::cell_iterator cell(
-      &triangulation, particle->first.first, particle->first.second);
     return cell;
   }
 
@@ -468,9 +530,9 @@ namespace Particles
   inline const ArrayView<double>
   ParticleAccessor<dim, spacedim>::get_properties()
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.get_properties();
+    return get_particle().get_properties();
   }
 
 
@@ -479,9 +541,9 @@ namespace Particles
   inline std::size_t
   ParticleAccessor<dim, spacedim>::serialized_size_in_bytes() const
   {
-    Assert(particle != map->end(), ExcInternalError());
+    Assert(state() == IteratorState::valid, ExcInternalError());
 
-    return particle->second.serialized_size_in_bytes();
+    return get_particle().serialized_size_in_bytes();
   }
 
 
@@ -490,8 +552,22 @@ namespace Particles
   inline void
   ParticleAccessor<dim, spacedim>::next()
   {
-    Assert(particle != map->end(), ExcInternalError());
-    ++particle;
+    Assert(state() == IteratorState::valid, ExcInternalError());
+
+    ++particle_index_within_cell;
+
+    if (particle_index_within_cell > (*particles)[active_cell_index].size() - 1)
+      {
+        particle_index_within_cell = 0;
+
+        do
+          {
+            ++cell;
+            ++active_cell_index;
+          }
+        while (cell.state() == IteratorState::valid &&
+               (*particles)[active_cell_index].size() == 0);
+      }
   }
 
 
@@ -500,8 +576,32 @@ namespace Particles
   inline void
   ParticleAccessor<dim, spacedim>::prev()
   {
-    Assert(particle != map->begin(), ExcInternalError());
-    --particle;
+    Assert(state() == IteratorState::valid, ExcInternalError());
+
+    if (particle_index_within_cell > 0)
+      --particle_index_within_cell;
+    else
+      {
+        do
+          {
+            --cell;
+            --active_cell_index;
+
+            if (cell.state() != IteratorState::valid)
+              {
+                active_cell_index = particles->size();
+                break;
+              }
+
+            if ((*particles)[active_cell_index].size() > 0)
+              {
+                particle_index_within_cell =
+                  (*particles)[active_cell_index].size() - 1;
+                break;
+              }
+          }
+        while ((*particles)[active_cell_index].size() == 0);
+      }
   }
 
 
@@ -511,7 +611,7 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::
   operator!=(const ParticleAccessor<dim, spacedim> &other) const
   {
-    return (map != other.map) || (particle != other.particle);
+    return !(*this == other);
   }
 
 
@@ -521,10 +621,47 @@ namespace Particles
   ParticleAccessor<dim, spacedim>::
   operator==(const ParticleAccessor<dim, spacedim> &other) const
   {
-    return (map == other.map) && (particle == other.particle);
+    return (particles == other.particles) && (cell == other.cell) &&
+           (particle_index_within_cell == other.particle_index_within_cell);
+  }
+
+
+
+  template <int dim, int spacedim>
+  inline IteratorState::IteratorStates
+  ParticleAccessor<dim, spacedim>::state() const
+  {
+    if (particles != nullptr && cell.state() == IteratorState::valid &&
+        particle_index_within_cell < (*particles)[active_cell_index].size())
+      return IteratorState::valid;
+    else if (particles != nullptr &&
+             cell.state() == IteratorState::past_the_end &&
+             particle_index_within_cell == 0)
+      return IteratorState::past_the_end;
+    else
+      return IteratorState::invalid;
+
+    return IteratorState::invalid;
+  }
+
+
+
+  template <int dim, int spacedim>
+  inline Particle<dim, spacedim> &
+  ParticleAccessor<dim, spacedim>::get_particle()
+  {
+    return (*particles)[active_cell_index][particle_index_within_cell];
   }
 
 
+
+  template <int dim, int spacedim>
+  inline const Particle<dim, spacedim> &
+  ParticleAccessor<dim, spacedim>::get_particle() const
+  {
+    return (*particles)[active_cell_index][particle_index_within_cell];
+  }
+
 } // namespace Particles
 
 DEAL_II_NAMESPACE_CLOSE
index e1954f513460fd23364413424aa87fe0caae7f9b..7a131922c3d661b2f24601a97653c2cfdffc3afe 100644 (file)
@@ -71,6 +71,12 @@ namespace Particles
      */
     using particle_iterator_range = boost::iterator_range<particle_iterator>;
 
+    /**
+     * A type for the storage container for particles.
+     */
+    using particle_container =
+      std::vector<std::vector<Particle<dim, spacedim>>>;
+
     /**
      * Default constructor.
      */
@@ -291,11 +297,22 @@ namespace Particles
       const;
 
     /**
-     * Remove a particle pointed to by the iterator.
+     * Remove a particle pointed to by the iterator. Afterwards the iterator
+     * will point to one of the remaining particles in the cell. Note that
+     * particle iterators that point to other particles in the same cell as @p particle
+     * may not longer be valid after this function call.
      */
     void
     remove_particle(const particle_iterator &particle);
 
+    /**
+     * Remove a vector of particles indicated by the particle iterators.
+     * The iterators and all other iterators are invalidated during the
+     * function call.
+     */
+    void
+    remove_particles(const std::vector<particle_iterator> &particles);
+
     /**
      * Insert a particle into the collection of particles. Return an iterator
      * to the new position of the particle. This function involves a copy of
@@ -833,14 +850,13 @@ namespace Particles
      * Set of particles currently living in the local domain, organized by
      * the level/index of the cell they are in.
      */
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>> particles;
+    particle_container particles;
 
     /**
-     * Set of particles that currently live in the ghost cells of the local
-     * domain, organized by the level/index of the cell they are in. These
-     * particles are equivalent to the ghost entries in distributed vectors.
+     * Set of particles currently living in the local domain, organized by
+     * the level/index of the cell they are in.
      */
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>> ghost_particles;
+    particle_container ghost_particles;
 
     /**
      * This variable stores how many particles are stored globally. It is
@@ -848,6 +864,12 @@ namespace Particles
      */
     types::particle_index global_number_of_particles;
 
+    /**
+     * This variable stores how many particles are stored locally. It is
+     * calculated by update_cached_numbers().
+     */
+    types::particle_index local_number_of_particles;
+
     /**
      * The maximum number of particles per cell in the global domain. This
      * variable is important to store and load particle data during
@@ -953,9 +975,8 @@ namespace Particles
     void
     send_recv_particles(
       const std::map<types::subdomain_id, std::vector<particle_iterator>>
-        &particles_to_send,
-      std::multimap<internal::LevelInd, Particle<dim, spacedim>>
-        &received_particles,
+        &                                                particles_to_send,
+      std::vector<std::vector<Particle<dim, spacedim>>> &received_particles,
       const std::map<
         types::subdomain_id,
         std::vector<
@@ -986,9 +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::multimap<internal::LevelInd, Particle<dim, spacedim>>
-        &received_particles);
+        &                                                particles_to_send,
+      std::vector<std::vector<Particle<dim, spacedim>>> &received_particles);
 
 
 #endif
@@ -1041,7 +1061,15 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::begin()
   {
-    return particle_iterator(particles, particles.begin());
+    if (particles.size() == 0)
+      return end();
+
+    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 end();
   }
 
 
@@ -1059,7 +1087,7 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end()
   {
-    return particle_iterator(particles, particles.end());
+    return particle_iterator(particles, triangulation->end(), 0);
   }
 
 
@@ -1077,7 +1105,15 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::begin_ghost()
   {
-    return particle_iterator(ghost_particles, ghost_particles.begin());
+    if (particles.size() == 0)
+      return end();
+
+    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 end_ghost();
   }
 
 
@@ -1095,7 +1131,7 @@ namespace Particles
   inline typename ParticleHandler<dim, spacedim>::particle_iterator
   ParticleHandler<dim, spacedim>::end_ghost()
   {
-    return particle_iterator(ghost_particles, ghost_particles.end());
+    return particle_iterator(ghost_particles, triangulation->end(), 0);
   }
 
 
index 9572e01a0bab328b06d1fe3f302dede1d0db82bc..b817880b25dfd1dd1aae23bdfdfbe920370ce7a2 100644 (file)
@@ -46,13 +46,13 @@ namespace Particles
 
     /**
      * Constructor of the iterator. Takes a reference to the particle
-     * container, and an iterator to the cell-particle pair.
+     * container, an iterator to the cell, and the particle index within that
+     * cell.
      */
     ParticleIterator(
-      const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
-      const typename std::multimap<internal::LevelInd,
-                                   Particle<dim, spacedim>>::iterator
-        &particle);
+      const std::vector<std::vector<Particle<dim, spacedim>>> &   particles,
+      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+      const unsigned int particle_index_within_cell);
 
     /**
      * Dereferencing operator, returns a reference to an accessor. Usage is thus
@@ -122,6 +122,12 @@ namespace Particles
     ParticleIterator
     operator--(int);
 
+    /**
+     * Return the state of the iterator.
+     */
+    IteratorState::IteratorStates
+    state() const;
+
     /**
      * Mark the class as bidirectional iterator and declare some alias which
      * are standard for iterators and are used by algorithms to enquire about
@@ -146,10 +152,10 @@ namespace Particles
 
   template <int dim, int spacedim>
   inline ParticleIterator<dim, spacedim>::ParticleIterator(
-    const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
-    const typename std::multimap<internal::LevelInd,
-                                 Particle<dim, spacedim>>::iterator & particle)
-    : accessor(map, particle)
+    const std::vector<std::vector<Particle<dim, spacedim>>> &   particles,
+    const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+    const unsigned int particle_index_within_cell)
+    : accessor(particles, cell, particle_index_within_cell)
   {}
 
 
@@ -253,6 +259,14 @@ namespace Particles
   }
 
 
+
+  template <int dim, int spacedim>
+  inline IteratorState::IteratorStates
+  ParticleIterator<dim, spacedim>::state() const
+  {
+    return accessor.state();
+  }
+
 } // namespace Particles
 
 DEAL_II_NAMESPACE_CLOSE
index 2a7ef41c5bf9b3ded19a040f599cf4c3ec853300..d17e7c0fc2777e51a9e946b2cc1dd23fbd5c1066 100644 (file)
@@ -94,9 +94,7 @@ namespace Particles
        * without clearing the multimap of ghost particles, thus greatly
        * reducing the cost of exchanging the ghost particles information.
        */
-      std::vector<typename std::multimap<internal::LevelInd,
-                                         Particle<dim, spacedim>>::iterator>
-        ghost_particles_iterators;
+      std::vector<particle_iterator> ghost_particles_iterators;
 
       /**
        * Temporary storage that holds the data of the particles to be sent
index 298a9021c790b5f46d99d009bfc3ac1c6f50de1c..ba099f690776036989f4cda142f2add2ebb716bb 100644 (file)
@@ -197,7 +197,6 @@ namespace Particles
           return; // nothing else to do here
         }
 
-      const auto &tria     = field_dh.get_triangulation();
       const auto &fe       = field_dh.get_fe();
       auto        particle = particle_handler.begin();
 
@@ -224,7 +223,7 @@ namespace Particles
 
       while (particle != particle_handler.end())
         {
-          const auto &cell = particle->get_surrounding_cell(tria);
+          const auto &cell = particle->get_surrounding_cell();
           const auto &dh_cell =
             typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &field_dh);
           dh_cell->get_dof_indices(dof_indices);
index d7dd259f8d83a4b68aeb0fbfd21612e8b9a4395c..d3aa57f0ed881bd2aad44773f74493d5457a2dda 100644 (file)
@@ -96,6 +96,7 @@ namespace Particles
     , particles()
     , ghost_particles()
     , global_number_of_particles(0)
+    , local_number_of_particles(0)
     , global_max_particles_per_cell(0)
     , next_free_particle_index(0)
     , size_callback()
@@ -114,9 +115,10 @@ namespace Particles
     : triangulation(&triangulation, typeid(*this).name())
     , mapping(&mapping, typeid(*this).name())
     , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
-    , particles()
-    , ghost_particles()
+    , particles(triangulation.n_active_cells())
+    , ghost_particles(triangulation.n_active_cells())
     , global_number_of_particles(0)
+    , local_number_of_particles(0)
     , global_max_particles_per_cell(0)
     , next_free_particle_index(0)
     , size_callback()
@@ -148,6 +150,9 @@ namespace Particles
     triangulation_cache =
       std::make_unique<GridTools::Cache<dim, spacedim>>(new_triangulation,
                                                         new_mapping);
+
+    particles.resize(triangulation->n_active_cells());
+    ghost_particles.resize(triangulation->n_active_cells());
   }
 
 
@@ -157,18 +162,20 @@ namespace Particles
   ParticleHandler<dim, spacedim>::copy_from(
     const ParticleHandler<dim, spacedim> &particle_handler)
   {
-    // clear and initialize this object before copying particles
+    // clear this object before copying particles
     clear();
+
     const unsigned int n_properties =
       particle_handler.property_pool->n_properties_per_slot();
     initialize(*particle_handler.triangulation,
                *particle_handler.mapping,
                n_properties);
-    property_pool->reserve(particle_handler.particles.size() +
-                           particle_handler.ghost_particles.size());
+
+    property_pool->reserve(particle_handler.n_locally_owned_particles());
 
     // copy static members
     global_number_of_particles = particle_handler.global_number_of_particles;
+    local_number_of_particles  = particle_handler.local_number_of_particles;
     global_max_particles_per_cell =
       particle_handler.global_max_particles_per_cell;
     next_free_particle_index = particle_handler.next_free_particle_index;
@@ -180,19 +187,11 @@ namespace Particles
     handle = particle_handler.handle;
 
     // copy dynamic properties
-    auto from_particle = particle_handler.begin();
     for (auto &particle : *this)
-      {
-        particle.set_property_pool(*property_pool);
-        ++from_particle;
-      }
+      particle.set_property_pool(*property_pool);
 
-    auto from_ghost = particle_handler.begin_ghost();
-    for (auto ghost = begin_ghost(); ghost != end_ghost();
-         ++ghost, ++from_ghost)
-      {
-        ghost->set_property_pool(*property_pool);
-      }
+    for (auto ghost = begin_ghost(); ghost != end_ghost(); ++ghost)
+      ghost->set_property_pool(*property_pool);
   }
 
 
@@ -203,6 +202,7 @@ namespace Particles
   {
     clear_particles();
     global_number_of_particles    = 0;
+    local_number_of_particles     = 0;
     next_free_particle_index      = 0;
     global_max_particles_per_cell = 0;
   }
@@ -214,7 +214,12 @@ namespace Particles
   ParticleHandler<dim, spacedim>::clear_particles()
   {
     particles.clear();
+    if (triangulation != nullptr)
+      particles.resize(triangulation->n_active_cells());
+
     ghost_particles.clear();
+    if (triangulation != nullptr)
+      ghost_particles.resize(triangulation->n_active_cells());
 
     // the particle properties have already been deleted by their destructor,
     // but the memory is still allocated. Return the memory as well.
@@ -227,26 +232,24 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::update_cached_numbers()
   {
-    types::particle_index locally_highest_index        = 0;
-    unsigned int          local_max_particles_per_cell = 0;
-    unsigned int          current_particles_per_cell   = 0;
-    typename Triangulation<dim, spacedim>::active_cell_iterator current_cell =
-      triangulation->begin_active();
+    local_number_of_particles                          = 0;
+    types::particle_index local_max_particle_index     = 0;
+    types::particle_index local_max_particles_per_cell = 0;
 
-    for (const auto &particle : *this)
+    for (const auto &particles_in_cell : particles)
       {
-        locally_highest_index =
-          std::max(locally_highest_index, particle.get_id());
+        const types::particle_index n_particles_in_cell =
+          particles_in_cell.size();
 
-        if (particle.get_surrounding_cell(*triangulation) != current_cell)
+        local_max_particles_per_cell =
+          std::max(local_max_particles_per_cell, n_particles_in_cell);
+        local_number_of_particles += n_particles_in_cell;
+
+        for (const auto &particle : particles_in_cell)
           {
-            current_particles_per_cell = 0;
-            current_cell = particle.get_surrounding_cell(*triangulation);
+            local_max_particle_index =
+              std::max(local_max_particle_index, particle.get_id());
           }
-
-        ++current_particles_per_cell;
-        local_max_particles_per_cell =
-          std::max(local_max_particles_per_cell, current_particles_per_cell);
       }
 
     if (const auto parallel_triangulation =
@@ -254,23 +257,33 @@ namespace Particles
             &*triangulation))
       {
         global_number_of_particles = dealii::Utilities::MPI::sum(
-          particles.size(), parallel_triangulation->get_communicator());
-        next_free_particle_index =
-          global_number_of_particles == 0 ?
-            0 :
-            dealii::Utilities::MPI::max(
-              locally_highest_index,
-              parallel_triangulation->get_communicator()) +
-              1;
-        global_max_particles_per_cell = dealii::Utilities::MPI::max(
-          local_max_particles_per_cell,
+          local_number_of_particles,
           parallel_triangulation->get_communicator());
+
+        if (global_number_of_particles == 0)
+          {
+            next_free_particle_index      = 0;
+            global_max_particles_per_cell = 0;
+          }
+        else
+          {
+            types::particle_index local_max_values[2] = {
+              local_max_particle_index, local_max_particles_per_cell};
+            types::particle_index global_max_values[2];
+
+            Utilities::MPI::max(local_max_values,
+                                parallel_triangulation->get_communicator(),
+                                global_max_values);
+
+            next_free_particle_index      = global_max_values[0] + 1;
+            global_max_particles_per_cell = global_max_values[1];
+          }
       }
     else
       {
-        global_number_of_particles = particles.size();
+        global_number_of_particles = local_number_of_particles;
         next_free_particle_index =
-          global_number_of_particles == 0 ? 0 : locally_highest_index + 1;
+          global_number_of_particles == 0 ? 0 : local_max_particle_index + 1;
         global_max_particles_per_cell = local_max_particles_per_cell;
       }
   }
@@ -283,13 +296,16 @@ namespace Particles
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
     const
   {
-    const internal::LevelInd found_cell =
-      std::make_pair(cell->level(), cell->index());
+    if (particles.size() == 0)
+      return 0;
 
-    if (cell->is_locally_owned())
-      return particles.count(found_cell);
-    else if (cell->is_ghost())
-      return ghost_particles.count(found_cell);
+    if (cell->is_artificial() == false)
+      {
+        if (cell->is_locally_owned() == true)
+          return particles[cell->active_cell_index()].size();
+        else
+          return ghost_particles[cell->active_cell_index()].size();
+      }
     else
       AssertThrow(false,
                   ExcMessage("You can't ask for the particles on an artificial "
@@ -318,22 +334,31 @@ namespace Particles
   ParticleHandler<dim, spacedim>::particles_in_cell(
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
   {
-    const internal::LevelInd level_index =
-      std::make_pair(cell->level(), cell->index());
+    const unsigned int active_cell_index = cell->active_cell_index();
 
-    if (cell->is_ghost())
+    if (cell->is_artificial() == false)
       {
-        const auto particles_in_cell = ghost_particles.equal_range(level_index);
-        return boost::make_iterator_range(
-          particle_iterator(ghost_particles, particles_in_cell.first),
-          particle_iterator(ghost_particles, particles_in_cell.second));
-      }
-    else if (cell->is_locally_owned())
-      {
-        const auto particles_in_cell = particles.equal_range(level_index);
-        return boost::make_iterator_range(
-          particle_iterator(particles, particles_in_cell.first),
-          particle_iterator(particles, particles_in_cell.second));
+        particle_container &container =
+          (cell->is_locally_owned() == true) ? particles : ghost_particles;
+
+        if (container[active_cell_index].size() == 0)
+          {
+            return boost::make_iterator_range(
+              particle_iterator(container, cell, 0),
+              particle_iterator(container, cell, 0));
+          }
+        else
+          {
+            particle_iterator begin(container, cell, 0);
+            particle_iterator end(container,
+                                  cell,
+                                  container[active_cell_index].size() - 1);
+            // end needs to point to the particle after the last one in the
+            // cell.
+            ++end;
+
+            return boost::make_iterator_range(begin, end);
+          }
       }
     else
       AssertThrow(false,
@@ -351,7 +376,70 @@ namespace Particles
   ParticleHandler<dim, spacedim>::remove_particle(
     const ParticleHandler<dim, spacedim>::particle_iterator &particle)
   {
-    particles.erase(particle->particle);
+    const unsigned int active_cell_index = particle->active_cell_index;
+
+    if (particles[active_cell_index].size() > 1)
+      {
+        particles[active_cell_index][particle->particle_index_within_cell] =
+          std::move(particles[active_cell_index].back());
+        particles[active_cell_index].resize(
+          particles[active_cell_index].size() - 1);
+      }
+    else
+      {
+        particles[active_cell_index].clear();
+      }
+
+    --local_number_of_particles;
+  }
+
+
+
+  template <int dim, int spacedim>
+  void
+  ParticleHandler<dim, spacedim>::remove_particles(
+    const std::vector<ParticleHandler<dim, spacedim>::particle_iterator>
+      &particles_to_remove)
+  {
+    unsigned int n_particles_removed = 0;
+
+    for (unsigned int cell_index = 0; cell_index < particles.size();
+         ++cell_index)
+      {
+        // If there is nothing left to remove, break and return
+        if (n_particles_removed == particles_to_remove.size())
+          break;
+
+        // Skip cells where there is nothing to remove
+        if (particles_to_remove[n_particles_removed]->active_cell_index !=
+            cell_index)
+          continue;
+
+        const unsigned int n_particles_in_cell = particles[cell_index].size();
+        unsigned int       move_to             = 0;
+        for (unsigned int move_from = 0; move_from < n_particles_in_cell;
+             ++move_from)
+          {
+            if (n_particles_removed != particles_to_remove.size() &&
+                particles_to_remove[n_particles_removed]->active_cell_index ==
+                  cell_index &&
+                particles_to_remove[n_particles_removed]
+                    ->particle_index_within_cell == move_from)
+              {
+                ++n_particles_removed;
+                continue;
+              }
+            else
+              {
+                particles[cell_index][move_to] =
+                  std::move(particles[cell_index][move_from]);
+                ++move_to;
+              }
+          }
+        particles[cell_index].resize(move_to);
+      }
+
+    update_cached_numbers();
   }
 
 
@@ -362,18 +450,20 @@ namespace Particles
     const Particle<dim, spacedim> &                                    particle,
     const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
   {
-    typename std::multimap<internal::LevelInd,
-                           Particle<dim, spacedim>>::iterator it =
-      particles.insert(
-        std::make_pair(internal::LevelInd(cell->level(), cell->index()),
-                       particle));
+    Assert(triangulation != nullptr, ExcInternalError());
+    Assert(cell.state() == IteratorState::valid, ExcInternalError());
 
-    particle_iterator particle_it(particles, it);
+    if (particles.size() == 0)
+      particles.resize(triangulation->n_active_cells());
+
+    const unsigned int active_cell_index = cell->active_cell_index();
+    particles[active_cell_index].push_back(particle);
+    particle_iterator particle_it(particles,
+                                  cell,
+                                  particles[active_cell_index].size() - 1);
     particle_it->set_property_pool(*property_pool);
 
-    if (particle.has_properties())
-      for (unsigned int n = 0; n < particle.get_properties().size(); ++n)
-        particle_it->get_properties()[n] = particle.get_properties()[n];
+    ++local_number_of_particles;
 
     return particle_it;
   }
@@ -387,16 +477,16 @@ namespace Particles
       typename Triangulation<dim, spacedim>::active_cell_iterator,
       Particle<dim, spacedim>> &new_particles)
   {
+    Assert(triangulation != nullptr, ExcInternalError());
+    if (particles.size() == 0)
+      particles.resize(triangulation->n_active_cells());
+
     for (const auto &particle : new_particles)
       {
-        // Insert the particle. Store an iterator to the newly
-        // inserted particle, and then set its property_pool.
-        auto it = particles.insert(
-          particles.end(),
-          std::make_pair(internal::LevelInd(particle.first->level(),
-                                            particle.first->index()),
-                         particle.second));
-        it->second.set_property_pool(*property_pool);
+        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);
       }
 
     update_cached_numbers();
@@ -409,6 +499,10 @@ namespace Particles
   ParticleHandler<dim, spacedim>::insert_particles(
     const std::vector<Point<spacedim>> &positions)
   {
+    Assert(triangulation != nullptr, ExcInternalError());
+    if (particles.size() == 0)
+      particles.resize(triangulation->n_active_cells());
+
     update_cached_numbers();
 
     // Determine the starting particle index of this process, which
@@ -456,22 +550,19 @@ namespace Particles
     if (cells.size() == 0)
       return;
 
-    auto hint =
-      particles.find(std::make_pair(cells[0]->level(), cells[0]->index()));
     for (unsigned int i = 0; i < cells.size(); ++i)
       {
-        internal::LevelInd current_cell(cells[i]->level(), cells[i]->index());
+        const unsigned int active_cell_index = cells[i]->active_cell_index();
+
         for (unsigned int p = 0; p < local_positions[i].size(); ++p)
           {
-            hint = particles.insert(
-              hint,
-              std::make_pair(current_cell,
-                             Particle<dim, spacedim>(positions[index_map[i][p]],
-                                                     local_positions[i][p],
-                                                     local_start_index +
-                                                       index_map[i][p])));
-
-            hint->second.set_property_pool(*property_pool);
+            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]));
+
+            particles[active_cell_index].back().set_property_pool(
+              *property_pool);
           }
       }
 
@@ -779,7 +870,7 @@ namespace Particles
   types::particle_index
   ParticleHandler<dim, spacedim>::n_locally_owned_particles() const
   {
-    return particles.size();
+    return local_number_of_particles;
   }
 
 
@@ -939,33 +1030,47 @@ namespace Particles
     // processes around (with an invalid cell).
 
     std::vector<particle_iterator> particles_out_of_cell;
-    particles_out_of_cell.reserve(n_locally_owned_particles());
+
+    particles_out_of_cell.reserve(n_locally_owned_particles() / 4);
 
     // Now update the reference locations of the moved particles
     std::vector<Point<spacedim>> real_locations;
     std::vector<Point<dim>>      reference_locations;
-    for (auto particle = begin(); particle != end();)
+    real_locations.reserve(global_max_particles_per_cell);
+    reference_locations.reserve(global_max_particles_per_cell);
+
+    for (const auto &cell : triangulation->active_cell_iterators())
       {
-        const auto cell = particle->get_surrounding_cell(*triangulation);
+        const unsigned int active_cell_index = cell->active_cell_index();
+        const unsigned int n_pic = particles[active_cell_index].size();
+
+        // Particles can be inserted into arbitrary cells, e.g. if their cell is
+        // not known. However, for artificial cells we can not check evaluate
+        // the reference position of particles and particles in ghost cells
+        // should not be owned by this process. Assume all particles that are in
+        // not locally owned cells have to be resorted or transferred.
+        if (cell->is_locally_owned() == false)
+          {
+            for (unsigned int i = 0; i < n_pic; ++i)
+              {
+                particles_out_of_cell.push_back(
+                  particle_iterator(particles, cell, i));
+              }
+            continue;
+          }
+
+        auto pic = particles_in_cell(cell);
+
         real_locations.clear();
+        for (const auto &particle : pic)
+          real_locations.push_back(particle.get_location());
 
-        // Since we might also work on artificial cells when we initialize the
-        // particles on a remote processor, we cannot use the
-        // particles_in_cell method. Thus, We instead simply go through the
-        // particles and check if the next one belongs to the same cell as the
-        // current one.
-        for (auto it = particle;
-             it != end() && it->get_surrounding_cell(*triangulation) == cell;
-             ++it)
-          real_locations.push_back(it->get_location());
-
-        reference_locations.resize(real_locations.size());
-        ArrayView<Point<dim>> reference(reference_locations.data(),
-                                        reference_locations.size());
+        reference_locations.resize(n_pic);
         mapping->transform_points_real_to_unit_cell(cell,
                                                     real_locations,
-                                                    reference);
+                                                    reference_locations);
 
+        auto particle = pic.begin();
         for (const auto &p_unit : reference_locations)
           {
             if (p_unit[0] == std::numeric_limits<double>::infinity() ||
@@ -973,18 +1078,17 @@ namespace Particles
               particles_out_of_cell.push_back(particle);
             else
               particle->set_reference_location(p_unit);
+
             ++particle;
           }
       }
 
     // There are three reasons why a particle is not in its old cell:
     // It moved to another cell, to another subdomain or it left the mesh.
-    // Particles that moved to another cell are updated and stored inside the
-    // sorted_particles vector, particles that moved to another domain are
+    // Particles that moved to another cell are updated and moved inside the
+    // particles vector, particles that moved to another domain are
     // collected in the moved_particles_domain vector. Particles that left
     // the mesh completely are ignored and removed.
-    std::vector<std::pair<internal::LevelInd, Particle<dim, spacedim>>>
-      sorted_particles;
     std::map<types::subdomain_id, std::vector<particle_iterator>>
       moved_particles;
     std::map<
@@ -998,8 +1102,6 @@ namespace Particles
     // automatic and relatively fast (compared to other parts of this
     // algorithm) re-allocation will happen.
     using vector_size = typename std::vector<particle_iterator>::size_type;
-    sorted_particles.reserve(
-      static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
 
     std::set<types::subdomain_id> ghost_owners;
     if (const auto parallel_triangulation =
@@ -1007,38 +1109,41 @@ namespace Particles
             &*triangulation))
       ghost_owners = parallel_triangulation->ghost_owners();
 
-    for (const auto ghost_owner : ghost_owners)
+    for (const auto &ghost_owner : ghost_owners)
       moved_particles[ghost_owner].reserve(
         static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
-    for (const auto ghost_owner : ghost_owners)
+    for (const auto &ghost_owner : ghost_owners)
       moved_cells[ghost_owner].reserve(
         static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
 
     {
       // Create a map from vertices to adjacent cells using grid cache
-      std::vector<
+      const std::vector<
         std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
-        vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
+        &vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
 
-      // Create a corresponding map of vectors from vertex to cell center using
-      // grid cache
-      std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers =
-        triangulation_cache->get_vertex_to_cell_centers_directions();
+      // Create a corresponding map of vectors from vertex to cell center
+      // using grid cache
+      const std::vector<std::vector<Tensor<1, spacedim>>>
+        &vertex_to_cell_centers =
+          triangulation_cache->get_vertex_to_cell_centers_directions();
 
       std::vector<unsigned int> neighbor_permutation;
 
       // Find the cells that the particles moved to.
       for (auto &out_particle : particles_out_of_cell)
         {
+          // make a copy of the current cell, since we will modify the
+          // variable current_cell in the following but we need the original in
+          // the case the particle is not found
+          auto current_cell = out_particle->get_surrounding_cell();
+
           // The cell the particle is in
           Point<dim> current_reference_position;
           bool       found_cell = false;
 
           // Check if the particle is in one of the old cell's neighbors
           // that are adjacent to the closest vertex
-          typename Triangulation<dim, spacedim>::active_cell_iterator
-            current_cell = out_particle->get_surrounding_cell(*triangulation);
-
           const unsigned int closest_vertex =
             GridTools::find_closest_vertex_of_cell<dim, spacedim>(
               current_cell, out_particle->get_location(), *mapping);
@@ -1055,7 +1160,7 @@ namespace Particles
           for (unsigned int i = 0; i < n_neighbor_cells; ++i)
             neighbor_permutation[i] = i;
 
-          const auto cell_centers =
+          const auto &cell_centers =
             vertex_to_cell_centers[closest_vertex_index];
           std::sort(neighbor_permutation.begin(),
                     neighbor_permutation.end(),
@@ -1067,11 +1172,6 @@ namespace Particles
                                                           cell_centers);
                     });
 
-          // make a copy of the current cell, since we will modify the variable
-          // current_cell in the following but we need a backup in the case
-          // the particle is not found
-          const auto previous_cell_of_particle = current_cell;
-
           // Search all of the cells adjacent to the closest vertex of the
           // previous cell Most likely we will find the particle in them.
           for (unsigned int i = 0; i < n_neighbor_cells; ++i)
@@ -1118,7 +1218,7 @@ namespace Particles
                   // domain due to an integration error or an open boundary.
                   // Signal the loss and move on.
                   signals.particle_lost(out_particle,
-                                        previous_cell_of_particle);
+                                        out_particle->get_surrounding_cell());
                   continue;
                 }
             }
@@ -1131,10 +1231,11 @@ namespace Particles
           // Mark it for MPI transfer otherwise
           if (current_cell->is_locally_owned())
             {
-              sorted_particles.push_back(
-                std::make_pair(internal::LevelInd(current_cell->level(),
-                                                  current_cell->index()),
-                               out_particle->particle->second));
+              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]));
             }
           else
             {
@@ -1145,11 +1246,6 @@ namespace Particles
         }
     }
 
-    // Sort the updated particles. This pre-sort speeds up inserting
-    // them into particles to O(N) complexity.
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>>
-      sorted_particles_map;
-
     // Exchange particles between processors if we have more than one process
 #ifdef DEAL_II_WITH_MPI
     if (const auto parallel_triangulation =
@@ -1158,21 +1254,13 @@ namespace Particles
       {
         if (dealii::Utilities::MPI::n_mpi_processes(
               parallel_triangulation->get_communicator()) > 1)
-          send_recv_particles(moved_particles,
-                              sorted_particles_map,
-                              moved_cells);
+          send_recv_particles(moved_particles, particles, moved_cells);
       }
 #endif
 
-    sorted_particles_map.insert(sorted_particles.begin(),
-                                sorted_particles.end());
-
-    for (unsigned int i = 0; i < particles_out_of_cell.size(); ++i)
-      remove_particle(particles_out_of_cell[i]);
-
-    particles.insert(sorted_particles_map.begin(), sorted_particles_map.end());
-    update_cached_numbers();
-  }
+    // remove_particles also calls update_cached_numbers()
+    remove_particles(particles_out_of_cell);
+  } // namespace Particles
 
 
 
@@ -1199,6 +1287,7 @@ namespace Particles
 #else
     // First clear the current ghost_particle information
     ghost_particles.clear();
+    ghost_particles.resize(triangulation->n_active_cells());
 
     // Clear ghost particles data structures and invalidate cache
     ghost_particles_cache.ghost_particles_by_domain.clear();
@@ -1210,7 +1299,7 @@ namespace Particles
     for (const auto ghost_owner : ghost_owners)
       ghost_particles_cache.ghost_particles_by_domain[ghost_owner].reserve(
         static_cast<typename std::vector<particle_iterator>::size_type>(
-          particles.size() * 0.25));
+          local_number_of_particles * 0.25));
 
     const std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain =
       triangulation_cache->get_vertex_to_neighbor_subdomain();
@@ -1275,10 +1364,10 @@ namespace Particles
 #ifdef DEAL_II_WITH_MPI
     // First clear the current ghost_particle information
     // ghost_particles.clear();
-    Assert(
-      ghost_particles_cache.valid,
-      ExcMessage(
-        "Ghost particles cannot be updated if they first have not been exchanged at least once with the cache enabled"));
+    Assert(ghost_particles_cache.valid,
+           ExcMessage(
+             "Ghost particles cannot be updated if they first have not been "
+             "exchanged at least once with the cache enabled"));
 
 
     send_recv_particles_properties_and_location(
@@ -1293,9 +1382,8 @@ namespace Particles
   void
   ParticleHandler<dim, spacedim>::send_recv_particles(
     const std::map<types::subdomain_id, std::vector<particle_iterator>>
-      &particles_to_send,
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>>
-      &received_particles,
+      &                 particles_to_send,
+    particle_container &received_particles,
     const std::map<
       types::subdomain_id,
       std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
@@ -1310,7 +1398,11 @@ namespace Particles
     Assert(
       parallel_triangulation,
       ExcMessage(
-        "This function is only implemented for parallel::TriangulationBase objects."));
+        "This function is only implemented for parallel::TriangulationBase "
+        "objects."));
+
+    if (received_particles.size() == 0)
+      received_particles.resize(triangulation->n_active_cells());
 
     // Determine the communication pattern
     const std::set<types::subdomain_id> ghost_owners =
@@ -1350,12 +1442,12 @@ namespace Particles
     std::vector<unsigned int> send_offsets(n_neighbors, 0);
     std::vector<char>         send_data;
 
+    Particle<dim, spacedim> test_particle;
+    test_particle.set_property_pool(*property_pool);
+
     const unsigned int individual_particle_data_size =
-      Utilities::MPI::max(n_send_particles > 0 ?
-                            ((begin()->serialized_size_in_bytes() +
-                              (size_callback ? size_callback() : 0))) :
-                            0,
-                          parallel_triangulation->get_communicator());
+      test_particle.serialized_size_in_bytes() +
+      (size_callback ? size_callback() : 0);
 
     const unsigned int individual_total_particle_data_size =
       individual_particle_data_size + cellid_size;
@@ -1385,17 +1477,18 @@ namespace Particles
                      static_cast<std::size_t>(
                        n_particles_to_send *
                        individual_total_particle_data_size),
-                   ExcMessage("Overflow when trying to send particle data"));
+                   ExcMessage("Overflow when trying to send particle "
+                              "data"));
 
             for (unsigned int j = 0; j < n_particles_to_send; ++j)
               {
-                // If no target cells are given, use the iterator information
+                // If no target cells are given, use the iterator
+                // information
                 typename Triangulation<dim, spacedim>::active_cell_iterator
                   cell;
                 if (send_cells.size() == 0)
-                  cell =
-                    particles_to_send.at(neighbors[i])[j]->get_surrounding_cell(
-                      *triangulation);
+                  cell = particles_to_send.at(neighbors[i])[j]
+                           ->get_surrounding_cell();
                 else
                   cell = send_cells.at(neighbors[i])[j];
 
@@ -1553,19 +1646,21 @@ namespace Particles
         const typename Triangulation<dim, spacedim>::active_cell_iterator cell =
           triangulation->create_cell_iterator(id);
 
-        typename std::multimap<internal::LevelInd,
-                               Particle<dim, spacedim>>::iterator
-          recv_particle = received_particles.insert(std::make_pair(
-            internal::LevelInd(cell->level(), cell->index()),
-            Particle<dim, spacedim>(recv_data_it, property_pool.get())));
+        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()));
+
+        particle_iterator particle_it(
+          received_particles,
+          cell,
+          received_particles[active_cell_index].size() - 1);
+        particle_it->set_property_pool(*property_pool);
 
         if (load_callback)
-          recv_data_it =
-            load_callback(particle_iterator(received_particles, recv_particle),
-                          recv_data_it);
+          recv_data_it = load_callback(particle_it, recv_data_it);
 
         if (build_cache) // TODO: is this safe?
-          ghost_particles_iterators.push_back(recv_particle);
+          ghost_particles_iterators.push_back(particle_it);
       }
 
     AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
@@ -1582,9 +1677,8 @@ 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,
-    std::multimap<internal::LevelInd, Particle<dim, spacedim>>
-      &updated_particles)
+      &                 particles_to_send,
+    particle_container &updated_particles)
   {
     const auto &neighbors     = ghost_particles_cache.neighbors;
     const auto &send_pointers = ghost_particles_cache.send_pointers;
@@ -1596,7 +1690,8 @@ namespace Particles
     Assert(
       parallel_triangulation,
       ExcMessage(
-        "This function is only implemented for parallel::TriangulationBase objects."));
+        "This function is only implemented for parallel::TriangulationBase "
+        "objects."));
 
     std::vector<char> &send_data = ghost_particles_cache.send_data;
 
@@ -1673,12 +1768,14 @@ namespace Particles
         // Update particle data using previously allocated memory space
         // for efficiency reasons
         recv_data_it =
-          recv_particle->second.read_particle_data_from_memory(recv_data_it);
+          recv_particle->read_particle_data_from_memory(recv_data_it);
 
         if (load_callback)
-          recv_data_it =
-            load_callback(particle_iterator(updated_particles, recv_particle),
-                          recv_data_it);
+          recv_data_it = load_callback(
+            particle_iterator(updated_particles,
+                              recv_particle->cell,
+                              recv_particle->particle_index_within_cell),
+            recv_data_it);
       }
 
     AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
@@ -1820,27 +1917,7 @@ namespace Particles
           // If the cell persist or is refined store all particles of the
           // current cell.
           {
-            unsigned int n_particles = 0;
-
-            const internal::LevelInd level_index = {cell->level(),
-                                                    cell->index()};
-            const auto               particles_in_cell =
-              (cell->is_ghost() ? ghost_particles.equal_range(level_index) :
-                                  particles.equal_range(level_index));
-
-            n_particles = n_particles_in_cell(cell);
-            stored_particles_on_cell.reserve(n_particles);
-
-            std::for_each(
-              particles_in_cell.first,
-              particles_in_cell.second,
-              [&stored_particles_on_cell](
-                const std::pair<internal::LevelInd, Particle<dim, spacedim>>
-                  &particle) {
-                stored_particles_on_cell.push_back(particle.second);
-              });
-
-            AssertDimension(n_particles, stored_particles_on_cell.size());
+            stored_particles_on_cell = particles[cell->active_cell_index()];
           }
           break;
 
@@ -1848,35 +1925,19 @@ namespace Particles
           // If this cell is the parent of children that will be coarsened,
           // collect the particles of all children.
           {
-            unsigned int n_particles = 0;
-
             for (const auto &child : cell->child_iterators())
               {
-                n_particles += n_particles_in_cell(child);
-              }
+                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(n_particles);
+                stored_particles_on_cell.reserve(
+                  stored_particles_on_cell.size() + n_particles);
 
-            for (const auto &child : cell->child_iterators())
-              {
-                const internal::LevelInd level_index = {child->level(),
-                                                        child->index()};
-                const auto               particles_in_cell =
-                  (child->is_ghost() ?
-                     ghost_particles.equal_range(level_index) :
-                     particles.equal_range(level_index));
-
-                std::for_each(
-                  particles_in_cell.first,
-                  particles_in_cell.second,
-                  [&stored_particles_on_cell](
-                    const std::pair<internal::LevelInd, Particle<dim, spacedim>>
-                      &particle) {
-                    stored_particles_on_cell.push_back(particle.second);
-                  });
+                for (unsigned int i = 0; i < n_particles; ++i)
+                  stored_particles_on_cell.push_back(
+                    particles[active_cell_index][i]);
               }
-
-            AssertDimension(n_particles, stored_particles_on_cell.size());
           }
           break;
 
@@ -1912,66 +1973,30 @@ namespace Particles
       {
         case parallel::TriangulationBase<dim, spacedim>::CELL_PERSIST:
           {
-            auto position_hint = particles.end();
-            for (const auto &particle : loaded_particles_on_cell)
-              {
-                // Use std::multimap::emplace_hint to speed up insertion of
-                // particles.
-                position_hint =
-                  particles.emplace_hint(position_hint,
-                                         std::make_pair(cell->level(),
-                                                        cell->index()),
-                                         std::move(particle));
-                // Move the hint position forward by one, i.e., for the next
-                // particle. The 'hint' position will thus be right after the
-                // one just inserted.
-                ++position_hint;
-              }
+            particles[cell->active_cell_index()].insert(
+              particles[cell->active_cell_index()].end(),
+              loaded_particles_on_cell.begin(),
+              loaded_particles_on_cell.end());
           }
           break;
 
         case parallel::TriangulationBase<dim, spacedim>::CELL_COARSEN:
           {
-            typename std::multimap<internal::LevelInd,
-                                   Particle<dim, spacedim>>::iterator
-              position_hint = particles.end();
+            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);
-                // Use std::multimap::emplace_hint to speed up insertion of
-                // particles.
-                position_hint =
-                  particles.emplace_hint(position_hint,
-                                         std::make_pair(cell->level(),
-                                                        cell->index()),
-                                         std::move(particle));
-                // Move the hint position forward by one, i.e., for the next
-                // particle. The 'hint' position will thus be right after the
-                // one just inserted.
-                ++position_hint;
+                particles[active_cell_index].emplace_back(std::move(particle));
               }
           }
           break;
 
         case parallel::TriangulationBase<dim, spacedim>::CELL_REFINE:
           {
-            std::vector<
-              typename std::multimap<internal::LevelInd,
-                                     Particle<dim, spacedim>>::iterator>
-              position_hints(GeometryInfo<dim>::max_children_per_cell);
-            for (unsigned int child_index = 0;
-                 child_index < GeometryInfo<dim>::max_children_per_cell;
-                 ++child_index)
-              {
-                const typename Triangulation<dim, spacedim>::cell_iterator
-                  child                     = cell->child(child_index);
-                position_hints[child_index] = particles.upper_bound(
-                  std::make_pair(child->level(), child->index()));
-              }
-
             for (auto &particle : loaded_particles_on_cell)
               {
                 for (unsigned int child_index = 0;
@@ -1979,7 +2004,9 @@ namespace Particles
                      ++child_index)
                   {
                     const typename Triangulation<dim, spacedim>::cell_iterator
-                      child = cell->child(child_index);
+                                       child = cell->child(child_index);
+                    const unsigned int active_cell_index =
+                      child->active_cell_index();
 
                     try
                       {
@@ -1991,15 +2018,8 @@ namespace Particles
                             particle.set_reference_location(p_unit);
                             // Use std::multimap::emplace_hint to speed up
                             // insertion of particles.
-                            position_hints[child_index] =
-                              particles.emplace_hint(
-                                position_hints[child_index],
-                                std::make_pair(child->level(), child->index()),
-                                std::move(particle));
-                            // Move the hint position forward by one, i.e.,
-                            // for the next particle. The 'hint' position will
-                            // thus be right after the one just inserted.
-                            ++position_hints[child_index];
+                            particles[active_cell_index].emplace_back(
+                              std::move(particle));
                             break;
                           }
                       }
index 30a4b8100169ea5ce4318c4e02338a012f0893ba..897abf99b951594dad63af1fde3ad5f2acb243ee 100644 (file)
@@ -39,8 +39,7 @@ namespace Particles
       if (particle_handler.n_locally_owned_particles() == 0)
         return; // nothing to do here
 
-      const auto &tria = space_dh.get_triangulation();
-      const auto &fe   = space_dh.get_fe();
+      const auto &fe = space_dh.get_fe();
       const auto  max_particles_per_cell =
         particle_handler.n_global_max_particles_per_cell();
 
@@ -80,7 +79,7 @@ namespace Particles
       auto particle = particle_handler.begin();
       while (particle != particle_handler.end())
         {
-          const auto &cell = particle->get_surrounding_cell(tria);
+          const auto &cell = particle->get_surrounding_cell();
           const auto  dh_cell =
             typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
           dh_cell->get_dof_indices(dof_indices);
@@ -127,8 +126,7 @@ namespace Particles
 
       AssertDimension(matrix.n(), space_dh.n_dofs());
 
-      const auto &tria = space_dh.get_triangulation();
-      const auto &fe   = space_dh.get_fe();
+      const auto &fe = space_dh.get_fe();
       const auto  max_particles_per_cell =
         particle_handler.n_global_max_particles_per_cell();
 
@@ -174,7 +172,7 @@ namespace Particles
       auto particle = particle_handler.begin();
       while (particle != particle_handler.end())
         {
-          const auto &cell = particle->get_surrounding_cell(tria);
+          const auto &cell = particle->get_surrounding_cell();
           const auto &dh_cell =
             typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
           dh_cell->get_dof_indices(dof_indices);
index 510a4e0f6189f33d3109e4256b45063e16506ba3..7c5dccce3d6e0bb955ce7671b6e97b21a120dcc5 100644 (file)
@@ -70,8 +70,11 @@ test()
                                                                       2,
                                                                       0);
 
-    auto particle_it = particle_handler.insert_particle(particle1, cell1);
-    particle_it->set_properties(properties);
+    if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
+      {
+        auto particle_it = particle_handler.insert_particle(particle1, cell1);
+        particle_it->set_properties(properties);
+      }
 
     std::vector<std::string> data_names;
     std::vector<DataComponentInterpretation::DataComponentInterpretation>
@@ -101,7 +104,7 @@ main(int argc, char *argv[])
 {
   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
 
-  initlog();
+  MPILogInitAll all;
   deallog.push("2d/2d");
   test<2, 2>();
   deallog.pop();
index 99cd3890803406854ac942450cc9ff45ac485f39..66fea8ec361e21cb99819d5ff2956b6bc1857a44 100644 (file)
@@ -8,7 +8,7 @@
 # <x> <y> <id> <property_coord_0> <property_coord_1> 
 0.125000 0.125000 0.00000 0.125000 0.125000 
 
-DEAL:2d/2d::OK
+DEAL:0:2d/2d::OK
 # This file was generated by the deal.II library.
 
 
@@ -18,7 +18,7 @@ DEAL:2d/2d::OK
 # <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2> 
 0.125000 0.125000 0.00000 0.00000 0.125000 0.125000 0.125000 
 
-DEAL:2d/3d::OK
+DEAL:0:2d/3d::OK
 # This file was generated by the deal.II library.
 
 
@@ -28,4 +28,9 @@ DEAL:2d/3d::OK
 # <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2> 
 0.125000 0.125000 0.125000 0.00000 0.125000 0.125000 0.125000 
 
-DEAL:3d/3d::OK
+DEAL:0:3d/3d::OK
+
+DEAL:1:2d/2d::OK
+DEAL:1:2d/3d::OK
+DEAL:1:3d/3d::OK
+
index cc99df90b7c2644e7042c5e9b2568e4484b77ed3..83fa81748f852eb858c2f31c903694f4c0b12f28 100644 (file)
@@ -39,8 +39,8 @@ test()
     tr.refine_global(2);
     MappingQ<dim, spacedim> mapping(1);
 
-    // both processes create a particle handler, but only the first creates
-    // particles
+    // both processes create a particle handler with a particle in a
+    // position that is in a ghost cell to the other process
     Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
 
     Point<spacedim> position;
@@ -57,6 +57,9 @@ test()
       position,
       reference_position,
       Utilities::MPI::this_mpi_process(tr.get_communicator()));
+
+    // We give a random cell hint to check that sorting and
+    // transferring ghost particles works.
     typename Triangulation<dim, spacedim>::active_cell_iterator cell =
       tr.begin_active();
     particle_handler.insert_particle(particle, cell);
index 1c9486191a057c83e38f74ceb6d73b904aadda9d..6b4962720aefd4aa96853757dba88f2cfbb17796 100644 (file)
@@ -1,65 +1,65 @@
 
 DEAL:2d/2d::Particle number: 10
-DEAL:2d/2d::Particle id 2 is in cell 2.0
-DEAL:2d/2d::     at location 0.250000 0.250000
-DEAL:2d/2d::     at reference location 1.00000 1.00000
-DEAL:2d/2d::Particle id 1 is in cell 2.0
-DEAL:2d/2d::     at location 0.150000 0.150000
-DEAL:2d/2d::     at reference location 0.600000 0.600000
 DEAL:2d/2d::Particle id 0 is in cell 2.0
 DEAL:2d/2d::     at location 0.0500000 0.0500000
 DEAL:2d/2d::     at reference location 0.200000 0.200000
-DEAL:2d/2d::Particle id 4 is in cell 2.3
-DEAL:2d/2d::     at location 0.450000 0.450000
-DEAL:2d/2d::     at reference location 0.800000 0.800000
+DEAL:2d/2d::Particle id 1 is in cell 2.0
+DEAL:2d/2d::     at location 0.150000 0.150000
+DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 2 is in cell 2.0
+DEAL:2d/2d::     at location 0.250000 0.250000
+DEAL:2d/2d::     at reference location 1.00000 1.00000
 DEAL:2d/2d::Particle id 3 is in cell 2.3
 DEAL:2d/2d::     at location 0.350000 0.350000
 DEAL:2d/2d::     at reference location 0.400000 0.400000
-DEAL:2d/2d::Particle id 7 is in cell 2.12
-DEAL:2d/2d::     at location 0.750000 0.750000
-DEAL:2d/2d::     at reference location 1.00000 1.00000
-DEAL:2d/2d::Particle id 6 is in cell 2.12
-DEAL:2d/2d::     at location 0.650000 0.650000
-DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 4 is in cell 2.3
+DEAL:2d/2d::     at location 0.450000 0.450000
+DEAL:2d/2d::     at reference location 0.800000 0.800000
 DEAL:2d/2d::Particle id 5 is in cell 2.12
 DEAL:2d/2d::     at location 0.550000 0.550000
 DEAL:2d/2d::     at reference location 0.200000 0.200000
-DEAL:2d/2d::Particle id 9 is in cell 2.15
-DEAL:2d/2d::     at location 0.950000 0.950000
-DEAL:2d/2d::     at reference location 0.800000 0.800000
+DEAL:2d/2d::Particle id 6 is in cell 2.12
+DEAL:2d/2d::     at location 0.650000 0.650000
+DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 7 is in cell 2.12
+DEAL:2d/2d::     at location 0.750000 0.750000
+DEAL:2d/2d::     at reference location 1.00000 1.00000
 DEAL:2d/2d::Particle id 8 is in cell 2.15
 DEAL:2d/2d::     at location 0.850000 0.850000
 DEAL:2d/2d::     at reference location 0.400000 0.400000
+DEAL:2d/2d::Particle id 9 is in cell 2.15
+DEAL:2d/2d::     at location 0.950000 0.950000
+DEAL:2d/2d::     at reference location 0.800000 0.800000
 DEAL:2d/2d::OK
 DEAL:3d/3d::Particle number: 10
-DEAL:3d/3d::Particle id 2 is in cell 2.0
-DEAL:3d/3d::     at location 0.250000 0.250000 0.250000
-DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
-DEAL:3d/3d::Particle id 1 is in cell 2.0
-DEAL:3d/3d::     at location 0.150000 0.150000 0.150000
-DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
 DEAL:3d/3d::Particle id 0 is in cell 2.0
 DEAL:3d/3d::     at location 0.0500000 0.0500000 0.0500000
 DEAL:3d/3d::     at reference location 0.200000 0.200000 0.200000
-DEAL:3d/3d::Particle id 4 is in cell 2.7
-DEAL:3d/3d::     at location 0.450000 0.450000 0.450000
-DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
+DEAL:3d/3d::Particle id 1 is in cell 2.0
+DEAL:3d/3d::     at location 0.150000 0.150000 0.150000
+DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 2 is in cell 2.0
+DEAL:3d/3d::     at location 0.250000 0.250000 0.250000
+DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
 DEAL:3d/3d::Particle id 3 is in cell 2.7
 DEAL:3d/3d::     at location 0.350000 0.350000 0.350000
 DEAL:3d/3d::     at reference location 0.400000 0.400000 0.400000
-DEAL:3d/3d::Particle id 7 is in cell 2.56
-DEAL:3d/3d::     at location 0.750000 0.750000 0.750000
-DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
-DEAL:3d/3d::Particle id 6 is in cell 2.56
-DEAL:3d/3d::     at location 0.650000 0.650000 0.650000
-DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 4 is in cell 2.7
+DEAL:3d/3d::     at location 0.450000 0.450000 0.450000
+DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
 DEAL:3d/3d::Particle id 5 is in cell 2.56
 DEAL:3d/3d::     at location 0.550000 0.550000 0.550000
 DEAL:3d/3d::     at reference location 0.200000 0.200000 0.200000
-DEAL:3d/3d::Particle id 9 is in cell 2.63
-DEAL:3d/3d::     at location 0.950000 0.950000 0.950000
-DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
+DEAL:3d/3d::Particle id 6 is in cell 2.56
+DEAL:3d/3d::     at location 0.650000 0.650000 0.650000
+DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 7 is in cell 2.56
+DEAL:3d/3d::     at location 0.750000 0.750000 0.750000
+DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
 DEAL:3d/3d::Particle id 8 is in cell 2.63
 DEAL:3d/3d::     at location 0.850000 0.850000 0.850000
 DEAL:3d/3d::     at reference location 0.400000 0.400000 0.400000
+DEAL:3d/3d::Particle id 9 is in cell 2.63
+DEAL:3d/3d::     at location 0.950000 0.950000 0.950000
+DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
 DEAL:3d/3d::OK
index 1c9486191a057c83e38f74ceb6d73b904aadda9d..6b4962720aefd4aa96853757dba88f2cfbb17796 100644 (file)
@@ -1,65 +1,65 @@
 
 DEAL:2d/2d::Particle number: 10
-DEAL:2d/2d::Particle id 2 is in cell 2.0
-DEAL:2d/2d::     at location 0.250000 0.250000
-DEAL:2d/2d::     at reference location 1.00000 1.00000
-DEAL:2d/2d::Particle id 1 is in cell 2.0
-DEAL:2d/2d::     at location 0.150000 0.150000
-DEAL:2d/2d::     at reference location 0.600000 0.600000
 DEAL:2d/2d::Particle id 0 is in cell 2.0
 DEAL:2d/2d::     at location 0.0500000 0.0500000
 DEAL:2d/2d::     at reference location 0.200000 0.200000
-DEAL:2d/2d::Particle id 4 is in cell 2.3
-DEAL:2d/2d::     at location 0.450000 0.450000
-DEAL:2d/2d::     at reference location 0.800000 0.800000
+DEAL:2d/2d::Particle id 1 is in cell 2.0
+DEAL:2d/2d::     at location 0.150000 0.150000
+DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 2 is in cell 2.0
+DEAL:2d/2d::     at location 0.250000 0.250000
+DEAL:2d/2d::     at reference location 1.00000 1.00000
 DEAL:2d/2d::Particle id 3 is in cell 2.3
 DEAL:2d/2d::     at location 0.350000 0.350000
 DEAL:2d/2d::     at reference location 0.400000 0.400000
-DEAL:2d/2d::Particle id 7 is in cell 2.12
-DEAL:2d/2d::     at location 0.750000 0.750000
-DEAL:2d/2d::     at reference location 1.00000 1.00000
-DEAL:2d/2d::Particle id 6 is in cell 2.12
-DEAL:2d/2d::     at location 0.650000 0.650000
-DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 4 is in cell 2.3
+DEAL:2d/2d::     at location 0.450000 0.450000
+DEAL:2d/2d::     at reference location 0.800000 0.800000
 DEAL:2d/2d::Particle id 5 is in cell 2.12
 DEAL:2d/2d::     at location 0.550000 0.550000
 DEAL:2d/2d::     at reference location 0.200000 0.200000
-DEAL:2d/2d::Particle id 9 is in cell 2.15
-DEAL:2d/2d::     at location 0.950000 0.950000
-DEAL:2d/2d::     at reference location 0.800000 0.800000
+DEAL:2d/2d::Particle id 6 is in cell 2.12
+DEAL:2d/2d::     at location 0.650000 0.650000
+DEAL:2d/2d::     at reference location 0.600000 0.600000
+DEAL:2d/2d::Particle id 7 is in cell 2.12
+DEAL:2d/2d::     at location 0.750000 0.750000
+DEAL:2d/2d::     at reference location 1.00000 1.00000
 DEAL:2d/2d::Particle id 8 is in cell 2.15
 DEAL:2d/2d::     at location 0.850000 0.850000
 DEAL:2d/2d::     at reference location 0.400000 0.400000
+DEAL:2d/2d::Particle id 9 is in cell 2.15
+DEAL:2d/2d::     at location 0.950000 0.950000
+DEAL:2d/2d::     at reference location 0.800000 0.800000
 DEAL:2d/2d::OK
 DEAL:3d/3d::Particle number: 10
-DEAL:3d/3d::Particle id 2 is in cell 2.0
-DEAL:3d/3d::     at location 0.250000 0.250000 0.250000
-DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
-DEAL:3d/3d::Particle id 1 is in cell 2.0
-DEAL:3d/3d::     at location 0.150000 0.150000 0.150000
-DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
 DEAL:3d/3d::Particle id 0 is in cell 2.0
 DEAL:3d/3d::     at location 0.0500000 0.0500000 0.0500000
 DEAL:3d/3d::     at reference location 0.200000 0.200000 0.200000
-DEAL:3d/3d::Particle id 4 is in cell 2.7
-DEAL:3d/3d::     at location 0.450000 0.450000 0.450000
-DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
+DEAL:3d/3d::Particle id 1 is in cell 2.0
+DEAL:3d/3d::     at location 0.150000 0.150000 0.150000
+DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 2 is in cell 2.0
+DEAL:3d/3d::     at location 0.250000 0.250000 0.250000
+DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
 DEAL:3d/3d::Particle id 3 is in cell 2.7
 DEAL:3d/3d::     at location 0.350000 0.350000 0.350000
 DEAL:3d/3d::     at reference location 0.400000 0.400000 0.400000
-DEAL:3d/3d::Particle id 7 is in cell 2.56
-DEAL:3d/3d::     at location 0.750000 0.750000 0.750000
-DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
-DEAL:3d/3d::Particle id 6 is in cell 2.56
-DEAL:3d/3d::     at location 0.650000 0.650000 0.650000
-DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 4 is in cell 2.7
+DEAL:3d/3d::     at location 0.450000 0.450000 0.450000
+DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
 DEAL:3d/3d::Particle id 5 is in cell 2.56
 DEAL:3d/3d::     at location 0.550000 0.550000 0.550000
 DEAL:3d/3d::     at reference location 0.200000 0.200000 0.200000
-DEAL:3d/3d::Particle id 9 is in cell 2.63
-DEAL:3d/3d::     at location 0.950000 0.950000 0.950000
-DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
+DEAL:3d/3d::Particle id 6 is in cell 2.56
+DEAL:3d/3d::     at location 0.650000 0.650000 0.650000
+DEAL:3d/3d::     at reference location 0.600000 0.600000 0.600000
+DEAL:3d/3d::Particle id 7 is in cell 2.56
+DEAL:3d/3d::     at location 0.750000 0.750000 0.750000
+DEAL:3d/3d::     at reference location 1.00000 1.00000 1.00000
 DEAL:3d/3d::Particle id 8 is in cell 2.63
 DEAL:3d/3d::     at location 0.850000 0.850000 0.850000
 DEAL:3d/3d::     at reference location 0.400000 0.400000 0.400000
+DEAL:3d/3d::Particle id 9 is in cell 2.63
+DEAL:3d/3d::     at location 0.950000 0.950000 0.950000
+DEAL:3d/3d::     at reference location 0.800000 0.800000 0.800000
 DEAL:3d/3d::OK
index 2162263a28c5d72c6922d4eef04e97d6aac60bda..d47d99c63942f3c876e90e820f5577c754b64bca 100644 (file)
 
 
 
-// Like particle_03, but tests the creation and use of a
+// Like particle_handler_serial_01, but tests the creation and use of a
 // particle iterator from the created particle.
 
 #include <deal.II/base/array_view.h>
 
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
 #include <deal.II/particles/particle.h>
 #include <deal.II/particles/particle_iterator.h>
 
@@ -31,6 +36,11 @@ void
 test()
 {
   {
+    Triangulation<dim> tr;
+
+    GridGenerator::hyper_cube(tr);
+    MappingQ<dim> mapping(1);
+
     const unsigned int           n_properties_per_particle = 3;
     Particles::PropertyPool<dim> pool(n_properties_per_particle);
 
@@ -55,16 +65,20 @@ test()
     particle.set_properties(
       ArrayView<double>(&properties[0], properties.size()));
 
-    std::multimap<Particles::internal::LevelInd, Particles::Particle<dim>> map;
+    std::vector<std::vector<Particles::Particle<dim>>> particle_container;
+    particle_container.resize(1);
 
-    Particles::internal::LevelInd level_index = std::make_pair(0, 0);
-    map.insert(std::make_pair(level_index, particle));
+    particle_container[0].push_back(particle);
 
     particle.get_properties()[0] = 0.05;
-    map.insert(std::make_pair(level_index, particle));
-
-    Particles::ParticleIterator<dim> particle_it(map, map.begin());
-    Particles::ParticleIterator<dim> particle_end(map, map.end());
+    particle_container[0].push_back(particle);
+
+    Particles::ParticleIterator<dim> particle_it(particle_container,
+                                                 tr.begin(),
+                                                 0);
+    Particles::ParticleIterator<dim> particle_end(particle_container,
+                                                  tr.end(),
+                                                  0);
 
     for (; particle_it != particle_end; ++particle_it)
       {
diff --git a/tests/particles/particle_iterator_02.cc b/tests/particles/particle_iterator_02.cc
new file mode 100644 (file)
index 0000000..fa4a113
--- /dev/null
@@ -0,0 +1,132 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Like particle_iterator_01, but tests state of the iterator.
+
+#include <deal.II/base/array_view.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/particles/particle.h>
+#include <deal.II/particles/particle_iterator.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test()
+{
+  {
+    Triangulation<dim> tr;
+
+    GridGenerator::hyper_cube(tr);
+    MappingQ<dim> mapping(1);
+
+    const unsigned int           n_properties_per_particle = 3;
+    Particles::PropertyPool<dim> pool(n_properties_per_particle);
+
+    Point<dim> position;
+    position(0) = 0.3;
+
+    if (dim > 1)
+      position(1) = 0.5;
+
+    Point<dim> reference_position;
+    reference_position(0) = 0.2;
+
+    if (dim > 1)
+      reference_position(1) = 0.4;
+
+    const types::particle_index index(7);
+
+    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);
+
+    Particles::ParticleIterator<dim> particle_begin(particle_container,
+                                                    tr.begin(),
+                                                    0);
+    Particles::ParticleIterator<dim> particle_end(particle_container,
+                                                  tr.end(),
+                                                  0);
+    Particles::ParticleIterator<dim> particle_nonexistent1(particle_container,
+                                                           tr.begin(),
+                                                           1);
+    Particles::ParticleIterator<dim> particle_nonexistent2(particle_container,
+                                                           tr.end(),
+                                                           1);
+    Particles::ParticleIterator<dim> particle_invalid;
+
+    Assert(particle_begin->state() == IteratorState::valid, ExcInternalError());
+    Assert(particle_end->state() == IteratorState::past_the_end,
+           ExcInternalError());
+    Assert(particle_nonexistent1->state() == IteratorState::invalid,
+           ExcInternalError());
+    Assert(particle_nonexistent2->state() == IteratorState::invalid,
+           ExcInternalError());
+    Assert(particle_invalid->state() == IteratorState::invalid,
+           ExcInternalError());
+
+    Particles::ParticleIterator<dim> particle_iterator = particle_begin;
+    ++particle_iterator;
+
+    Assert(particle_iterator->state() == IteratorState::past_the_end,
+           ExcInternalError());
+
+    particle_iterator = particle_begin;
+    --particle_iterator;
+
+    Assert(particle_iterator->state() == IteratorState::past_the_end,
+           ExcInternalError());
+
+    for (particle_iterator = particle_begin; particle_iterator != particle_end;
+         ++particle_iterator)
+      {
+        deallog << "Particle position: " << particle_iterator->get_location()
+                << ". Iterator valid: "
+                << (particle_iterator->state() == IteratorState::valid)
+                << std::endl;
+      }
+
+    Assert(particle_iterator->state() == IteratorState::past_the_end,
+           ExcInternalError());
+  }
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+  initlog();
+  test<2>();
+  test<3>();
+}
diff --git a/tests/particles/particle_iterator_02.output b/tests/particles/particle_iterator_02.output
new file mode 100644 (file)
index 0000000..f7d1ea9
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL::Particle position: 0.300000 0.500000. Iterator valid: 1
+DEAL::OK
+DEAL::Particle position: 0.300000 0.500000 0.00000. Iterator valid: 1
+DEAL::OK
index 5403e3fac772e23204f76fc330d7a672d41bc935..a514f6094d6f38b07c84d367642188c988520d4e 100644 (file)
@@ -1,31 +1,31 @@
 
-DEAL::Particle 5: 0.500000 
-DEAL::Particle 4: 0.400000 
-DEAL::Particle 3: 0.300000 
-DEAL::Particle 2: 0.200000 
-DEAL::Particle 1: 0.100000 
 DEAL::Particle 0: 0.00000 
-DEAL::Particle 9: 0.900000 
-DEAL::Particle 8: 0.800000 
-DEAL::Particle 7: 0.700000 
-DEAL::Particle 6: 0.600000 
-DEAL::Particle 5: 0.500000 
-DEAL::Particle 4: 0.400000 
-DEAL::Particle 3: 0.300000 
-DEAL::Particle 2: 0.200000 
 DEAL::Particle 1: 0.100000 
-DEAL::Particle 0: 0.00000 
-DEAL::Particle 9: 0.900000 
-DEAL::Particle 8: 0.800000 
-DEAL::Particle 7: 0.700000 
-DEAL::Particle 6: 0.600000 
-DEAL::Particle 5: 0.500000 
-DEAL::Particle 4: 0.400000 
-DEAL::Particle 3: 0.300000 
 DEAL::Particle 2: 0.200000 
-DEAL::Particle 1: 0.100000 
-DEAL::Particle 0: 0.00000 
-DEAL::Particle 9: 0.900000 
+DEAL::Particle 3: 0.300000 
+DEAL::Particle 4: 0.400000 
+DEAL::Particle 5: 0.500000 
+DEAL::Particle 6: 0.600000 
+DEAL::Particle 7: 0.700000 
 DEAL::Particle 8: 0.800000 
+DEAL::Particle 9: 0.900000 
+DEAL::Particle 0: 0.00000 
+DEAL::Particle 1: 0.100000 
+DEAL::Particle 2: 0.200000 
+DEAL::Particle 3: 0.300000 
+DEAL::Particle 4: 0.400000 
+DEAL::Particle 5: 0.500000 
+DEAL::Particle 6: 0.600000 
 DEAL::Particle 7: 0.700000 
+DEAL::Particle 8: 0.800000 
+DEAL::Particle 9: 0.900000 
+DEAL::Particle 0: 0.00000 
+DEAL::Particle 1: 0.100000 
+DEAL::Particle 2: 0.200000 
+DEAL::Particle 3: 0.300000 
+DEAL::Particle 4: 0.400000 
+DEAL::Particle 5: 0.500000 
 DEAL::Particle 6: 0.600000 
+DEAL::Particle 7: 0.700000 
+DEAL::Particle 8: 0.800000 
+DEAL::Particle 9: 0.900000 

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.