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