auto particle = particle_handler.begin();
while (particle != particle_handler.end())
{
- const auto cell =
- particle->get_surrounding_cell(background_triangulation);
+ const auto cell = particle->get_surrounding_cell();
const auto dh_cell =
typename DoFHandler<dim>::cell_iterator(*cell, &fluid_dh);
// the particle itself. We can then assemble the additional
// terms in the system matrix and the right hand side as we would
// normally.
- const auto &cell = particle->get_surrounding_cell(fluid_tria);
+ const auto &cell = particle->get_surrounding_cell();
const auto &dh_cell =
typename DoFHandler<spacedim>::cell_iterator(*cell, &fluid_dh);
dh_cell->get_dof_indices(fluid_dof_indices);
#include <deal.II/base/array_view.h>
#include <deal.II/grid/tria.h>
+#include <deal.II/grid/tria_iterator_base.h>
#include <deal.II/particles/particle.h>
class ParticleAccessor
{
public:
+ /**
+ * A type for the storage container for particles.
+ */
+ using particle_container =
+ std::vector<std::vector<Particle<dim, spacedim>>>;
+
/**
* @copydoc Particle::write_particle_data_to_memory
*/
* but the triangulation itself is not stored in the particle this
* operation requires a reference to the triangulation.
*/
- typename Triangulation<dim, spacedim>::cell_iterator
+ const typename Triangulation<dim, spacedim>::cell_iterator &
+ get_surrounding_cell() const;
+
+ /**
+ * @deprecated: Deprecated version of the function with the same
+ * name above.
+ */
+ DEAL_II_DEPRECATED
+ const typename Triangulation<dim, spacedim>::cell_iterator &
get_surrounding_cell(
const Triangulation<dim, spacedim> &triangulation) const;
bool
operator==(const ParticleAccessor<dim, spacedim> &other) const;
+ /**
+ * Return the state of the accessor.
+ */
+ IteratorState::IteratorStates
+ state() const;
+
private:
/**
* Construct an invalid accessor. Such an object is not usable.
ParticleAccessor();
/**
- * Construct an accessor from a reference to a map and an iterator to the
- * map. This constructor is `private` so that it can only be accessed by
+ * Construct an accessor from a reference to a container, an iterator to the
+ * current cell, and the particle index within that cell.
+ * This constructor is `private` so that it can only be accessed by
* friend classes.
*/
ParticleAccessor(
- const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
- const typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator
- &particle);
+ const particle_container &particles,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+ const unsigned int particle_index_within_cell);
+
+ /**
+ * Returns a reference to the current Particle. Because the
+ * internal structure may change this is not intended for public use and
+ * only a convenience function for internal purposes.
+ */
+ const Particle<dim, spacedim> &
+ get_particle() const;
+
+ /**
+ * Non-@p const version of the function with the same name above.
+ */
+ Particle<dim, spacedim> &
+ get_particle();
- private:
/**
* A pointer to the container that stores the particles. Obviously,
* this accessor is invalidated if the container changes.
*/
- std::multimap<internal::LevelInd, Particle<dim, spacedim>> *map;
+ particle_container *particles;
/**
- * An iterator into the container of particles. Obviously,
- * this accessor is invalidated if the container changes.
+ * Cell iterator to the cell of the current particle.
+ */
+ typename Triangulation<dim, spacedim>::active_cell_iterator cell;
+
+ /**
+ * Index to the cell this particle is stored in at the moment. This could be
+ * read from the member variable cell, but is used in many performance
+ * critical functions and is therefore stored and updated separately.
+ */
+ unsigned int active_cell_index;
+
+ /**
+ * Local index of the particle within its current cell.
*/
- typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator particle;
+ unsigned int particle_index_within_cell;
// Make ParticleIterator a friend to allow it constructing
// ParticleAccessors.
ParticleAccessor<dim, spacedim>::serialize(Archive & ar,
const unsigned int version)
{
- return particle->second.serialize(ar, version);
+ return get_particle().serialize(ar, version);
}
template <int dim, int spacedim>
inline ParticleAccessor<dim, spacedim>::ParticleAccessor()
- : map(nullptr)
- , particle()
+ : particles(nullptr)
+ , cell()
+ , active_cell_index(numbers::invalid_unsigned_int)
+ , particle_index_within_cell(numbers::invalid_unsigned_int)
{}
template <int dim, int spacedim>
inline ParticleAccessor<dim, spacedim>::ParticleAccessor(
- const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
- const typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator & particle)
- : map(const_cast<
- std::multimap<internal::LevelInd, Particle<dim, spacedim>> *>(&map))
- , particle(particle)
- {}
+ const particle_container &particles,
+ const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
+ const unsigned int particle_index_within_cell)
+ : particles(const_cast<particle_container *>(&particles))
+ , cell(cell)
+ , particle_index_within_cell(particle_index_within_cell)
+ {
+ if (cell.state() == IteratorState::valid)
+ active_cell_index = cell->active_cell_index();
+ else if (cell.state() == IteratorState::past_the_end)
+ active_cell_index = particles.size();
+ else
+ active_cell_index = numbers::invalid_unsigned_int;
+ }
ParticleAccessor<dim, spacedim>::read_particle_data_from_memory(
const void *data)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.read_particle_data_from_memory(data);
+ return get_particle().read_particle_data_from_memory(data);
}
ParticleAccessor<dim, spacedim>::write_particle_data_to_memory(
void *data) const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.write_particle_data_to_memory(data);
+ return get_particle().write_particle_data_to_memory(data);
}
inline void
ParticleAccessor<dim, spacedim>::set_location(const Point<spacedim> &new_loc)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- particle->second.set_location(new_loc);
+ get_particle().set_location(new_loc);
}
inline const Point<spacedim> &
ParticleAccessor<dim, spacedim>::get_location() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.get_location();
+ return get_particle().get_location();
}
ParticleAccessor<dim, spacedim>::set_reference_location(
const Point<dim> &new_loc)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- particle->second.set_reference_location(new_loc);
+ get_particle().set_reference_location(new_loc);
}
inline const Point<dim> &
ParticleAccessor<dim, spacedim>::get_reference_location() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.get_reference_location();
+ return get_particle().get_reference_location();
}
inline types::particle_index
ParticleAccessor<dim, spacedim>::get_id() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.get_id();
+ return get_particle().get_id();
}
ParticleAccessor<dim, spacedim>::set_property_pool(
PropertyPool<dim, spacedim> &new_property_pool)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- particle->second.set_property_pool(new_property_pool);
+ get_particle().set_property_pool(new_property_pool);
}
inline bool
ParticleAccessor<dim, spacedim>::has_properties() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.has_properties();
+ return get_particle().has_properties();
}
ParticleAccessor<dim, spacedim>::set_properties(
const std::vector<double> &new_properties)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- particle->second.set_properties(new_properties);
+ get_particle().set_properties(new_properties);
}
ParticleAccessor<dim, spacedim>::set_properties(
const ArrayView<const double> &new_properties)
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- particle->second.set_properties(new_properties);
+ get_particle().set_properties(new_properties);
}
inline const ArrayView<const double>
ParticleAccessor<dim, spacedim>::get_properties() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
+
+ return get_particle().get_properties();
+ }
+
+
+
+ template <int dim, int spacedim>
+ inline const typename Triangulation<dim, spacedim>::cell_iterator &
+ ParticleAccessor<dim, spacedim>::get_surrounding_cell() const
+ {
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.get_properties();
+ return cell;
}
template <int dim, int spacedim>
- inline typename Triangulation<dim, spacedim>::cell_iterator
+ inline const typename Triangulation<dim, spacedim>::cell_iterator &
ParticleAccessor<dim, spacedim>::get_surrounding_cell(
- const Triangulation<dim, spacedim> &triangulation) const
+ const Triangulation<dim, spacedim> & /*triangulation*/) const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- const typename Triangulation<dim, spacedim>::cell_iterator cell(
- &triangulation, particle->first.first, particle->first.second);
return cell;
}
inline const ArrayView<double>
ParticleAccessor<dim, spacedim>::get_properties()
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.get_properties();
+ return get_particle().get_properties();
}
inline std::size_t
ParticleAccessor<dim, spacedim>::serialized_size_in_bytes() const
{
- Assert(particle != map->end(), ExcInternalError());
+ Assert(state() == IteratorState::valid, ExcInternalError());
- return particle->second.serialized_size_in_bytes();
+ return get_particle().serialized_size_in_bytes();
}
inline void
ParticleAccessor<dim, spacedim>::next()
{
- Assert(particle != map->end(), ExcInternalError());
- ++particle;
+ Assert(state() == IteratorState::valid, ExcInternalError());
+
+ ++particle_index_within_cell;
+
+ if (particle_index_within_cell > (*particles)[active_cell_index].size() - 1)
+ {
+ particle_index_within_cell = 0;
+
+ do
+ {
+ ++cell;
+ ++active_cell_index;
+ }
+ while (cell.state() == IteratorState::valid &&
+ (*particles)[active_cell_index].size() == 0);
+ }
}
inline void
ParticleAccessor<dim, spacedim>::prev()
{
- Assert(particle != map->begin(), ExcInternalError());
- --particle;
+ Assert(state() == IteratorState::valid, ExcInternalError());
+
+ if (particle_index_within_cell > 0)
+ --particle_index_within_cell;
+ else
+ {
+ do
+ {
+ --cell;
+ --active_cell_index;
+
+ if (cell.state() != IteratorState::valid)
+ {
+ active_cell_index = particles->size();
+ break;
+ }
+
+ if ((*particles)[active_cell_index].size() > 0)
+ {
+ particle_index_within_cell =
+ (*particles)[active_cell_index].size() - 1;
+ break;
+ }
+ }
+ while ((*particles)[active_cell_index].size() == 0);
+ }
}
ParticleAccessor<dim, spacedim>::
operator!=(const ParticleAccessor<dim, spacedim> &other) const
{
- return (map != other.map) || (particle != other.particle);
+ return !(*this == other);
}
ParticleAccessor<dim, spacedim>::
operator==(const ParticleAccessor<dim, spacedim> &other) const
{
- return (map == other.map) && (particle == other.particle);
+ return (particles == other.particles) && (cell == other.cell) &&
+ (particle_index_within_cell == other.particle_index_within_cell);
+ }
+
+
+
+ template <int dim, int spacedim>
+ inline IteratorState::IteratorStates
+ ParticleAccessor<dim, spacedim>::state() const
+ {
+ if (particles != nullptr && cell.state() == IteratorState::valid &&
+ particle_index_within_cell < (*particles)[active_cell_index].size())
+ return IteratorState::valid;
+ else if (particles != nullptr &&
+ cell.state() == IteratorState::past_the_end &&
+ particle_index_within_cell == 0)
+ return IteratorState::past_the_end;
+ else
+ return IteratorState::invalid;
+
+ return IteratorState::invalid;
+ }
+
+
+
+ template <int dim, int spacedim>
+ inline Particle<dim, spacedim> &
+ ParticleAccessor<dim, spacedim>::get_particle()
+ {
+ return (*particles)[active_cell_index][particle_index_within_cell];
}
+
+ template <int dim, int spacedim>
+ inline const Particle<dim, spacedim> &
+ ParticleAccessor<dim, spacedim>::get_particle() const
+ {
+ return (*particles)[active_cell_index][particle_index_within_cell];
+ }
+
} // namespace Particles
DEAL_II_NAMESPACE_CLOSE
*/
using particle_iterator_range = boost::iterator_range<particle_iterator>;
+ /**
+ * A type for the storage container for particles.
+ */
+ using particle_container =
+ std::vector<std::vector<Particle<dim, spacedim>>>;
+
/**
* Default constructor.
*/
const;
/**
- * Remove a particle pointed to by the iterator.
+ * Remove a particle pointed to by the iterator. Afterwards the iterator
+ * will point to one of the remaining particles in the cell. Note that
+ * particle iterators that point to other particles in the same cell as @p particle
+ * may not longer be valid after this function call.
*/
void
remove_particle(const particle_iterator &particle);
+ /**
+ * Remove a vector of particles indicated by the particle iterators.
+ * The iterators and all other iterators are invalidated during the
+ * function call.
+ */
+ void
+ remove_particles(const std::vector<particle_iterator> &particles);
+
/**
* Insert a particle into the collection of particles. Return an iterator
* to the new position of the particle. This function involves a copy of
* Set of particles currently living in the local domain, organized by
* the level/index of the cell they are in.
*/
- std::multimap<internal::LevelInd, Particle<dim, spacedim>> particles;
+ particle_container particles;
/**
- * Set of particles that currently live in the ghost cells of the local
- * domain, organized by the level/index of the cell they are in. These
- * particles are equivalent to the ghost entries in distributed vectors.
+ * Set of particles currently living in the local domain, organized by
+ * the level/index of the cell they are in.
*/
- std::multimap<internal::LevelInd, Particle<dim, spacedim>> ghost_particles;
+ particle_container ghost_particles;
/**
* This variable stores how many particles are stored globally. It is
*/
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
void
send_recv_particles(
const std::map<types::subdomain_id, std::vector<particle_iterator>>
- &particles_to_send,
- std::multimap<internal::LevelInd, Particle<dim, spacedim>>
- &received_particles,
+ & particles_to_send,
+ std::vector<std::vector<Particle<dim, spacedim>>> &received_particles,
const std::map<
types::subdomain_id,
std::vector<
void
send_recv_particles_properties_and_location(
const std::map<types::subdomain_id, std::vector<particle_iterator>>
- &particles_to_send,
- std::multimap<internal::LevelInd, Particle<dim, spacedim>>
- &received_particles);
+ & particles_to_send,
+ std::vector<std::vector<Particle<dim, spacedim>>> &received_particles);
#endif
inline typename ParticleHandler<dim, spacedim>::particle_iterator
ParticleHandler<dim, spacedim>::begin()
{
- return particle_iterator(particles, particles.begin());
+ if (particles.size() == 0)
+ return end();
+
+ for (const auto &cell : triangulation->active_cell_iterators())
+ if (cell->is_locally_owned() &&
+ particles[cell->active_cell_index()].size() != 0)
+ return particle_iterator(particles, cell, 0);
+
+ return end();
}
inline typename ParticleHandler<dim, spacedim>::particle_iterator
ParticleHandler<dim, spacedim>::end()
{
- return particle_iterator(particles, particles.end());
+ return particle_iterator(particles, triangulation->end(), 0);
}
inline typename ParticleHandler<dim, spacedim>::particle_iterator
ParticleHandler<dim, spacedim>::begin_ghost()
{
- return particle_iterator(ghost_particles, ghost_particles.begin());
+ if (particles.size() == 0)
+ return end();
+
+ for (const auto &cell : triangulation->active_cell_iterators())
+ if (cell->is_locally_owned() == false &&
+ ghost_particles[cell->active_cell_index()].size() != 0)
+ return particle_iterator(ghost_particles, cell, 0);
+
+ return end_ghost();
}
inline typename ParticleHandler<dim, spacedim>::particle_iterator
ParticleHandler<dim, spacedim>::end_ghost()
{
- return particle_iterator(ghost_particles, ghost_particles.end());
+ return particle_iterator(ghost_particles, triangulation->end(), 0);
}
/**
* Constructor of the iterator. Takes a reference to the particle
- * container, and an iterator to the cell-particle pair.
+ * container, an iterator to the cell, and the particle index within that
+ * cell.
*/
ParticleIterator(
- const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
- const typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator
- &particle);
+ const std::vector<std::vector<Particle<dim, spacedim>>> & particles,
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int particle_index_within_cell);
/**
* Dereferencing operator, returns a reference to an accessor. Usage is thus
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
template <int dim, int spacedim>
inline ParticleIterator<dim, spacedim>::ParticleIterator(
- const std::multimap<internal::LevelInd, Particle<dim, spacedim>> &map,
- const typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator & particle)
- : accessor(map, particle)
+ const std::vector<std::vector<Particle<dim, spacedim>>> & particles,
+ const typename Triangulation<dim, spacedim>::cell_iterator &cell,
+ const unsigned int particle_index_within_cell)
+ : accessor(particles, cell, particle_index_within_cell)
{}
}
+
+ template <int dim, int spacedim>
+ inline IteratorState::IteratorStates
+ ParticleIterator<dim, spacedim>::state() const
+ {
+ return accessor.state();
+ }
+
} // namespace Particles
DEAL_II_NAMESPACE_CLOSE
* without clearing the multimap of ghost particles, thus greatly
* reducing the cost of exchanging the ghost particles information.
*/
- std::vector<typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator>
- ghost_particles_iterators;
+ std::vector<particle_iterator> ghost_particles_iterators;
/**
* Temporary storage that holds the data of the particles to be sent
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();
while (particle != particle_handler.end())
{
- const auto &cell = particle->get_surrounding_cell(tria);
+ const auto &cell = particle->get_surrounding_cell();
const auto &dh_cell =
typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &field_dh);
dh_cell->get_dof_indices(dof_indices);
, 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()
: triangulation(&triangulation, typeid(*this).name())
, mapping(&mapping, typeid(*this).name())
, property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
- , particles()
- , ghost_particles()
+ , particles(triangulation.n_active_cells())
+ , ghost_particles(triangulation.n_active_cells())
, global_number_of_particles(0)
+ , local_number_of_particles(0)
, global_max_particles_per_cell(0)
, next_free_particle_index(0)
, size_callback()
triangulation_cache =
std::make_unique<GridTools::Cache<dim, spacedim>>(new_triangulation,
new_mapping);
+
+ particles.resize(triangulation->n_active_cells());
+ ghost_particles.resize(triangulation->n_active_cells());
}
ParticleHandler<dim, spacedim>::copy_from(
const ParticleHandler<dim, spacedim> &particle_handler)
{
- // clear and initialize this object before copying particles
+ // clear this object before copying particles
clear();
+
const unsigned int n_properties =
particle_handler.property_pool->n_properties_per_slot();
initialize(*particle_handler.triangulation,
*particle_handler.mapping,
n_properties);
- property_pool->reserve(particle_handler.particles.size() +
- particle_handler.ghost_particles.size());
+
+ property_pool->reserve(particle_handler.n_locally_owned_particles());
// copy static members
global_number_of_particles = particle_handler.global_number_of_particles;
+ local_number_of_particles = particle_handler.local_number_of_particles;
global_max_particles_per_cell =
particle_handler.global_max_particles_per_cell;
next_free_particle_index = particle_handler.next_free_particle_index;
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);
}
{
clear_particles();
global_number_of_particles = 0;
+ local_number_of_particles = 0;
next_free_particle_index = 0;
global_max_particles_per_cell = 0;
}
ParticleHandler<dim, spacedim>::clear_particles()
{
particles.clear();
+ if (triangulation != nullptr)
+ particles.resize(triangulation->n_active_cells());
+
ghost_particles.clear();
+ if (triangulation != nullptr)
+ ghost_particles.resize(triangulation->n_active_cells());
// the particle properties have already been deleted by their destructor,
// but the memory is still allocated. Return the memory as well.
void
ParticleHandler<dim, spacedim>::update_cached_numbers()
{
- types::particle_index locally_highest_index = 0;
- unsigned int local_max_particles_per_cell = 0;
- unsigned int current_particles_per_cell = 0;
- typename Triangulation<dim, spacedim>::active_cell_iterator current_cell =
- triangulation->begin_active();
+ local_number_of_particles = 0;
+ types::particle_index local_max_particle_index = 0;
+ types::particle_index local_max_particles_per_cell = 0;
- for (const auto &particle : *this)
+ for (const auto &particles_in_cell : particles)
{
- locally_highest_index =
- std::max(locally_highest_index, particle.get_id());
+ const types::particle_index n_particles_in_cell =
+ particles_in_cell.size();
- if (particle.get_surrounding_cell(*triangulation) != current_cell)
+ local_max_particles_per_cell =
+ std::max(local_max_particles_per_cell, n_particles_in_cell);
+ local_number_of_particles += n_particles_in_cell;
+
+ for (const auto &particle : particles_in_cell)
{
- current_particles_per_cell = 0;
- current_cell = particle.get_surrounding_cell(*triangulation);
+ local_max_particle_index =
+ std::max(local_max_particle_index, particle.get_id());
}
-
- ++current_particles_per_cell;
- local_max_particles_per_cell =
- std::max(local_max_particles_per_cell, current_particles_per_cell);
}
if (const auto parallel_triangulation =
&*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;
}
}
const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
const
{
- const internal::LevelInd found_cell =
- std::make_pair(cell->level(), cell->index());
+ if (particles.size() == 0)
+ return 0;
- if (cell->is_locally_owned())
- return particles.count(found_cell);
- else if (cell->is_ghost())
- return ghost_particles.count(found_cell);
+ if (cell->is_artificial() == false)
+ {
+ if (cell->is_locally_owned() == true)
+ return particles[cell->active_cell_index()].size();
+ else
+ return ghost_particles[cell->active_cell_index()].size();
+ }
else
AssertThrow(false,
ExcMessage("You can't ask for the particles on an artificial "
ParticleHandler<dim, spacedim>::particles_in_cell(
const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
{
- const internal::LevelInd level_index =
- std::make_pair(cell->level(), cell->index());
+ const unsigned int active_cell_index = cell->active_cell_index();
- if (cell->is_ghost())
+ if (cell->is_artificial() == false)
{
- const auto particles_in_cell = ghost_particles.equal_range(level_index);
- return boost::make_iterator_range(
- particle_iterator(ghost_particles, particles_in_cell.first),
- particle_iterator(ghost_particles, particles_in_cell.second));
- }
- else if (cell->is_locally_owned())
- {
- const auto particles_in_cell = particles.equal_range(level_index);
- return boost::make_iterator_range(
- particle_iterator(particles, particles_in_cell.first),
- particle_iterator(particles, particles_in_cell.second));
+ particle_container &container =
+ (cell->is_locally_owned() == true) ? particles : ghost_particles;
+
+ if (container[active_cell_index].size() == 0)
+ {
+ return boost::make_iterator_range(
+ particle_iterator(container, cell, 0),
+ particle_iterator(container, cell, 0));
+ }
+ else
+ {
+ particle_iterator begin(container, cell, 0);
+ particle_iterator end(container,
+ cell,
+ container[active_cell_index].size() - 1);
+ // end needs to point to the particle after the last one in the
+ // cell.
+ ++end;
+
+ return boost::make_iterator_range(begin, end);
+ }
}
else
AssertThrow(false,
ParticleHandler<dim, spacedim>::remove_particle(
const ParticleHandler<dim, spacedim>::particle_iterator &particle)
{
- particles.erase(particle->particle);
+ const unsigned int active_cell_index = particle->active_cell_index;
+
+ if (particles[active_cell_index].size() > 1)
+ {
+ particles[active_cell_index][particle->particle_index_within_cell] =
+ std::move(particles[active_cell_index].back());
+ particles[active_cell_index].resize(
+ particles[active_cell_index].size() - 1);
+ }
+ else
+ {
+ particles[active_cell_index].clear();
+ }
+
+ --local_number_of_particles;
+ }
+
+
+
+ template <int dim, int spacedim>
+ void
+ ParticleHandler<dim, spacedim>::remove_particles(
+ const std::vector<ParticleHandler<dim, spacedim>::particle_iterator>
+ &particles_to_remove)
+ {
+ unsigned int n_particles_removed = 0;
+
+ for (unsigned int cell_index = 0; cell_index < particles.size();
+ ++cell_index)
+ {
+ // If there is nothing left to remove, break and return
+ if (n_particles_removed == particles_to_remove.size())
+ break;
+
+ // Skip cells where there is nothing to remove
+ if (particles_to_remove[n_particles_removed]->active_cell_index !=
+ cell_index)
+ continue;
+
+ const unsigned int n_particles_in_cell = particles[cell_index].size();
+ unsigned int move_to = 0;
+ for (unsigned int move_from = 0; move_from < n_particles_in_cell;
+ ++move_from)
+ {
+ if (n_particles_removed != particles_to_remove.size() &&
+ particles_to_remove[n_particles_removed]->active_cell_index ==
+ cell_index &&
+ particles_to_remove[n_particles_removed]
+ ->particle_index_within_cell == move_from)
+ {
+ ++n_particles_removed;
+ continue;
+ }
+ else
+ {
+ particles[cell_index][move_to] =
+ std::move(particles[cell_index][move_from]);
+ ++move_to;
+ }
+ }
+ particles[cell_index].resize(move_to);
+ }
+
+ update_cached_numbers();
}
const Particle<dim, spacedim> & particle,
const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
{
- typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator it =
- particles.insert(
- std::make_pair(internal::LevelInd(cell->level(), cell->index()),
- particle));
+ Assert(triangulation != nullptr, ExcInternalError());
+ Assert(cell.state() == IteratorState::valid, ExcInternalError());
- particle_iterator particle_it(particles, it);
+ if (particles.size() == 0)
+ particles.resize(triangulation->n_active_cells());
+
+ const unsigned int active_cell_index = cell->active_cell_index();
+ particles[active_cell_index].push_back(particle);
+ particle_iterator particle_it(particles,
+ cell,
+ particles[active_cell_index].size() - 1);
particle_it->set_property_pool(*property_pool);
- if (particle.has_properties())
- for (unsigned int n = 0; n < particle.get_properties().size(); ++n)
- particle_it->get_properties()[n] = particle.get_properties()[n];
+ ++local_number_of_particles;
return particle_it;
}
typename Triangulation<dim, spacedim>::active_cell_iterator,
Particle<dim, spacedim>> &new_particles)
{
+ Assert(triangulation != nullptr, ExcInternalError());
+ if (particles.size() == 0)
+ particles.resize(triangulation->n_active_cells());
+
for (const auto &particle : new_particles)
{
- // Insert the particle. Store an iterator to the newly
- // inserted particle, and then set its property_pool.
- auto it = particles.insert(
- particles.end(),
- std::make_pair(internal::LevelInd(particle.first->level(),
- particle.first->index()),
- particle.second));
- it->second.set_property_pool(*property_pool);
+ const unsigned int active_cell_index =
+ particle.first->active_cell_index();
+ particles[active_cell_index].push_back(particle.second);
+ particles[active_cell_index].back().set_property_pool(*property_pool);
}
update_cached_numbers();
ParticleHandler<dim, spacedim>::insert_particles(
const std::vector<Point<spacedim>> &positions)
{
+ Assert(triangulation != nullptr, ExcInternalError());
+ if (particles.size() == 0)
+ particles.resize(triangulation->n_active_cells());
+
update_cached_numbers();
// Determine the starting particle index of this process, which
if (cells.size() == 0)
return;
- auto hint =
- particles.find(std::make_pair(cells[0]->level(), cells[0]->index()));
for (unsigned int i = 0; i < cells.size(); ++i)
{
- internal::LevelInd current_cell(cells[i]->level(), cells[i]->index());
+ const unsigned int active_cell_index = cells[i]->active_cell_index();
+
for (unsigned int p = 0; p < local_positions[i].size(); ++p)
{
- hint = particles.insert(
- hint,
- std::make_pair(current_cell,
- Particle<dim, spacedim>(positions[index_map[i][p]],
- local_positions[i][p],
- local_start_index +
- index_map[i][p])));
-
- hint->second.set_property_pool(*property_pool);
+ particles[active_cell_index].push_back(
+ Particle<dim, spacedim>(positions[index_map[i][p]],
+ local_positions[i][p],
+ local_start_index + index_map[i][p]));
+
+ particles[active_cell_index].back().set_property_pool(
+ *property_pool);
}
}
types::particle_index
ParticleHandler<dim, spacedim>::n_locally_owned_particles() const
{
- return particles.size();
+ return local_number_of_particles;
}
// processes around (with an invalid cell).
std::vector<particle_iterator> particles_out_of_cell;
- particles_out_of_cell.reserve(n_locally_owned_particles());
+
+ particles_out_of_cell.reserve(n_locally_owned_particles() / 4);
// Now update the reference locations of the moved particles
std::vector<Point<spacedim>> real_locations;
std::vector<Point<dim>> reference_locations;
- for (auto particle = begin(); particle != end();)
+ real_locations.reserve(global_max_particles_per_cell);
+ reference_locations.reserve(global_max_particles_per_cell);
+
+ for (const auto &cell : triangulation->active_cell_iterators())
{
- const auto cell = particle->get_surrounding_cell(*triangulation);
+ const unsigned int active_cell_index = cell->active_cell_index();
+ const unsigned int n_pic = particles[active_cell_index].size();
+
+ // Particles can be inserted into arbitrary cells, e.g. if their cell is
+ // not known. However, for artificial cells we can not check evaluate
+ // the reference position of particles and particles in ghost cells
+ // should not be owned by this process. Assume all particles that are in
+ // not locally owned cells have to be resorted or transferred.
+ if (cell->is_locally_owned() == false)
+ {
+ for (unsigned int i = 0; i < n_pic; ++i)
+ {
+ particles_out_of_cell.push_back(
+ particle_iterator(particles, cell, i));
+ }
+ continue;
+ }
+
+ auto pic = particles_in_cell(cell);
+
real_locations.clear();
+ for (const auto &particle : pic)
+ real_locations.push_back(particle.get_location());
- // Since we might also work on artificial cells when we initialize the
- // particles on a remote processor, we cannot use the
- // particles_in_cell method. Thus, We instead simply go through the
- // particles and check if the next one belongs to the same cell as the
- // current one.
- for (auto it = particle;
- it != end() && it->get_surrounding_cell(*triangulation) == cell;
- ++it)
- real_locations.push_back(it->get_location());
-
- reference_locations.resize(real_locations.size());
- ArrayView<Point<dim>> reference(reference_locations.data(),
- reference_locations.size());
+ reference_locations.resize(n_pic);
mapping->transform_points_real_to_unit_cell(cell,
real_locations,
- reference);
+ reference_locations);
+ auto particle = pic.begin();
for (const auto &p_unit : reference_locations)
{
if (p_unit[0] == std::numeric_limits<double>::infinity() ||
particles_out_of_cell.push_back(particle);
else
particle->set_reference_location(p_unit);
+
++particle;
}
}
// There are three reasons why a particle is not in its old cell:
// It moved to another cell, to another subdomain or it left the mesh.
- // Particles that moved to another cell are updated and stored inside the
- // sorted_particles vector, particles that moved to another domain are
+ // Particles that moved to another cell are updated and moved inside the
+ // particles vector, particles that moved to another domain are
// collected in the moved_particles_domain vector. Particles that left
// the mesh completely are ignored and removed.
- std::vector<std::pair<internal::LevelInd, Particle<dim, spacedim>>>
- sorted_particles;
std::map<types::subdomain_id, std::vector<particle_iterator>>
moved_particles;
std::map<
// automatic and relatively fast (compared to other parts of this
// algorithm) re-allocation will happen.
using vector_size = typename std::vector<particle_iterator>::size_type;
- sorted_particles.reserve(
- static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
std::set<types::subdomain_id> ghost_owners;
if (const auto parallel_triangulation =
&*triangulation))
ghost_owners = parallel_triangulation->ghost_owners();
- for (const auto ghost_owner : ghost_owners)
+ for (const auto &ghost_owner : ghost_owners)
moved_particles[ghost_owner].reserve(
static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
- for (const auto ghost_owner : ghost_owners)
+ for (const auto &ghost_owner : ghost_owners)
moved_cells[ghost_owner].reserve(
static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
{
// Create a map from vertices to adjacent cells using grid cache
- std::vector<
+ const std::vector<
std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
- vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
+ &vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
- // Create a corresponding map of vectors from vertex to cell center using
- // grid cache
- std::vector<std::vector<Tensor<1, spacedim>>> vertex_to_cell_centers =
- triangulation_cache->get_vertex_to_cell_centers_directions();
+ // Create a corresponding map of vectors from vertex to cell center
+ // using grid cache
+ const std::vector<std::vector<Tensor<1, spacedim>>>
+ &vertex_to_cell_centers =
+ triangulation_cache->get_vertex_to_cell_centers_directions();
std::vector<unsigned int> neighbor_permutation;
// Find the cells that the particles moved to.
for (auto &out_particle : particles_out_of_cell)
{
+ // make a copy of the current cell, since we will modify the
+ // variable current_cell in the following but we need the original in
+ // the case the particle is not found
+ auto current_cell = out_particle->get_surrounding_cell();
+
// The cell the particle is in
Point<dim> current_reference_position;
bool found_cell = false;
// Check if the particle is in one of the old cell's neighbors
// that are adjacent to the closest vertex
- typename Triangulation<dim, spacedim>::active_cell_iterator
- current_cell = out_particle->get_surrounding_cell(*triangulation);
-
const unsigned int closest_vertex =
GridTools::find_closest_vertex_of_cell<dim, spacedim>(
current_cell, out_particle->get_location(), *mapping);
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(),
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)
// 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;
}
}
// 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
{
}
}
- // Sort the updated particles. This pre-sort speeds up inserting
- // them into particles to O(N) complexity.
- std::multimap<internal::LevelInd, Particle<dim, spacedim>>
- sorted_particles_map;
-
// Exchange particles between processors if we have more than one process
#ifdef DEAL_II_WITH_MPI
if (const auto parallel_triangulation =
{
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
#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();
for (const auto ghost_owner : ghost_owners)
ghost_particles_cache.ghost_particles_by_domain[ghost_owner].reserve(
static_cast<typename std::vector<particle_iterator>::size_type>(
- particles.size() * 0.25));
+ local_number_of_particles * 0.25));
const std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain =
triangulation_cache->get_vertex_to_neighbor_subdomain();
#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(
void
ParticleHandler<dim, spacedim>::send_recv_particles(
const std::map<types::subdomain_id, std::vector<particle_iterator>>
- &particles_to_send,
- std::multimap<internal::LevelInd, Particle<dim, spacedim>>
- &received_particles,
+ & particles_to_send,
+ particle_container &received_particles,
const std::map<
types::subdomain_id,
std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
Assert(
parallel_triangulation,
ExcMessage(
- "This function is only implemented for parallel::TriangulationBase objects."));
+ "This function is only implemented for parallel::TriangulationBase "
+ "objects."));
+
+ if (received_particles.size() == 0)
+ received_particles.resize(triangulation->n_active_cells());
// Determine the communication pattern
const std::set<types::subdomain_id> ghost_owners =
std::vector<unsigned int> send_offsets(n_neighbors, 0);
std::vector<char> send_data;
+ Particle<dim, spacedim> test_particle;
+ test_particle.set_property_pool(*property_pool);
+
const unsigned int individual_particle_data_size =
- Utilities::MPI::max(n_send_particles > 0 ?
- ((begin()->serialized_size_in_bytes() +
- (size_callback ? size_callback() : 0))) :
- 0,
- parallel_triangulation->get_communicator());
+ test_particle.serialized_size_in_bytes() +
+ (size_callback ? size_callback() : 0);
const unsigned int individual_total_particle_data_size =
individual_particle_data_size + cellid_size;
static_cast<std::size_t>(
n_particles_to_send *
individual_total_particle_data_size),
- ExcMessage("Overflow when trying to send particle data"));
+ ExcMessage("Overflow when trying to send particle "
+ "data"));
for (unsigned int j = 0; j < n_particles_to_send; ++j)
{
- // If no target cells are given, use the iterator information
+ // If no target cells are given, use the iterator
+ // information
typename Triangulation<dim, spacedim>::active_cell_iterator
cell;
if (send_cells.size() == 0)
- cell =
- particles_to_send.at(neighbors[i])[j]->get_surrounding_cell(
- *triangulation);
+ cell = particles_to_send.at(neighbors[i])[j]
+ ->get_surrounding_cell();
else
cell = send_cells.at(neighbors[i])[j];
const typename Triangulation<dim, spacedim>::active_cell_iterator cell =
triangulation->create_cell_iterator(id);
- typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator
- recv_particle = received_particles.insert(std::make_pair(
- internal::LevelInd(cell->level(), cell->index()),
- Particle<dim, spacedim>(recv_data_it, property_pool.get())));
+ const unsigned int active_cell_index = cell->active_cell_index();
+ received_particles[active_cell_index].push_back(
+ Particle<dim, spacedim>(recv_data_it, property_pool.get()));
+
+ particle_iterator particle_it(
+ received_particles,
+ cell,
+ received_particles[active_cell_index].size() - 1);
+ particle_it->set_property_pool(*property_pool);
if (load_callback)
- recv_data_it =
- load_callback(particle_iterator(received_particles, recv_particle),
- recv_data_it);
+ recv_data_it = load_callback(particle_it, recv_data_it);
if (build_cache) // TODO: is this safe?
- ghost_particles_iterators.push_back(recv_particle);
+ ghost_particles_iterators.push_back(particle_it);
}
AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
void
ParticleHandler<dim, spacedim>::send_recv_particles_properties_and_location(
const std::map<types::subdomain_id, std::vector<particle_iterator>>
- &particles_to_send,
- std::multimap<internal::LevelInd, Particle<dim, spacedim>>
- &updated_particles)
+ & particles_to_send,
+ particle_container &updated_particles)
{
const auto &neighbors = ghost_particles_cache.neighbors;
const auto &send_pointers = ghost_particles_cache.send_pointers;
Assert(
parallel_triangulation,
ExcMessage(
- "This function is only implemented for parallel::TriangulationBase objects."));
+ "This function is only implemented for parallel::TriangulationBase "
+ "objects."));
std::vector<char> &send_data = ghost_particles_cache.send_data;
// 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(),
// If the cell persist or is refined store all particles of the
// current cell.
{
- unsigned int n_particles = 0;
-
- const internal::LevelInd level_index = {cell->level(),
- cell->index()};
- const auto particles_in_cell =
- (cell->is_ghost() ? ghost_particles.equal_range(level_index) :
- particles.equal_range(level_index));
-
- n_particles = n_particles_in_cell(cell);
- stored_particles_on_cell.reserve(n_particles);
-
- std::for_each(
- particles_in_cell.first,
- particles_in_cell.second,
- [&stored_particles_on_cell](
- const std::pair<internal::LevelInd, Particle<dim, spacedim>>
- &particle) {
- stored_particles_on_cell.push_back(particle.second);
- });
-
- AssertDimension(n_particles, stored_particles_on_cell.size());
+ stored_particles_on_cell = particles[cell->active_cell_index()];
}
break;
// If this cell is the parent of children that will be coarsened,
// collect the particles of all children.
{
- unsigned int n_particles = 0;
-
for (const auto &child : cell->child_iterators())
{
- n_particles += n_particles_in_cell(child);
- }
+ const unsigned int active_cell_index =
+ child->active_cell_index();
+ const unsigned int n_particles = n_particles_in_cell(child);
- stored_particles_on_cell.reserve(n_particles);
+ stored_particles_on_cell.reserve(
+ stored_particles_on_cell.size() + n_particles);
- for (const auto &child : cell->child_iterators())
- {
- const internal::LevelInd level_index = {child->level(),
- child->index()};
- const auto particles_in_cell =
- (child->is_ghost() ?
- ghost_particles.equal_range(level_index) :
- particles.equal_range(level_index));
-
- std::for_each(
- particles_in_cell.first,
- particles_in_cell.second,
- [&stored_particles_on_cell](
- const std::pair<internal::LevelInd, Particle<dim, spacedim>>
- &particle) {
- stored_particles_on_cell.push_back(particle.second);
- });
+ for (unsigned int i = 0; i < n_particles; ++i)
+ stored_particles_on_cell.push_back(
+ particles[active_cell_index][i]);
}
-
- AssertDimension(n_particles, stored_particles_on_cell.size());
}
break;
{
case parallel::TriangulationBase<dim, spacedim>::CELL_PERSIST:
{
- auto position_hint = particles.end();
- for (const auto &particle : loaded_particles_on_cell)
- {
- // Use std::multimap::emplace_hint to speed up insertion of
- // particles.
- position_hint =
- particles.emplace_hint(position_hint,
- std::make_pair(cell->level(),
- cell->index()),
- std::move(particle));
- // Move the hint position forward by one, i.e., for the next
- // particle. The 'hint' position will thus be right after the
- // one just inserted.
- ++position_hint;
- }
+ particles[cell->active_cell_index()].insert(
+ particles[cell->active_cell_index()].end(),
+ loaded_particles_on_cell.begin(),
+ loaded_particles_on_cell.end());
}
break;
case parallel::TriangulationBase<dim, spacedim>::CELL_COARSEN:
{
- typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator
- position_hint = particles.end();
+ const unsigned int active_cell_index = cell->active_cell_index();
+
for (auto &particle : loaded_particles_on_cell)
{
const Point<dim> p_unit =
mapping->transform_real_to_unit_cell(cell,
particle.get_location());
particle.set_reference_location(p_unit);
- // Use std::multimap::emplace_hint to speed up insertion of
- // particles.
- position_hint =
- particles.emplace_hint(position_hint,
- std::make_pair(cell->level(),
- cell->index()),
- std::move(particle));
- // Move the hint position forward by one, i.e., for the next
- // particle. The 'hint' position will thus be right after the
- // one just inserted.
- ++position_hint;
+ particles[active_cell_index].emplace_back(std::move(particle));
}
}
break;
case parallel::TriangulationBase<dim, spacedim>::CELL_REFINE:
{
- std::vector<
- typename std::multimap<internal::LevelInd,
- Particle<dim, spacedim>>::iterator>
- position_hints(GeometryInfo<dim>::max_children_per_cell);
- for (unsigned int child_index = 0;
- child_index < GeometryInfo<dim>::max_children_per_cell;
- ++child_index)
- {
- const typename Triangulation<dim, spacedim>::cell_iterator
- child = cell->child(child_index);
- position_hints[child_index] = particles.upper_bound(
- std::make_pair(child->level(), child->index()));
- }
-
for (auto &particle : loaded_particles_on_cell)
{
for (unsigned int child_index = 0;
++child_index)
{
const typename Triangulation<dim, spacedim>::cell_iterator
- child = cell->child(child_index);
+ child = cell->child(child_index);
+ const unsigned int active_cell_index =
+ child->active_cell_index();
try
{
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;
}
}
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();
auto particle = particle_handler.begin();
while (particle != particle_handler.end())
{
- const auto &cell = particle->get_surrounding_cell(tria);
+ const auto &cell = particle->get_surrounding_cell();
const auto dh_cell =
typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
dh_cell->get_dof_indices(dof_indices);
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();
auto particle = particle_handler.begin();
while (particle != particle_handler.end())
{
- const auto &cell = particle->get_surrounding_cell(tria);
+ const auto &cell = particle->get_surrounding_cell();
const auto &dh_cell =
typename DoFHandler<dim, spacedim>::cell_iterator(*cell, &space_dh);
dh_cell->get_dof_indices(dof_indices);
2,
0);
- auto particle_it = particle_handler.insert_particle(particle1, cell1);
- particle_it->set_properties(properties);
+ if (Utilities::MPI::this_mpi_process(tr.get_communicator()) == 0)
+ {
+ auto particle_it = particle_handler.insert_particle(particle1, cell1);
+ particle_it->set_properties(properties);
+ }
std::vector<std::string> data_names;
std::vector<DataComponentInterpretation::DataComponentInterpretation>
{
Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
- initlog();
+ MPILogInitAll all;
deallog.push("2d/2d");
test<2, 2>();
deallog.pop();
# <x> <y> <id> <property_coord_0> <property_coord_1>
0.125000 0.125000 0.00000 0.125000 0.125000
-DEAL:2d/2d::OK
+DEAL:0:2d/2d::OK
# This file was generated by the deal.II library.
# <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2>
0.125000 0.125000 0.00000 0.00000 0.125000 0.125000 0.125000
-DEAL:2d/3d::OK
+DEAL:0:2d/3d::OK
# This file was generated by the deal.II library.
# <x> <y> <z> <id> <property_coord_0> <property_coord_1> <property_coord_2>
0.125000 0.125000 0.125000 0.00000 0.125000 0.125000 0.125000
-DEAL:3d/3d::OK
+DEAL:0:3d/3d::OK
+
+DEAL:1:2d/2d::OK
+DEAL:1:2d/3d::OK
+DEAL:1:3d/3d::OK
+
tr.refine_global(2);
MappingQ<dim, spacedim> mapping(1);
- // both processes create a particle handler, but only the first creates
- // particles
+ // both processes create a particle handler with a particle in a
+ // position that is in a ghost cell to the other process
Particles::ParticleHandler<dim, spacedim> particle_handler(tr, mapping);
Point<spacedim> position;
position,
reference_position,
Utilities::MPI::this_mpi_process(tr.get_communicator()));
+
+ // We give a random cell hint to check that sorting and
+ // transferring ghost particles works.
typename Triangulation<dim, spacedim>::active_cell_iterator cell =
tr.begin_active();
particle_handler.insert_particle(particle, cell);
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
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
-// Like particle_03, but tests the creation and use of a
+// Like particle_handler_serial_01, but tests the creation and use of a
// particle iterator from the created particle.
#include <deal.II/base/array_view.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
#include <deal.II/particles/particle.h>
#include <deal.II/particles/particle_iterator.h>
test()
{
{
+ Triangulation<dim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ MappingQ<dim> mapping(1);
+
const unsigned int n_properties_per_particle = 3;
Particles::PropertyPool<dim> pool(n_properties_per_particle);
particle.set_properties(
ArrayView<double>(&properties[0], properties.size()));
- std::multimap<Particles::internal::LevelInd, Particles::Particle<dim>> map;
+ std::vector<std::vector<Particles::Particle<dim>>> particle_container;
+ particle_container.resize(1);
- Particles::internal::LevelInd level_index = std::make_pair(0, 0);
- map.insert(std::make_pair(level_index, particle));
+ particle_container[0].push_back(particle);
particle.get_properties()[0] = 0.05;
- map.insert(std::make_pair(level_index, particle));
-
- Particles::ParticleIterator<dim> particle_it(map, map.begin());
- Particles::ParticleIterator<dim> particle_end(map, map.end());
+ particle_container[0].push_back(particle);
+
+ Particles::ParticleIterator<dim> particle_it(particle_container,
+ tr.begin(),
+ 0);
+ Particles::ParticleIterator<dim> particle_end(particle_container,
+ tr.end(),
+ 0);
for (; particle_it != particle_end; ++particle_it)
{
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 - 2020 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+
+// Like particle_iterator_01, but tests state of the iterator.
+
+#include <deal.II/base/array_view.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+
+#include <deal.II/particles/particle.h>
+#include <deal.II/particles/particle_iterator.h>
+
+#include "../tests.h"
+
+
+template <int dim>
+void
+test()
+{
+ {
+ Triangulation<dim> tr;
+
+ GridGenerator::hyper_cube(tr);
+ MappingQ<dim> mapping(1);
+
+ const unsigned int n_properties_per_particle = 3;
+ Particles::PropertyPool<dim> pool(n_properties_per_particle);
+
+ Point<dim> position;
+ position(0) = 0.3;
+
+ if (dim > 1)
+ position(1) = 0.5;
+
+ Point<dim> reference_position;
+ reference_position(0) = 0.2;
+
+ if (dim > 1)
+ reference_position(1) = 0.4;
+
+ const types::particle_index index(7);
+
+ std::vector<double> properties = {0.15, 0.45, 0.75};
+
+ Particles::Particle<dim> particle(position, reference_position, index);
+ particle.set_property_pool(pool);
+ particle.set_properties(
+ ArrayView<double>(&properties[0], properties.size()));
+
+ std::vector<std::vector<Particles::Particle<dim>>> particle_container;
+ particle_container.resize(1);
+
+ particle_container[0].push_back(particle);
+
+ Particles::ParticleIterator<dim> particle_begin(particle_container,
+ tr.begin(),
+ 0);
+ Particles::ParticleIterator<dim> particle_end(particle_container,
+ tr.end(),
+ 0);
+ Particles::ParticleIterator<dim> particle_nonexistent1(particle_container,
+ tr.begin(),
+ 1);
+ Particles::ParticleIterator<dim> particle_nonexistent2(particle_container,
+ tr.end(),
+ 1);
+ Particles::ParticleIterator<dim> particle_invalid;
+
+ Assert(particle_begin->state() == IteratorState::valid, ExcInternalError());
+ Assert(particle_end->state() == IteratorState::past_the_end,
+ ExcInternalError());
+ Assert(particle_nonexistent1->state() == IteratorState::invalid,
+ ExcInternalError());
+ Assert(particle_nonexistent2->state() == IteratorState::invalid,
+ ExcInternalError());
+ Assert(particle_invalid->state() == IteratorState::invalid,
+ ExcInternalError());
+
+ Particles::ParticleIterator<dim> particle_iterator = particle_begin;
+ ++particle_iterator;
+
+ Assert(particle_iterator->state() == IteratorState::past_the_end,
+ ExcInternalError());
+
+ particle_iterator = particle_begin;
+ --particle_iterator;
+
+ Assert(particle_iterator->state() == IteratorState::past_the_end,
+ ExcInternalError());
+
+ for (particle_iterator = particle_begin; particle_iterator != particle_end;
+ ++particle_iterator)
+ {
+ deallog << "Particle position: " << particle_iterator->get_location()
+ << ". Iterator valid: "
+ << (particle_iterator->state() == IteratorState::valid)
+ << std::endl;
+ }
+
+ Assert(particle_iterator->state() == IteratorState::past_the_end,
+ ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int
+main()
+{
+ initlog();
+ test<2>();
+ test<3>();
+}
--- /dev/null
+
+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
-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
projects / dealii.git / commitdiff