From b57ad5c7a2888342daf6f3a42ee12514e88f7de0 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 3 Nov 2021 16:19:49 +0100 Subject: [PATCH] ParticleHandler: Store vector of particles per cell in std::list --- include/deal.II/particles/particle_accessor.h | 114 ++--- include/deal.II/particles/particle_handler.h | 70 ++- include/deal.II/particles/particle_iterator.h | 23 +- source/particles/particle_handler.cc | 458 ++++++++++-------- source/particles/property_pool.cc | 8 +- 5 files changed, 358 insertions(+), 315 deletions(-) diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h index e9228cc533..d09ed9cd58 100644 --- a/include/deal.II/particles/particle_accessor.h +++ b/include/deal.II/particles/particle_accessor.h @@ -48,8 +48,9 @@ namespace Particles /** * A type for the storage container for particles. */ - using particle_container = - std::vector::Handle>>; + using particle_container = std::list< + std::pair::Handle>, + typename Triangulation::active_cell_iterator>>; /** * @copydoc Particle::write_particle_data_to_memory @@ -358,10 +359,10 @@ namespace Particles * friend classes. */ ParticleAccessor( - const particle_container & particles, - const PropertyPool &property_pool, - const typename Triangulation::active_cell_iterator &cell, - const unsigned int particle_index_within_cell); + const particle_container & particles, + const typename particle_container::iterator particles_in_cell, + const PropertyPool & property_pool, + const unsigned int particle_index_within_cell); /** * Returns a reference to the current Particle. Because the @@ -378,20 +379,20 @@ namespace Particles get_handle(); /** - * A pointer to the container that stores the particles. Obviously, - * this accessor becomes invalid if the container changes. + * A pointer to the container that stores the particles. */ - typename particle_container::iterator particles_on_cell; + const particle_container *particles; /** - * A pointer to the property pool that stores the actual particle data. + * An iterator to the currently active particles within the particles + * object. */ - PropertyPool *property_pool; + typename particle_container::iterator particles_in_cell; /** - * Cell iterator to the cell of the current particle. + * A pointer to the property pool that stores the actual particle data. */ - typename Triangulation::active_cell_iterator cell; + PropertyPool *property_pool; /** * Local index of the particle within its current cell. @@ -468,9 +469,9 @@ namespace Particles template inline ParticleAccessor::ParticleAccessor() - : particles_on_cell(typename particle_container::iterator()) + : particles(nullptr) + , particles_in_cell(typename particle_container::iterator()) , property_pool(nullptr) - , cell() , particle_index_within_cell(numbers::invalid_unsigned_int) {} @@ -478,21 +479,15 @@ namespace Particles template inline ParticleAccessor::ParticleAccessor( - const particle_container & particles, - const PropertyPool &property_pool, - const typename Triangulation::active_cell_iterator &cell, - const unsigned int particle_index_within_cell) - : property_pool(const_cast *>(&property_pool)) - , cell(cell) + 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) + , property_pool(const_cast *>(&property_pool)) , particle_index_within_cell(particle_index_within_cell) - { - if (cell.state() == IteratorState::valid) - particles_on_cell = - const_cast(&particles)->begin() + - cell->active_cell_index(); - else - particles_on_cell = typename particle_container::iterator(); - } + {} @@ -749,8 +744,10 @@ namespace Particles ParticleAccessor::get_surrounding_cell() const { Assert(state() == IteratorState::valid, ExcInternalError()); + Assert(particles_in_cell->second.state() == IteratorState::valid, + ExcInternalError()); - return cell; + return particles_in_cell->second; } @@ -761,8 +758,10 @@ namespace Particles const Triangulation & /*triangulation*/) const { Assert(state() == IteratorState::valid, ExcInternalError()); + Assert(particles_in_cell->second.state() == IteratorState::valid, + ExcInternalError()); - return cell; + return particles_in_cell->second; } @@ -804,20 +803,10 @@ namespace Particles ++particle_index_within_cell; - if (particle_index_within_cell >= particles_on_cell->size()) + if (particle_index_within_cell >= particles_in_cell->first.size()) { - const bool particle_is_locally_owned = cell->is_locally_owned(); - particle_index_within_cell = 0; - - do - { - ++cell; - ++particles_on_cell; - } - while (cell.state() == IteratorState::valid && - (particles_on_cell->empty() || - cell->is_locally_owned() != particle_is_locally_owned)); + ++particles_in_cell; } } @@ -833,28 +822,9 @@ namespace Particles --particle_index_within_cell; else { - const bool particle_is_locally_owned = cell->is_locally_owned(); - - do - { - --cell; - --particles_on_cell; - - if (cell.state() != IteratorState::valid) - { - particles_on_cell = {}; - break; - } - - if (!particles_on_cell->empty() && - cell->is_locally_owned() == particle_is_locally_owned) - { - particle_index_within_cell = particles_on_cell->size() - 1; - break; - } - } - while (particles_on_cell->size() == 0 || - cell->is_locally_owned() != particle_is_locally_owned); + --particles_in_cell; + if (particles_in_cell != particles->end()) + particle_index_within_cell = particles_in_cell->first.size() - 1; } } @@ -875,7 +845,8 @@ namespace Particles ParticleAccessor::operator==( const ParticleAccessor &other) const { - return (property_pool == other.property_pool) && (cell == other.cell) && + return (property_pool == other.property_pool) && + (particles_in_cell == other.particles_in_cell) && (particle_index_within_cell == other.particle_index_within_cell); } @@ -885,10 +856,11 @@ namespace Particles inline IteratorState::IteratorStates ParticleAccessor::state() const { - if (property_pool != nullptr && cell.state() == IteratorState::valid && - particle_index_within_cell < particles_on_cell->size()) + if (particles != nullptr && property_pool != nullptr && + particles_in_cell != particles->end() && + particle_index_within_cell < particles_in_cell->first.size()) return IteratorState::valid; - else if (cell.state() == IteratorState::past_the_end && + else if (particles_in_cell == particles->end() && particle_index_within_cell == 0) return IteratorState::past_the_end; else @@ -903,7 +875,7 @@ namespace Particles inline typename PropertyPool::Handle & ParticleAccessor::get_handle() { - return (*particles_on_cell)[particle_index_within_cell]; + return particles_in_cell->first[particle_index_within_cell]; } @@ -912,7 +884,7 @@ namespace Particles inline const typename PropertyPool::Handle & ParticleAccessor::get_handle() const { - return (*particles_on_cell)[particle_index_within_cell]; + return particles_in_cell->first[particle_index_within_cell]; } } // namespace Particles diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index f2f8dd5b0f..7bbd290e1a 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -79,8 +79,9 @@ namespace Particles /** * A type for the storage container for particles. */ - using particle_container = - std::vector::Handle>>; + using particle_container = std::list< + std::pair::Handle>, + typename Triangulation::active_cell_iterator>>; /** * Default constructor. @@ -878,6 +879,16 @@ namespace Particles const typename Triangulation::active_cell_iterator &cell, const ArrayView &properties = {}); + /** + * Perform the local insertion operation into the particle container. This + * function is used during the higher-level functions inserting particles. + */ + void + insert_particle( + const typename PropertyPool::Handle handle, + const typename Triangulation::cell_iterator &cell, + particle_container & particles); + /** * Address of the triangulation to work on. */ @@ -902,10 +913,21 @@ namespace Particles std::unique_ptr> property_pool; /** - * Set of particles currently living in the local domain including ghost - * cells, organized by the active cell index of the cell they are in. + * Set of particles currently living in the locally owned cells. */ - particle_container particles; + particle_container owned_particles; + + /** + * Set of particles currently living in the locally owned cells. + */ + particle_container ghost_particles; + + /** + * List from the active cells on the present MPI process to positions in + * either `owned_particles` or `ghost_particles` for fast $\mathcal O(1)$ + * access to the particles of a cell. + */ + std::vector cells_to_particle_cache; /** * This variable stores how many particles are stored globally. It is @@ -1150,15 +1172,13 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::begin() { - if (particles.size() == 0) + if (owned_particles.empty()) 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, *property_pool, cell, 0); - - return end(); + else + return particle_iterator(owned_particles, + owned_particles.begin(), + *property_pool, + 0); } @@ -1176,9 +1196,9 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::end() { - return particle_iterator(particles, + return particle_iterator(owned_particles, + owned_particles.end(), *property_pool, - triangulation->end(), 0); } @@ -1197,15 +1217,13 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::begin_ghost() { - if (particles.size() == 0) - return end(); - - for (const auto &cell : triangulation->active_cell_iterators()) - if (cell->is_locally_owned() == false && - particles[cell->active_cell_index()].size() != 0) - return particle_iterator(particles, *property_pool, cell, 0); - - return end_ghost(); + if (ghost_particles.empty()) + return end_ghost(); + else + return particle_iterator(ghost_particles, + ghost_particles.begin(), + *property_pool, + 0); } @@ -1223,9 +1241,9 @@ namespace Particles inline typename ParticleHandler::particle_iterator ParticleHandler::end_ghost() { - return particle_iterator(particles, + return particle_iterator(ghost_particles, + ghost_particles.end(), *property_pool, - triangulation->end(), 0); } diff --git a/include/deal.II/particles/particle_iterator.h b/include/deal.II/particles/particle_iterator.h index 0c3b415df2..ea89a4368a 100644 --- a/include/deal.II/particles/particle_iterator.h +++ b/include/deal.II/particles/particle_iterator.h @@ -43,7 +43,7 @@ namespace Particles * A type for the storage container for particles. */ using particle_container = - std::vector::Handle>>; + typename ParticleAccessor::particle_container; /** * Empty constructor. Such an object is not usable! @@ -56,10 +56,10 @@ namespace Particles * cell. */ ParticleIterator( - const particle_container & particles, - const PropertyPool & property_pool, - const typename Triangulation::cell_iterator &cell, - const unsigned int particle_index_within_cell); + const particle_container & particles, + const typename particle_container::iterator particles_in_cell, + const PropertyPool & property_pool, + const unsigned int particle_index_within_cell); /** * Dereferencing operator, returns a reference to an accessor. Usage is thus @@ -163,11 +163,14 @@ namespace Particles template inline ParticleIterator::ParticleIterator( - const particle_container & particles, - const PropertyPool & property_pool, - const typename Triangulation::cell_iterator &cell, - const unsigned int particle_index_within_cell) - : accessor(particles, property_pool, cell, particle_index_within_cell) + 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) {} diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index f86f0ceabf..148c599d25 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -56,7 +56,6 @@ namespace Particles : triangulation() , mapping() , property_pool(std::make_unique>(0)) - , particles() , global_number_of_particles(0) , number_of_locally_owned_particles(0) , global_max_particles_per_cell(0) @@ -78,7 +77,7 @@ namespace Particles : triangulation(&triangulation, typeid(*this).name()) , mapping(&mapping, typeid(*this).name()) , property_pool(std::make_unique>(n_properties)) - , particles(triangulation.n_active_cells()) + , cells_to_particle_cache(triangulation.n_active_cells()) , global_number_of_particles(0) , number_of_locally_owned_particles(0) , global_max_particles_per_cell(0) @@ -129,7 +128,7 @@ namespace Particles std::make_unique>(new_triangulation, new_mapping); - particles.resize(triangulation->n_active_cells()); + cells_to_particle_cache.resize(triangulation->n_active_cells()); connect_to_triangulation_signals(); } @@ -161,7 +160,9 @@ namespace Particles global_max_particles_per_cell = particle_handler.global_max_particles_per_cell; next_free_particle_index = particle_handler.next_free_particle_index; - particles = particle_handler.particles; + cells_to_particle_cache = particle_handler.cells_to_particle_cache; + owned_particles = particle_handler.owned_particles; + ghost_particles = particle_handler.ghost_particles; ghost_particles_cache.ghost_particles_by_domain = particle_handler.ghost_particles_cache.ghost_particles_by_domain; @@ -187,14 +188,21 @@ namespace Particles void ParticleHandler::clear_particles() { - for (auto &particles_in_cell : particles) - for (auto &particle : particles_in_cell) + for (auto &particles_in_cell : owned_particles) + for (auto &particle : particles_in_cell.first) if (particle != PropertyPool::invalid_handle) property_pool->deregister_particle(particle); - particles.clear(); + for (auto &particles_in_cell : ghost_particles) + for (auto &particle : particles_in_cell.first) + if (particle != PropertyPool::invalid_handle) + property_pool->deregister_particle(particle); + + cells_to_particle_cache.clear(); + owned_particles.clear(); + ghost_particles.clear(); if (triangulation != nullptr) - particles.resize(triangulation->n_active_cells()); + cells_to_particle_cache.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. @@ -203,68 +211,76 @@ namespace Particles - template - void - ParticleHandler::update_cached_numbers() + namespace { - number_of_locally_owned_particles = 0; - types::particle_index local_max_particle_index = 0; - types::particle_index local_max_particles_per_cell = 0; + 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.second < b.second; + }); + + std::array result = {}; + for (const auto &particles_in_cell : particles) + { + const types::particle_index n_particles_in_cell = + particles_in_cell.first.size(); - for (const auto &cell : triangulation->active_cell_iterators()) - { - const auto &particles_in_cell = particles[cell->active_cell_index()]; + // local_max_particles_per_cell + result[0] = std::max(result[0], n_particles_in_cell); - const types::particle_index n_particles_in_cell = - particles_in_cell.size(); + // number of locally owned particles + result[2] += n_particles_in_cell; - local_max_particles_per_cell = - std::max(local_max_particles_per_cell, n_particles_in_cell); + // local_max_particle_index + for (const auto &particle : particles_in_cell.first) + { + result[1] = std::max(result[1], property_pool.get_id(particle)); + } + } - if (cell->is_locally_owned()) - number_of_locally_owned_particles += n_particles_in_cell; + return result; + } + } // namespace - for (const auto &particle : particles_in_cell) - { - local_max_particle_index = - std::max(local_max_particle_index, - property_pool->get_id(particle)); - } - } - if (const auto parallel_triangulation = - dynamic_cast *>( - &*triangulation)) + 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); + + 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]); + + global_number_of_particles = + dealii::Utilities::MPI::sum(number_of_locally_owned_particles, + triangulation->get_communicator()); + + if (global_number_of_particles == 0) { - global_number_of_particles = dealii::Utilities::MPI::sum( - number_of_locally_owned_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]; - } + next_free_particle_index = 0; + global_max_particles_per_cell = 0; } else { - global_number_of_particles = number_of_locally_owned_particles; - next_free_particle_index = - global_number_of_particles == 0 ? 0 : local_max_particle_index + 1; - global_max_particles_per_cell = local_max_particles_per_cell; + Utilities::MPI::max( + ArrayView(result_owned.data(), 2), + triangulation->get_communicator(), + make_array_view(result_owned.data(), result_owned.data() + 2)); + + next_free_particle_index = result_owned[1] + 1; + global_max_particles_per_cell = result_owned[0]; } } @@ -276,12 +292,16 @@ namespace Particles const typename Triangulation::active_cell_iterator &cell) const { - if (particles.size() == 0) + if (cells_to_particle_cache.size() == 0) return 0; if (cell->is_artificial() == false) { - return particles[cell->active_cell_index()].size(); + return cells_to_particle_cache[cell->active_cell_index()] != + typename particle_container::iterator() ? + cells_to_particle_cache[cell->active_cell_index()] + ->first.size() : + 0; } else AssertThrow(false, @@ -315,24 +335,35 @@ namespace Particles if (cell->is_artificial() == false) { - if (particles[active_cell_index].size() == 0) + if (cells_to_particle_cache[active_cell_index] == + typename particle_container::iterator()) { return boost::make_iterator_range( - particle_iterator(particles, *property_pool, cell, 0), - particle_iterator(particles, *property_pool, cell, 0)); + particle_iterator( + owned_particles, owned_particles.end(), *property_pool, 0), + particle_iterator( + owned_particles, owned_particles.end(), *property_pool, 0)); } else { - particle_iterator begin(particles, *property_pool, cell, 0); - particle_iterator end(particles, - *property_pool, - cell, - particles[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); + typename particle_container::iterator particles = + cells_to_particle_cache[active_cell_index]; + if (cell->is_locally_owned()) + { + return boost::make_iterator_range( + particle_iterator( + owned_particles, particles, *property_pool, 0), + particle_iterator( + owned_particles, ++particles, *property_pool, 0)); + } + else + { + return boost::make_iterator_range( + particle_iterator( + ghost_particles, particles, *property_pool, 0), + particle_iterator( + ghost_particles, ++particles, *property_pool, 0)); + } } } else @@ -351,21 +382,21 @@ namespace Particles ParticleHandler::remove_particle( const ParticleHandler::particle_iterator &particle) { - auto &particles_on_cell = *(particle->particles_on_cell); + auto &particles_on_cell = particle->particles_in_cell->first; // 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); - } + property_pool->deregister_particle(handle); // need to reduce the cached number before deleting, because the iterator // may be invalid after removing the particle even if only // accessing the cell - if (particle->get_surrounding_cell()->is_locally_owned()) + const auto cell = particle->get_surrounding_cell(); + const bool owned_cell = cell->is_locally_owned(); + if (owned_cell) --number_of_locally_owned_particles; if (particles_on_cell.size() > 1) @@ -376,7 +407,12 @@ namespace Particles } else { - particles_on_cell.clear(); + 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(); } } @@ -388,52 +424,22 @@ namespace Particles const std::vector::particle_iterator> &particles_to_remove) { - unsigned int n_particles_removed = 0; - - for (auto particles_on_cell = particles.begin(); - particles_on_cell != particles.end(); - ++particles_on_cell) - { - // 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]->particles_on_cell != - particles_on_cell) - continue; - - const unsigned int n_particles_in_cell = particles_on_cell->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]->particles_on_cell == - particles_on_cell && - 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; - } - else - { - (*particles_on_cell)[move_to] = - std::move((*particles_on_cell)[move_from]); - ++move_to; - } - } - particles_on_cell->resize(move_to); - } + // 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->second < b->particles_in_cell->second || + (a->particles_in_cell->second == b->particles_in_cell->second && + a->particle_index_within_cell > b->particle_index_within_cell); + }); + + for (const auto &particle : copy) + remove_particle(particle); update_cached_numbers(); } @@ -455,6 +461,30 @@ namespace Particles + template + void + ParticleHandler::insert_particle( + const typename PropertyPool::Handle handle, + const typename Triangulation::cell_iterator &cell, + particle_container & particles) + { + 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); + else + { + cache->first.push_back(handle); + Assert(cache->second == cell, ExcInternalError()); + } + } + + + template typename ParticleHandler::particle_iterator ParticleHandler::insert_particle( @@ -462,19 +492,19 @@ namespace Particles const typename Triangulation::active_cell_iterator &cell) { Assert(triangulation != nullptr, ExcInternalError()); - Assert(particles.size() == triangulation->n_active_cells(), + Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(), ExcInternalError()); - Assert( - cell->is_locally_owned(), - ExcMessage( - "You can't insert particles in a cell that is not locally owned.")); - - const unsigned int active_cell_index = cell->active_cell_index(); - particles[active_cell_index].push_back(property_pool->register_particle()); - particle_iterator particle_it(particles, + Assert(cell->is_locally_owned(), + 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, - cell, - particles[active_cell_index].size() - 1); + cache->first.size() - 1); data = particle_it->read_particle_data_from_memory(data); @@ -495,20 +525,20 @@ namespace Particles const ArrayView &properties) { Assert(triangulation != nullptr, ExcInternalError()); - Assert(particles.size() == triangulation->n_active_cells(), + Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(), ExcInternalError()); Assert(cell.state() == IteratorState::valid, ExcInternalError()); - Assert( - cell->is_locally_owned(), - ExcMessage( - "You tried to insert a particle into a cell that is not locally owned. This is not supported.")); - - const unsigned int active_cell_index = cell->active_cell_index(); - particles[active_cell_index].push_back(property_pool->register_particle()); - particle_iterator particle_it(particles, + Assert(cell->is_locally_owned(), + 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, - cell, - particles[active_cell_index].size() - 1); + cache->first.size() - 1); particle_it->set_location(position); particle_it->set_reference_location(reference_position); @@ -1053,7 +1083,7 @@ namespace Particles ParticleHandler::sort_particles_into_subdomains_and_cells() { Assert(triangulation != nullptr, ExcInternalError()); - Assert(particles.size() == triangulation->n_active_cells(), + Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(), ExcInternalError()); // TODO: The current algorithm only works for particles that are in @@ -1079,9 +1109,6 @@ namespace Particles for (const auto &cell : triangulation->active_cell_iterators()) { - 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 evaluate // the reference position of particles. Do not sort particles that are @@ -1092,7 +1119,8 @@ namespace Particles continue; } - auto pic = particles_in_cell(cell); + const unsigned int n_pic = n_particles_in_cell(cell); + auto pic = particles_in_cell(cell); real_locations.clear(); for (const auto &particle : pic) @@ -1265,17 +1293,12 @@ namespace Particles // Mark it for MPI transfer otherwise if (current_cell->is_locally_owned()) { - const unsigned int active_cell_index = - current_cell->active_cell_index(); - - particles[active_cell_index].push_back( - (*out_particle->particles_on_cell) - [out_particle->particle_index_within_cell]); + auto &old = out_particle->particles_in_cell + ->first[out_particle->particle_index_within_cell]; + insert_particle(old, current_cell, owned_particles); // Avoid deallocating the memory of this particle - (*out_particle->particles_on_cell) - [out_particle->particle_index_within_cell] = - PropertyPool::invalid_handle; + old = PropertyPool::invalid_handle; } else { @@ -1294,7 +1317,7 @@ namespace Particles { if (dealii::Utilities::MPI::n_mpi_processes( parallel_triangulation->get_communicator()) > 1) - send_recv_particles(moved_particles, particles, moved_cells); + send_recv_particles(moved_particles, owned_particles, moved_cells); } #endif @@ -1306,8 +1329,14 @@ namespace Particles unsorted_handles.reserve(property_pool->n_registered_slots()); typename PropertyPool::Handle sorted_handle = 0; - for (auto &particles_in_cell : particles) - for (auto &particle : particles_in_cell) + for (auto &particles_in_cell : owned_particles) + for (auto &particle : particles_in_cell.first) + { + unsorted_handles.push_back(particle); + particle = sorted_handle++; + } + for (auto &particles_in_cell : ghost_particles) + for (auto &particle : particles_in_cell.first) { unsorted_handles.push_back(particle); particle = sorted_handle++; @@ -1342,21 +1371,29 @@ namespace Particles #else // Clear ghost particles and their properties for (const auto &cell : triangulation->active_cell_iterators()) - if (cell->is_ghost()) + if (cell->is_ghost() && + cells_to_particle_cache[cell->active_cell_index()] != + typename particle_container::iterator()) { + Assert(cells_to_particle_cache[cell->active_cell_index()]->second == + cell, + ExcInternalError()); // Clear particle properties - for (auto &ghost_particle : particles[cell->active_cell_index()]) + for (auto &ghost_particle : + cells_to_particle_cache[cell->active_cell_index()]->first) property_pool->deregister_particle(ghost_particle); // Clear particles themselves - particles[cell->active_cell_index()].clear(); + ghost_particles.erase( + cells_to_particle_cache[cell->active_cell_index()]); + cells_to_particle_cache[cell->active_cell_index()] = + typename particle_container::iterator(); } // Clear ghost particles cache and invalidate it ghost_particles_cache.ghost_particles_by_domain.clear(); ghost_particles_cache.valid = false; - const std::set ghost_owners = parallel_triangulation->ghost_owners(); for (const auto ghost_owner : ghost_owners) @@ -1398,7 +1435,7 @@ namespace Particles send_recv_particles( ghost_particles_cache.ghost_particles_by_domain, - particles, + ghost_particles, std::map< types::subdomain_id, std::vector< @@ -1407,6 +1444,8 @@ namespace Particles #endif } + + template void ParticleHandler::update_ghost_particles() @@ -1433,7 +1472,7 @@ namespace Particles send_recv_particles_properties_and_location( - ghost_particles_cache.ghost_particles_by_domain, particles); + ghost_particles_cache.ghost_particles_by_domain, ghost_particles); #endif } @@ -1453,7 +1492,7 @@ namespace Particles const bool build_cache) { Assert(triangulation != nullptr, ExcInternalError()); - Assert(received_particles.size() == triangulation->n_active_cells(), + Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(), ExcInternalError()); ghost_particles_cache.valid = build_cache; @@ -1461,11 +1500,9 @@ namespace Particles const auto parallel_triangulation = dynamic_cast *>( &*triangulation); - Assert( - parallel_triangulation, - ExcMessage( - "This function is only implemented for parallel::TriangulationBase " - "objects.")); + Assert(parallel_triangulation, + ExcMessage("This function is only implemented for " + "parallel::TriangulationBase objects.")); // Determine the communication pattern const std::set ghost_owners = @@ -1709,14 +1746,17 @@ namespace Particles const typename Triangulation::active_cell_iterator cell = triangulation->create_cell_iterator(id); - const unsigned int active_cell_index = cell->active_cell_index(); - 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); + insert_particle(property_pool->register_particle(), + cell, + received_particles); + const typename particle_container::iterator &cache = + cells_to_particle_cache[cell->active_cell_index()]; + Assert(cache->second == cell, ExcInternalError()); + + particle_iterator particle_it(received_particles, + cache, + *property_pool, + cache->first.size() - 1); recv_data_it = particle_it->read_particle_data_from_memory(recv_data_it); @@ -1754,9 +1794,6 @@ namespace Particles "This function is only implemented for parallel::TriangulationBase " "objects.")); - Assert(updated_particles.size() == parallel_triangulation->n_active_cells(), - ExcInternalError()); - const auto &neighbors = ghost_particles_cache.neighbors; const auto &send_pointers = ghost_particles_cache.send_pointers; const auto &recv_pointers = ghost_particles_cache.recv_pointers; @@ -1838,11 +1875,14 @@ namespace Particles recv_data_it = recv_particle->read_particle_data_from_memory(recv_data_it); + Assert(recv_particle->particles_in_cell->second->is_ghost(), + ExcInternalError()); + if (load_callback) recv_data_it = load_callback( particle_iterator(updated_particles, + recv_particle->particles_in_cell, *property_pool, - recv_particle->cell, recv_particle->particle_index_within_cell), recv_data_it); } @@ -1932,7 +1972,7 @@ namespace Particles // Resize the container if it is possible without // transferring particles if (number_of_locally_owned_particles == 0) - particles.resize(triangulation->n_active_cells()); + cells_to_particle_cache.resize(triangulation->n_active_cells()); } @@ -2099,13 +2139,15 @@ namespace Particles // If the cell persist or is refined store all particles of the // current cell. { - const unsigned int n_particles = - particles[cell->active_cell_index()].size(); + const unsigned int n_particles = n_particles_in_cell(cell); 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)); + stored_particles_on_cell.push_back(particle_iterator( + owned_particles, + cells_to_particle_cache[cell->active_cell_index()], + *property_pool, + i)); } break; @@ -2120,9 +2162,11 @@ namespace Particles stored_particles_on_cell.reserve( stored_particles_on_cell.size() + n_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(particles, *property_pool, child, i)); + stored_particles_on_cell.push_back(particle_iterator( + owned_particles, cache, *property_pool, i)); } } break; @@ -2193,13 +2237,17 @@ namespace Particles { // we need to find the correct child to store the particles and // their reference location has changed - const unsigned int stored_index = - cell_to_store_particles->active_cell_index(); + const typename particle_container::iterator &cache = + cells_to_particle_cache[cell_to_store_particles + ->active_cell_index()]; // 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 < particles[stored_index].size();) + for (unsigned int i = 0; i < cache->first.size();) { for (unsigned int child_index = 0; child_index < GeometryInfo::max_children_per_cell; @@ -2207,6 +2255,7 @@ namespace Particles { const typename Triangulation::cell_iterator child = cell->child(child_index); + Assert(child->is_locally_owned(), ExcInternalError()); try { @@ -2222,13 +2271,12 @@ namespace Particles // redo the loop; otherwise move on to next particle if (child_index != 0) { - particles[child->active_cell_index()].push_back( - particles[stored_index][i]); + insert_particle(cache->first[i], + child, + owned_particles); - particles[stored_index][i] = - particles[stored_index].back(); - particles[stored_index].resize( - particles[stored_index].size() - 1); + cache->first[i] = cache->first.back(); + cache->first.resize(cache->first.size() - 1); } else { diff --git a/source/particles/property_pool.cc b/source/particles/property_pool.cc index 9a03f3195f..2aa595f232 100644 --- a/source/particles/property_pool.cc +++ b/source/particles/property_pool.cc @@ -221,9 +221,11 @@ namespace Particles Assert(sorted_locations.size() == locations.size() - currently_available_handles.size(), - ExcMessage( - "Number of sorted property handles is not equal to number " - "of currently registered handles.")); + ExcMessage("Number of sorted property handles is not equal to " + "number of currently registered handles: " + + std::to_string(sorted_locations.size()) + " vs " + + std::to_string(locations.size()) + " - " + + std::to_string(currently_available_handles.size()))); locations = std::move(sorted_locations); reference_locations = std::move(sorted_reference_locations); -- 2.39.5