From 5e60400564ff37a0f88f604741eb50d0b5df1920 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Fri, 5 Nov 2021 08:42:10 +0100 Subject: [PATCH] Keep both owned/ghost particles in a single container This is needed to keep _GLIBCXX_DEBUG happy when comparing iterators from cache. This commit also addresses the review comments --- include/deal.II/particles/particle_accessor.h | 87 ++-- include/deal.II/particles/particle_handler.h | 93 ++-- include/deal.II/particles/particle_iterator.h | 7 +- source/particles/particle_handler.cc | 477 ++++++++++-------- tests/particles/particle_iterator_02.cc | 17 +- 5 files changed, 381 insertions(+), 300 deletions(-) diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h index 4b3ff0910c..d5a3b3b3ae 100644 --- a/include/deal.II/particles/particle_accessor.h +++ b/include/deal.II/particles/particle_accessor.h @@ -46,8 +46,37 @@ namespace Particles { public: /** - * Data structure to describes the particles in a given cell. This is used - * inside an std::list in `particle_container`. + * Data structure to describe the particles in a given cell. This is used + * inside an std::list in `particle_container`. The storage of this field + * is typically handled by ParticleHandler, using an std::list of this + * structure. + * + * There are four main reasons for the present design: + * */ struct ParticlesInCell { @@ -62,10 +91,9 @@ namespace Particles ParticlesInCell( const std::vector::Handle> &particles, - const typename Triangulation::active_cell_iterator - &cell_iterator) + const typename Triangulation::active_cell_iterator &cell) : particles(particles) - , cell_iterator(cell_iterator) + , cell(cell) {} /** @@ -76,7 +104,7 @@ namespace Particles /** * The underlying cell. */ - typename Triangulation::active_cell_iterator cell_iterator; + typename Triangulation::active_cell_iterator cell; }; /** @@ -391,7 +419,6 @@ namespace Particles * friend classes. */ ParticleAccessor( - const particle_container & particles, const typename particle_container::iterator particles_in_cell, const PropertyPool & property_pool, const unsigned int particle_index_within_cell); @@ -411,13 +438,8 @@ namespace Particles get_handle(); /** - * A pointer to the container that stores the particles. - */ - const particle_container *particles; - - /** - * An iterator to the currently active particles within the particles - * object. + * An iterator to the particles in the current cell within the + * particle_container object. */ typename particle_container::iterator particles_in_cell; @@ -501,8 +523,7 @@ namespace Particles template inline ParticleAccessor::ParticleAccessor() - : particles(nullptr) - , particles_in_cell(typename particle_container::iterator()) + : particles_in_cell(typename particle_container::iterator()) , property_pool(nullptr) , particle_index_within_cell(numbers::invalid_unsigned_int) {} @@ -511,12 +532,10 @@ namespace Particles template inline ParticleAccessor::ParticleAccessor( - const particle_container & particles, const typename particle_container::iterator particles_in_cell, const PropertyPool & property_pool, const unsigned int particle_index_within_cell) - : particles(&particles) - , particles_in_cell(particles_in_cell) + : particles_in_cell(particles_in_cell) , property_pool(const_cast *>(&property_pool)) , particle_index_within_cell(particle_index_within_cell) {} @@ -776,10 +795,10 @@ namespace Particles ParticleAccessor::get_surrounding_cell() const { Assert(state() == IteratorState::valid, ExcInternalError()); - Assert(particles_in_cell->cell_iterator.state() == IteratorState::valid, + Assert(particles_in_cell->cell.state() == IteratorState::valid, ExcInternalError()); - return particles_in_cell->cell_iterator; + return particles_in_cell->cell; } @@ -790,10 +809,10 @@ namespace Particles const Triangulation & /*triangulation*/) const { Assert(state() == IteratorState::valid, ExcInternalError()); - Assert(particles_in_cell->cell_iterator.state() == IteratorState::valid, + Assert(particles_in_cell->cell.state() == IteratorState::valid, ExcInternalError()); - return particles_in_cell->cell_iterator; + return particles_in_cell->cell; } @@ -854,17 +873,10 @@ namespace Particles --particle_index_within_cell; else { - if (particles_in_cell == particles->begin()) - { - particles_in_cell = - const_cast(particles)->end(); - } - else - { - --particles_in_cell; - particle_index_within_cell = - particles_in_cell->particles.size() - 1; - } + --particles_in_cell; + particle_index_within_cell = particles_in_cell->particles.empty() ? + 0 : + particles_in_cell->particles.size() - 1; } } @@ -896,11 +908,12 @@ namespace Particles inline IteratorState::IteratorStates ParticleAccessor::state() const { - if (particles != nullptr && property_pool != nullptr && - particles_in_cell != particles->end() && + if (property_pool != nullptr && + particles_in_cell->cell.state() == IteratorState::valid && particle_index_within_cell < particles_in_cell->particles.size()) return IteratorState::valid; - else if (particles != nullptr && particles_in_cell == particles->end() && + else if (property_pool != nullptr && + particles_in_cell->cell.state() == IteratorState::past_the_end && particle_index_within_cell == 0) return IteratorState::past_the_end; else diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index ddc2859621..1b6e5b5f45 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -880,13 +880,25 @@ namespace Particles /** * Perform the local insertion operation into the particle container. This - * function is used during the higher-level functions inserting particles. + * function is used in the higher-level functions inserting particles. */ - void + particle_iterator insert_particle( const typename PropertyPool::Handle handle, - const typename Triangulation::cell_iterator &cell, - particle_container & particles); + const typename Triangulation::cell_iterator &cell); + + /** + * Delete all entries in the particles container, and then set the three + * anchor entries used to distinguish between owned and ghost cells: We + * add one item to the front, one between and one as the last element in + * the particle container to be able to iterate across particles without + * `if` statements, solely relying on ParticleAccessor::operator== to + * terminate operations, and using the `cell_iterator` inside the particle + * container to check for valid states. + */ + void + reset_particle_container(const Triangulation *trianulation, + particle_container & particles); /** * Address of the triangulation to work on. @@ -912,14 +924,15 @@ namespace Particles std::unique_ptr> property_pool; /** - * Set of particles currently living in the locally owned cells. + * Set of particles currently living in the locally owned or ghost cells. */ - particle_container owned_particles; + particle_container particles; /** - * Set of particles currently living in the ghost cells. + * Iterator to the start of the end of the list elements of + * particle_container which belong to locally owned elements. */ - particle_container ghost_particles; + typename particle_container::iterator owned_particles_end; /** * List from the active cells on the present MPI process to positions in @@ -1018,15 +1031,11 @@ namespace Particles * Transfer particles that have crossed subdomain boundaries to other * processors. * All received particles and their new cells will be appended to the - * @p received_particles vector. + * class variable `particles` at the right slot. * * @param [in] particles_to_send All particles that should be sent and * their new subdomain_ids are in this map. * - * @param [in,out] received_particles Particle container that stores all received - * particles. Note that it is not required nor checked that the container - * is empty, received particles are simply inserted into the container. - * * @param [in] new_cells_for_particles Optional vector of cell * iterators with the same structure as @p particles_to_send. If this * parameter is given it should contain the cell iterator for every @@ -1044,8 +1053,7 @@ namespace Particles void send_recv_particles( const std::map> - & particles_to_send, - particle_container &received_particles, + &particles_to_send, const std::map< types::subdomain_id, std::vector< @@ -1057,27 +1065,20 @@ namespace Particles const bool enable_cache = false); /** - * Transfer ghost particles' position and properties assuming that - * the particles have not changed cells. This routine uses the + * Transfer ghost particles' position and properties assuming that the + * particles have not changed cells. This routine uses the * GhostParticlePartitioner as a caching structure to know which particles - * are ghost to other processes, and where they need to be sent. - * It inherently assumes that particles cannot have changed cell. + * are ghost to other processes, and where they need to be sent. It + * inherently assumes that particles cannot have changed cell, and writes + * the result back to the `particles` member variable. * * @param [in] particles_to_send All particles for which information * should be sent and their new subdomain_ids are in this map. - * - * @param [in,out] received_particles A map with all received - * particles. Note that it is not required nor checked that the container - * is empty, received particles are simply inserted into - * the container. - * */ void send_recv_particles_properties_and_location( const std::map> - & particles_to_send, - particle_container &received_particles); - + &particles_to_send); #endif @@ -1171,13 +1172,10 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::begin() { - if (owned_particles.empty()) - return end(); - else - return particle_iterator(owned_particles, - owned_particles.begin(), - *property_pool, - 0); + Assert(!particles.empty(), ExcInternalError()); + + // jump over the first entry which is used for book-keeping + return particle_iterator(++particles.begin(), *property_pool, 0); } @@ -1195,10 +1193,7 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::end() { - return particle_iterator(owned_particles, - owned_particles.end(), - *property_pool, - 0); + return particle_iterator(owned_particles_end, *property_pool, 0); } @@ -1216,13 +1211,13 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::begin_ghost() { - if (ghost_particles.empty()) - return end_ghost(); - else - return particle_iterator(ghost_particles, - ghost_particles.begin(), - *property_pool, - 0); + Assert(!particles.empty(), ExcInternalError()); + // find the start of ghost particles as one past the end of owned indices; + // the index in between is used for book-keeping + typename particle_container::iterator ghost_particles_begin = + owned_particles_end; + ++ghost_particles_begin; + return particle_iterator(ghost_particles_begin, *property_pool, 0); } @@ -1240,10 +1235,8 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::end_ghost() { - return particle_iterator(ghost_particles, - ghost_particles.end(), - *property_pool, - 0); + // skip the last entry which is used for book-keeping + return particle_iterator(--particles.end(), *property_pool, 0); } diff --git a/include/deal.II/particles/particle_iterator.h b/include/deal.II/particles/particle_iterator.h index ea89a4368a..75fe59acb3 100644 --- a/include/deal.II/particles/particle_iterator.h +++ b/include/deal.II/particles/particle_iterator.h @@ -56,7 +56,6 @@ namespace Particles * cell. */ ParticleIterator( - const particle_container & particles, const typename particle_container::iterator particles_in_cell, const PropertyPool & property_pool, const unsigned int particle_index_within_cell); @@ -163,14 +162,10 @@ namespace Particles template inline ParticleIterator::ParticleIterator( - const particle_container & particles, const typename particle_container::iterator particles_in_cell, const PropertyPool & property_pool, const unsigned int particle_index_within_cell) - : accessor(particles, - particles_in_cell, - property_pool, - particle_index_within_cell) + : accessor(particles_in_cell, property_pool, particle_index_within_cell) {} diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index a75d95a981..c00f560faf 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -51,6 +51,8 @@ namespace Particles } } // namespace + + template ParticleHandler::ParticleHandler() : triangulation() @@ -65,7 +67,9 @@ namespace Particles , load_callback() , handle(numbers::invalid_unsigned_int) , tria_listeners() - {} + { + reset_particle_container(nullptr, particles); + } @@ -77,7 +81,7 @@ namespace Particles : triangulation(&triangulation, typeid(*this).name()) , mapping(&mapping, typeid(*this).name()) , property_pool(std::make_unique>(n_properties)) - , cells_to_particle_cache(triangulation.n_active_cells()) + , cells_to_particle_cache(triangulation.n_active_cells(), particles.end()) , global_number_of_particles(0) , number_of_locally_owned_particles(0) , global_max_particles_per_cell(0) @@ -91,6 +95,7 @@ namespace Particles mapping)) , tria_listeners() { + reset_particle_container(&triangulation, particles); connect_to_triangulation_signals(); } @@ -119,6 +124,8 @@ namespace Particles triangulation = &new_triangulation; mapping = &new_mapping; + reset_particle_container(&new_triangulation, particles); + // Create the memory pool that will store all particle properties property_pool = std::make_unique>(n_properties); @@ -128,7 +135,8 @@ namespace Particles std::make_unique>(new_triangulation, new_mapping); - cells_to_particle_cache.resize(triangulation->n_active_cells()); + cells_to_particle_cache.resize(triangulation->n_active_cells(), + particles.end()); connect_to_triangulation_signals(); } @@ -140,9 +148,6 @@ namespace Particles ParticleHandler::copy_from( const ParticleHandler &particle_handler) { - // clear this object before copying particles - clear(); - const unsigned int n_properties = particle_handler.property_pool->n_properties_per_slot(); initialize(*particle_handler.triangulation, @@ -160,9 +165,16 @@ namespace Particles global_max_particles_per_cell = particle_handler.global_max_particles_per_cell; next_free_particle_index = particle_handler.next_free_particle_index; - cells_to_particle_cache = particle_handler.cells_to_particle_cache; - owned_particles = particle_handler.owned_particles; - ghost_particles = particle_handler.ghost_particles; + particles = particle_handler.particles; + owned_particles_end = particles.begin(); + for (auto it = particle_handler.particles.begin(); + it != particle_handler.owned_particles_end; + ++it) + ++owned_particles_end; + + for (auto it = particles.begin(); it != particles.end(); ++it) + if (!it->particles.empty()) + cells_to_particle_cache[it->cell->active_cell_index()] = it; ghost_particles_cache.ghost_particles_by_domain = particle_handler.ghost_particles_cache.ghost_particles_by_domain; @@ -188,21 +200,20 @@ namespace Particles void ParticleHandler::clear_particles() { - for (auto &particles_in_cell : owned_particles) - for (auto &particle : particles_in_cell.particles) - if (particle != PropertyPool::invalid_handle) - property_pool->deregister_particle(particle); - - for (auto &particles_in_cell : ghost_particles) + for (auto &particles_in_cell : particles) for (auto &particle : particles_in_cell.particles) if (particle != PropertyPool::invalid_handle) property_pool->deregister_particle(particle); cells_to_particle_cache.clear(); - owned_particles.clear(); - ghost_particles.clear(); if (triangulation != nullptr) - cells_to_particle_cache.resize(triangulation->n_active_cells()); + { + cells_to_particle_cache.resize(triangulation->n_active_cells(), + particles.end()); + reset_particle_container(&*triangulation, particles); + } + else + reset_particle_container(nullptr, particles); // the particle properties have already been deleted by their destructor, // but the memory is still allocated. Return the memory as well. @@ -211,69 +222,144 @@ namespace Particles - namespace + template + void + ParticleHandler::reset_particle_container( + const Triangulation *triangulation, + particle_container & given_particles) { - template - std::array - my_update_cached_numbers(const PropertyPool &property_pool, - ParticleContainer & particles) - { - // sort the particles by the active cell index - particles.sort([](const typename ParticleContainer::value_type &a, - const typename ParticleContainer::value_type &b) { - return a.cell_iterator < b.cell_iterator; - }); - - std::array result = {}; - for (const auto &particles_in_cell : particles) - { - const types::particle_index n_particles_in_cell = - particles_in_cell.particles.size(); - - // local_max_particles_per_cell - result[0] = std::max(result[0], n_particles_in_cell); + TriaActiveIterator> past_the_end_iterator( + triangulation, -1, -1); + given_particles.clear(); + for (unsigned int i = 0; i < 3; ++i) + given_particles.emplace_back( + std::vector::Handle>(), + past_the_end_iterator); + + // Set the end of owned particles to the middle of the three elements + owned_particles_end = ++given_particles.begin(); + } - // number of locally owned particles - result[2] += n_particles_in_cell; - // local_max_particle_index - for (const auto &particle : particles_in_cell.particles) + template + void + ParticleHandler::update_cached_numbers() + { + // first sort the owned particles by the active cell index + bool sort_is_necessary = false; + { + auto previous = ++particles.begin(); + for (auto next = previous; next != owned_particles_end; ++next) + { + if (previous->cell.state() == IteratorState::valid && + next->cell.state() == IteratorState::valid && + previous->cell > next->cell) { - result[1] = std::max(result[1], property_pool.get_id(particle)); + sort_is_necessary = true; + break; } + previous = next; } - - return result; } - } // namespace - + if (sort_is_necessary) + { + // we could have tried to call std::list::sort with a custom + // comparator to get things sorted, but things are complicated by the + // three anchor entries that we do not want to move, and we would + // hence pay an O(N log(N)) algorithm with a large constant in front + // of it (on the order of 20+ instructions with many + // difficult-to-predict branches). Therefore, we simply copy the list + // into a new one (keeping alive the possibly large vectors with + // particles on cells) into a new container. + particle_container sorted_particles; + reset_particle_container(&*triangulation, sorted_particles); + + // iterate over cells and insert the entries in the new order + for (const auto &cell : triangulation->active_cell_iterators()) + if (!cell->is_artificial()) + if (cells_to_particle_cache[cell->active_cell_index()] != + particles.end()) + { + // observe that owned_particles_end was already set to point + // into the new sorted_particles array + typename particle_container::iterator insert_position = + cell->is_locally_owned() ? owned_particles_end : + --sorted_particles.end(); + typename particle_container::iterator new_entry = + sorted_particles.insert( + insert_position, typename particle_container::value_type()); + new_entry->cell = cell; + new_entry->particles = + std::move(cells_to_particle_cache[cell->active_cell_index()] + ->particles); + } + particles = std::move(sorted_particles); + + // refresh cells_to_particle_cache + cells_to_particle_cache.clear(); + cells_to_particle_cache.resize(triangulation->n_active_cells(), + particles.end()); + for (auto it = particles.begin(); it != particles.end(); ++it) + if (!it->particles.empty()) + cells_to_particle_cache[it->cell->active_cell_index()] = it; + } - template - void - ParticleHandler::update_cached_numbers() - { - // first compute local result with the function above and then compute the - // collective results - const std::array result_ghosted = - my_update_cached_numbers(*property_pool, ghost_particles); - std::array result_owned = - my_update_cached_numbers(*property_pool, owned_particles); + // Ensure that we did not accidentally modify the anchor entries with + // special purpose. + Assert(particles.front().cell.state() == IteratorState::past_the_end && + particles.front().particles.empty() && + particles.back().cell.state() == IteratorState::past_the_end && + particles.back().particles.empty() && + owned_particles_end->cell.state() == IteratorState::past_the_end && + owned_particles_end->particles.empty(), + ExcInternalError()); #ifdef DEBUG + // check that no cache element hits the three anchor states in the list of + // particles + for (const auto &it : cells_to_particle_cache) + Assert(it != particles.begin() && it != owned_particles_end && + it != --(particles.end()), + ExcInternalError()); + + // check that we only have locally owned particles in the first region of + // cells; note that we skip the very first anchor element + for (auto it = ++particles.begin(); it != owned_particles_end; ++it) + Assert(it->cell->is_locally_owned(), ExcInternalError()); + + // check that the cache is consistent with the iterators std::vector verify_cache( - triangulation->n_active_cells()); - for (auto it = owned_particles.begin(); it != owned_particles.end(); ++it) - verify_cache[it->cell_iterator->active_cell_index()] = it; - for (auto it = ghost_particles.begin(); it != ghost_particles.end(); ++it) - verify_cache[it->cell_iterator->active_cell_index()] = it; + triangulation->n_active_cells(), particles.end()); + for (auto it = particles.begin(); it != particles.end(); ++it) + if (!it->particles.empty()) + verify_cache[it->cell->active_cell_index()] = it; for (unsigned int i = 0; i < verify_cache.size(); ++i) Assert(verify_cache[i] == cells_to_particle_cache[i], ExcInternalError()); #endif - number_of_locally_owned_particles = result_owned[2]; - result_owned[0] = std::max(result_owned[0], result_ghosted[0]); - result_owned[1] = std::max(result_owned[1], result_ghosted[1]); + // now compute local result with the function above and then compute the + // collective results + number_of_locally_owned_particles = 0; + + types::particle_index result[2] = {}; + for (const auto &particles_in_cell : particles) + { + const types::particle_index n_particles_in_cell = + particles_in_cell.particles.size(); + + // local_max_particles_per_cell + result[0] = std::max(result[0], n_particles_in_cell); + + // number of locally owned particles + if (n_particles_in_cell > 0 && + particles_in_cell.cell->is_locally_owned()) + number_of_locally_owned_particles += n_particles_in_cell; + + // local_max_particle_index + for (const auto &particle : particles_in_cell.particles) + result[1] = std::max(result[1], property_pool->get_id(particle)); + } global_number_of_particles = dealii::Utilities::MPI::sum(number_of_locally_owned_particles, @@ -286,13 +372,10 @@ namespace Particles } else { - Utilities::MPI::max( - ArrayView(result_owned.data(), 2), - triangulation->get_communicator(), - make_array_view(result_owned.data(), result_owned.data() + 2)); + Utilities::MPI::max(result, triangulation->get_communicator(), result); - next_free_particle_index = result_owned[1] + 1; - global_max_particles_per_cell = result_owned[0]; + next_free_particle_index = result[1] + 1; + global_max_particles_per_cell = result[0]; } } @@ -310,7 +393,7 @@ namespace Particles if (cell->is_artificial() == false) { return cells_to_particle_cache[cell->active_cell_index()] != - typename particle_container::iterator() ? + particles.end() ? cells_to_particle_cache[cell->active_cell_index()] ->particles.size() : 0; @@ -347,35 +430,23 @@ namespace Particles if (cell->is_artificial() == false) { - if (cells_to_particle_cache[active_cell_index] == - typename particle_container::iterator()) + if (cells_to_particle_cache[active_cell_index] == particles.end()) { return boost::make_iterator_range( - particle_iterator( - owned_particles, owned_particles.end(), *property_pool, 0), - particle_iterator( - owned_particles, owned_particles.end(), *property_pool, 0)); + particle_iterator(particles.begin(), *property_pool, 0), + particle_iterator(particles.begin(), *property_pool, 0)); } else { - const typename particle_container::iterator particles = - cells_to_particle_cache[active_cell_index]; - typename particle_container::iterator next = particles; - ++next; - if (cell->is_locally_owned()) - { - return boost::make_iterator_range( - particle_iterator( - owned_particles, particles, *property_pool, 0), - particle_iterator(owned_particles, next, *property_pool, 0)); - } - else - { - return boost::make_iterator_range( - particle_iterator( - ghost_particles, particles, *property_pool, 0), - particle_iterator(ghost_particles, next, *property_pool, 0)); - } + const typename particle_container::iterator + particles_in_current_cell = + cells_to_particle_cache[active_cell_index]; + typename particle_container::iterator particles_in_next_cell = + particles_in_current_cell; + ++particles_in_next_cell; + return boost::make_iterator_range( + particle_iterator(particles_in_current_cell, *property_pool, 0), + particle_iterator(particles_in_next_cell, *property_pool, 0)); } } else @@ -419,12 +490,8 @@ namespace Particles } else { - if (owned_cell) - owned_particles.erase(particle->particles_in_cell); - else - ghost_particles.erase(particle->particles_in_cell); - cells_to_particle_cache[cell->active_cell_index()] = - typename particle_container::iterator(); + particles.erase(particle->particles_in_cell); + cells_to_particle_cache[cell->active_cell_index()] = particles.end(); } } @@ -436,24 +503,47 @@ namespace Particles const std::vector::particle_iterator> &particles_to_remove) { - // sort particles in reverse order on each cell to ensure the removal with - // the other function is valid - std::vector::particle_iterator> copy( - particles_to_remove); - std::sort(copy.begin(), - copy.end(), - [&](const ParticleHandler::particle_iterator &a, - const ParticleHandler::particle_iterator &b) { - return a->particles_in_cell->cell_iterator < - b->particles_in_cell->cell_iterator || - (a->particles_in_cell->cell_iterator == - b->particles_in_cell->cell_iterator && - a->particle_index_within_cell > - b->particle_index_within_cell); - }); - - for (const auto &particle : copy) - remove_particle(particle); + // We need to remove particles backwards on each cell to keep the particle + // iterators alive as we keep removing particles on the same cell. To + // ensure that this is safe, we either check if we already have sorted + // iterators or if we need to manually sort + const auto check_greater = [](const particle_iterator &a, + const particle_iterator &b) { + return a->particles_in_cell->cell > b->particles_in_cell->cell || + (a->particles_in_cell->cell == b->particles_in_cell->cell && + a->particle_index_within_cell > b->particle_index_within_cell); + }; + + bool particles_are_sorted = true; + auto previous = particles_to_remove.begin(); + for (auto next = previous; next != particles_to_remove.end(); ++next) + { + if (check_greater(*previous, *next)) + { + particles_are_sorted = false; + break; + } + previous = next; + } + if (particles_are_sorted) + { + // pass along backwards in array + for (auto it = particles_to_remove.rbegin(); + it != particles_to_remove.rend(); + ++it) + remove_particle(*it); + } + else + { + std::vector::particle_iterator> + sorted_particles(particles_to_remove); + std::sort(sorted_particles.begin(), + sorted_particles.end(), + check_greater); + + for (const auto &particle : sorted_particles) + remove_particle(particle); + } update_cached_numbers(); } @@ -476,25 +566,31 @@ namespace Particles template - void + typename ParticleHandler::particle_iterator ParticleHandler::insert_particle( const typename PropertyPool::Handle handle, - const typename Triangulation::cell_iterator &cell, - particle_container & particles) + const typename Triangulation::cell_iterator &cell) { const unsigned int active_cell_index = cell->active_cell_index(); typename particle_container::iterator &cache = cells_to_particle_cache[active_cell_index]; - if (cache == typename particle_container::iterator()) - cache = particles.emplace( - particles.end(), - std::vector::Handle>{handle}, - cell); + if (cache == particles.end()) + { + const typename particle_container::iterator insert_position = + cell->is_locally_owned() ? owned_particles_end : --particles.end(); + cache = particles.emplace( + insert_position, + std::vector::Handle>{handle}, + cell); + } else { cache->particles.push_back(handle); - Assert(cache->cell_iterator == cell, ExcInternalError()); + Assert(cache->cell == cell, ExcInternalError()); } + return particle_iterator(cache, + *property_pool, + cache->particles.size() - 1); } @@ -512,13 +608,8 @@ namespace Particles ExcMessage("You tried to insert particles into a cell that is not " "locally owned. This is not supported.")); - insert_particle(property_pool->register_particle(), cell, owned_particles); - const typename particle_container::iterator &cache = - cells_to_particle_cache[cell->active_cell_index()]; - particle_iterator particle_it(owned_particles, - cache, - *property_pool, - cache->particles.size() - 1); + particle_iterator particle_it = + insert_particle(property_pool->register_particle(), cell); data = particle_it->read_particle_data_from_memory(data); @@ -546,13 +637,8 @@ namespace Particles ExcMessage("You tried to insert particles into a cell that is not " "locally owned. This is not supported.")); - insert_particle(property_pool->register_particle(), cell, owned_particles); - const typename particle_container::iterator &cache = - cells_to_particle_cache[cell->active_cell_index()]; - particle_iterator particle_it(owned_particles, - cache, - *property_pool, - cache->particles.size() - 1); + particle_iterator particle_it = + insert_particle(property_pool->register_particle(), cell); particle_it->set_location(position); particle_it->set_reference_location(reference_position); @@ -1310,7 +1396,7 @@ namespace Particles auto &old = out_particle->particles_in_cell ->particles[out_particle->particle_index_within_cell]; - insert_particle(old, current_cell, owned_particles); + insert_particle(old, current_cell); // Avoid deallocating the memory of this particle old = PropertyPool::invalid_handle; @@ -1332,7 +1418,7 @@ namespace Particles { if (dealii::Utilities::MPI::n_mpi_processes( parallel_triangulation->get_communicator()) > 1) - send_recv_particles(moved_particles, owned_particles, moved_cells); + send_recv_particles(moved_particles, moved_cells); } #endif @@ -1344,13 +1430,7 @@ namespace Particles unsorted_handles.reserve(property_pool->n_registered_slots()); typename PropertyPool::Handle sorted_handle = 0; - for (auto &particles_in_cell : owned_particles) - for (auto &particle : particles_in_cell.particles) - { - unsorted_handles.push_back(particle); - particle = sorted_handle++; - } - for (auto &particles_in_cell : ghost_particles) + for (auto &particles_in_cell : particles) for (auto &particle : particles_in_cell.particles) { unsorted_handles.push_back(particle); @@ -1387,11 +1467,10 @@ namespace Particles // Clear ghost particles and their properties for (const auto &cell : triangulation->active_cell_iterators()) if (cell->is_ghost() && - cells_to_particle_cache[cell->active_cell_index()] != - typename particle_container::iterator()) + cells_to_particle_cache[cell->active_cell_index()] != particles.end()) { - Assert(cells_to_particle_cache[cell->active_cell_index()] - ->cell_iterator == cell, + Assert(cells_to_particle_cache[cell->active_cell_index()]->cell == + cell, ExcInternalError()); // Clear particle properties for (auto &ghost_particle : @@ -1399,10 +1478,8 @@ namespace Particles property_pool->deregister_particle(ghost_particle); // Clear particles themselves - ghost_particles.erase( - cells_to_particle_cache[cell->active_cell_index()]); - cells_to_particle_cache[cell->active_cell_index()] = - typename particle_container::iterator(); + particles.erase(cells_to_particle_cache[cell->active_cell_index()]); + cells_to_particle_cache[cell->active_cell_index()] = particles.end(); } // Clear ghost particles cache and invalidate it @@ -1450,7 +1527,6 @@ namespace Particles send_recv_particles( ghost_particles_cache.ghost_particles_by_domain, - ghost_particles, std::map< types::subdomain_id, std::vector< @@ -1487,7 +1563,7 @@ namespace Particles send_recv_particles_properties_and_location( - ghost_particles_cache.ghost_particles_by_domain, ghost_particles); + ghost_particles_cache.ghost_particles_by_domain); #endif } @@ -1498,8 +1574,7 @@ namespace Particles void ParticleHandler::send_recv_particles( const std::map> - & particles_to_send, - particle_container &received_particles, + &particles_to_send, const std::map< types::subdomain_id, std::vector::active_cell_iterator>> @@ -1761,15 +1836,12 @@ namespace Particles const typename Triangulation::active_cell_iterator cell = triangulation->create_cell_iterator(id); - insert_particle(property_pool->register_particle(), - cell, - received_particles); + insert_particle(property_pool->register_particle(), cell); const typename particle_container::iterator &cache = cells_to_particle_cache[cell->active_cell_index()]; - Assert(cache->cell_iterator == cell, ExcInternalError()); + Assert(cache->cell == cell, ExcInternalError()); - particle_iterator particle_it(received_particles, - cache, + particle_iterator particle_it(cache, *property_pool, cache->particles.size() - 1); @@ -1797,8 +1869,7 @@ namespace Particles void ParticleHandler::send_recv_particles_properties_and_location( const std::map> - & particles_to_send, - particle_container &updated_particles) + &particles_to_send) { const auto parallel_triangulation = dynamic_cast *>( @@ -1890,13 +1961,12 @@ namespace Particles recv_data_it = recv_particle->read_particle_data_from_memory(recv_data_it); - Assert(recv_particle->particles_in_cell->cell_iterator->is_ghost(), + Assert(recv_particle->particles_in_cell->cell->is_ghost(), ExcInternalError()); if (load_callback) recv_data_it = load_callback( - particle_iterator(updated_particles, - recv_particle->particles_in_cell, + particle_iterator(recv_particle->particles_in_cell, *property_pool, recv_particle->particle_index_within_cell), recv_data_it); @@ -1987,7 +2057,8 @@ namespace Particles // Resize the container if it is possible without // transferring particles if (number_of_locally_owned_particles == 0) - cells_to_particle_cache.resize(triangulation->n_active_cells()); + cells_to_particle_cache.resize(triangulation->n_active_cells(), + particles.end()); } @@ -2159,7 +2230,6 @@ namespace Particles for (unsigned int i = 0; i < n_particles; ++i) stored_particles_on_cell.push_back(particle_iterator( - owned_particles, cells_to_particle_cache[cell->active_cell_index()], *property_pool, i)); @@ -2180,8 +2250,8 @@ namespace Particles const typename particle_container::iterator &cache = cells_to_particle_cache[child->active_cell_index()]; for (unsigned int i = 0; i < n_particles; ++i) - stored_particles_on_cell.push_back(particle_iterator( - owned_particles, cache, *property_pool, i)); + stored_particles_on_cell.push_back( + particle_iterator(cache, *property_pool, i)); } } break; @@ -2203,25 +2273,30 @@ namespace Particles const typename Triangulation::CellStatus status, const boost::iterator_range::const_iterator> &data_range) { + if (data_range.begin() == data_range.end()) + return; + const auto cell_to_store_particles = (status != parallel::TriangulationBase::CELL_REFINE) ? cell : cell->child(0); // deserialize particles and insert into local storage - { - const void *data = static_cast(&(*data_range.begin())); - - while (data < &(*data_range.end())) - insert_particle(data, cell_to_store_particles); - - Assert( - data == &(*data_range.end()), - ExcMessage( - "The particle data could not be deserialized successfully. " - "Check that when deserializing the particles you expect the same " - "number of properties that were serialized.")); - } + if (data_range.begin() != data_range.end()) + { + const void *data = static_cast(&(*data_range.begin())); + const void *end = static_cast( + &(*data_range.begin()) + (data_range.end() - data_range.begin())); + + while (data < end) + insert_particle(data, cell_to_store_particles); + + Assert(data == end, + ExcMessage( + "The particle data could not be deserialized successfully. " + "Check that when deserializing the particles you expect " + "the same number of properties that were serialized.")); + } auto loaded_particles_on_cell = particles_in_cell(cell_to_store_particles); @@ -2252,15 +2327,15 @@ namespace Particles { // we need to find the correct child to store the particles and // their reference location has changed - const typename particle_container::iterator &cache = + typename particle_container::iterator &cache = cells_to_particle_cache[cell_to_store_particles ->active_cell_index()]; + // make sure that the call above has inserted an entry + Assert(cache != particles.end(), ExcInternalError()); + // Cannot use range-based loop, because number of particles in cell // is going to change - if (cache == typename particle_container::iterator()) - break; - auto particle = loaded_particles_on_cell.begin(); for (unsigned int i = 0; i < cache->particles.size();) { @@ -2286,9 +2361,7 @@ namespace Particles // redo the loop; otherwise move on to next particle if (child_index != 0) { - insert_particle(cache->particles[i], - child, - owned_particles); + insert_particle(cache->particles[i], child); cache->particles[i] = cache->particles.back(); cache->particles.resize( @@ -2306,6 +2379,12 @@ namespace Particles {} } } + // clean up in case child 0 has no particle left + if (cache->particles.empty()) + { + particles.erase(cache); + cache = particles.end(); + } } break; diff --git a/tests/particles/particle_iterator_02.cc b/tests/particles/particle_iterator_02.cc index 17cc1c082f..0fce6bd3e6 100644 --- a/tests/particles/particle_iterator_02.cc +++ b/tests/particles/particle_iterator_02.cc @@ -68,6 +68,10 @@ test() particle_container.emplace_back( std::vector::Handle>(1, particle), tr.begin_active()); + particle_container.emplace_front( + std::vector::Handle>(), tr.end()); + particle_container.emplace_back( + std::vector::Handle>(), tr.end()); pool.set_location(particle, position); pool.set_reference_location(particle, reference_position); @@ -78,18 +82,15 @@ test() properties.end(), particle_properties.begin()); - Particles::ParticleIterator particle_begin(particle_container, - particle_container.begin(), - pool, - 0); - Particles::ParticleIterator particle_end(particle_container, - particle_container.end(), + Particles::ParticleIterator particle_begin( + ++particle_container.begin(), pool, 0); + Particles::ParticleIterator particle_end(--particle_container.end(), pool, 0); Particles::ParticleIterator particle_nonexistent1( - particle_container, particle_container.begin(), pool, 1); + ++particle_container.begin(), pool, 1); Particles::ParticleIterator particle_nonexistent2( - particle_container, particle_container.end(), pool, 1); + --particle_container.end(), pool, 1); Particles::ParticleIterator particle_invalid; Assert(particle_begin->state() == IteratorState::valid, ExcInternalError()); -- 2.39.5