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