*/
namespace types
{
- /**
- * Typedef of cell level/index pair. TODO: replace this by the
- * active_cell_index.
- */
- typedef std::pair<int, int> LevelInd;
-
/* Type definitions */
#ifdef DEAL_II_WITH_64BIT_INDICES
*/
namespace Particles
{
+ namespace internal
+ {
+ /**
+ * Internal typedef of cell level/index pair.
+ */
+ typedef std::pair<int, int> LevelInd;
+ }
+
/**
* Base class of particles - represents a particle with position,
* an ID number and a variable number of properties. This class
* This constructor is protected so that it can only be accessed by friend
* classes.
*/
- ParticleAccessor (const std::multimap<types::LevelInd, Particle<dim,spacedim> > &map,
- const typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator &particle);
+ ParticleAccessor (const std::multimap<internal::LevelInd, Particle<dim,spacedim> > &map,
+ const typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator &particle);
private:
/**
* A pointer to the container that stores the particles. Obviously,
* this accessor is invalidated if the container changes.
*/
- std::multimap<types::LevelInd, Particle<dim,spacedim> > *map;
+ std::multimap<internal::LevelInd, Particle<dim,spacedim> > *map;
/**
* An iterator into the container of particles. Obviously,
* this accessor is invalidated if the container changes.
*/
- typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator particle;
+ typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator particle;
/**
* Make ParticleIterator a friend to allow it constructing ParticleAccessors.
*/
template <int, int> friend class ParticleIterator;
template <int, int> friend class ParticleHandler;
+
};
}
const typename Triangulation<dim,spacedim>::active_cell_iterator &cell);
/**
- * Insert a number of particle into the collection of particles.
+ * Insert a number of particles into the collection of particles.
* This function involves a copy of the particles and their properties.
* Note that this function is of O(n_existing_particles + n_particles) complexity.
*/
void
- insert_particles(const std::multimap<types::LevelInd, Particle<dim,spacedim> > &particles);
+ insert_particles(const std::multimap<typename Triangulation<dim,spacedim>::active_cell_iterator,
+ Particle<dim,spacedim> > &particles);
+
+ /**
+ * Insert a number of particles into the collection of particles.
+ * This function takes a list of positions and creates a set of particles
+ * at these positions, which are then added to the local particle
+ * collection. Note that this function assumes all positions are within
+ * the local part of the triangulation, if one of them is not in the
+ * local domain this function will throw an exception.
+ */
+ void
+ insert_particles(const std::vector<Point<spacedim> > &positions);
/**
* This function allows to register three additional functions that are
/**
* Return the total number of particles that were managed by this class
- * the last time the update_n_global_particles() function was called.
+ * the last time the update_cached_numbers() function was called.
* The actual number of particles may have changed since then if
* particles have been added or removed.
*
/**
* Return the maximum number of particles per cell the last
- * time the update_n_global_particles() function was called.
+ * time the update_cached_numbers() function was called.
*
* @return Maximum number of particles in one cell in simulation.
*/
types::particle_index n_locally_owned_particles() const;
/**
- * Return the number of particles in the local part of the
- * triangulation.
+ * Return the next free particle index in the global set
+ * of particles the last
+ * time the update_cached_numbers() function was called.
*/
types::particle_index get_next_free_particle_index() const;
* Set of particles currently living in the local domain, organized by
* the level/index of the cell they are in.
*/
- std::multimap<types::LevelInd, Particle<dim,spacedim> > particles;
+ std::multimap<internal::LevelInd, Particle<dim,spacedim> > 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.
*/
- std::multimap<types::LevelInd, Particle<dim,spacedim> > ghost_particles;
+ std::multimap<internal::LevelInd, Particle<dim,spacedim> > ghost_particles;
/**
* This variable stores how many particles are stored globally. It is
- * calculated by update_n_global_particles().
+ * calculated by update_cached_numbers().
*/
types::particle_index global_number_of_particles;
*/
void
send_recv_particles(const std::map<types::subdomain_id, std::vector<particle_iterator> > &particles_to_send,
- std::multimap<types::LevelInd,Particle <dim,spacedim> > &received_particles,
+ std::multimap<internal::LevelInd,Particle <dim,spacedim> > &received_particles,
const std::map<types::subdomain_id, std::vector<typename Triangulation<dim,spacedim>::active_cell_iterator> > &new_cells_for_particles =
std::map<types::subdomain_id, std::vector<typename Triangulation<dim,spacedim>::active_cell_iterator> > ());
#ifndef dealii_particles_particle_iterator_h
#define dealii_particles_particle_iterator_h
+#include <deal.II/base/config.h>
#include <deal.II/particles/particle_accessor.h>
DEAL_II_NAMESPACE_OPEN
namespace Particles
{
+ template <int, int> class ParticleHandler;
+
/**
* A class that is used to iterate over particles. Together with the
* ParticleAccessor class this is used to hide the internal implementation
/**
* Constructor of the iterator. Takes a reference to the particle
- * container, and an iterator the the cell-particle pair.
+ * container, and an iterator to the cell-particle pair.
*/
- ParticleIterator (const std::multimap<types::LevelInd, Particle<dim,spacedim> > &map,
- const typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator &particle);
+ ParticleIterator (const std::multimap<internal::LevelInd, Particle<dim,spacedim> > &map,
+ const typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator &particle);
/**
* Dereferencing operator, returns a reference to an accessor. Usage is thus
template <int dim, int spacedim>
- ParticleAccessor<dim,spacedim>::ParticleAccessor (const std::multimap<types::LevelInd, Particle<dim,spacedim> > &map,
- const typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator &particle)
+ 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<types::LevelInd, Particle<dim,spacedim> > *> (&map)),
+ map (const_cast<std::multimap<internal::LevelInd, Particle<dim,spacedim> > *> (&map)),
particle (particle)
{}
#include <deal.II/particles/particle_handler.h>
#include <deal.II/grid/grid_tools.h>
+#include <deal.II/grid/grid_tools_cache.h>
+
#include <utility>
DEAL_II_NAMESPACE_OPEN
typename ParticleHandler<dim,spacedim>::particle_iterator_range
ParticleHandler<dim,spacedim>::particles_in_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell)
{
- const types::LevelInd level_index = std::make_pair<int, int> (cell->level(),cell->index());
+ const internal::LevelInd level_index = std::make_pair<int, int> (cell->level(),cell->index());
- std::pair<typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator,
- typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator> particles_in_cell;
+ std::pair<typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator,
+ typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator> particles_in_cell;
if (!cell->is_ghost())
particles_in_cell = particles.equal_range(level_index);
ParticleHandler<dim,spacedim>::insert_particle(const Particle<dim,spacedim> &particle,
const typename Triangulation<dim,spacedim>::active_cell_iterator &cell)
{
- typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator it =
- particles.insert(std::make_pair(types::LevelInd(cell->level(),cell->index()),particle));
+ typename std::multimap<internal::LevelInd, Particle<dim,spacedim> >::iterator it =
+ particles.insert(std::make_pair(internal::LevelInd(cell->level(),cell->index()),particle));
particle_iterator particle_it (particles,it);
particle_it->set_property_pool(*property_pool);
template <int dim,int spacedim>
void
- ParticleHandler<dim,spacedim>::insert_particles(const std::multimap<types::LevelInd, Particle<dim,spacedim> > &new_particles)
+ ParticleHandler<dim,spacedim>::insert_particles(const std::multimap<typename Triangulation<dim,spacedim>::active_cell_iterator,
+ Particle<dim,spacedim> > &new_particles)
+ {
+ for (auto particle = new_particles.begin(); particle != new_particles.end(); ++particle)
+ particles.insert(particles.end(),
+ std::make_pair(internal::LevelInd(particle->first->level(),particle->first->index()),
+ particle->second));
+
+ update_cached_numbers();
+ }
+
+
+
+ template <int dim,int spacedim>
+ void
+ ParticleHandler<dim,spacedim>::insert_particles(const std::vector<Point<spacedim> > &positions)
{
- particles.insert(new_particles.begin(),new_particles.end());
+ update_cached_numbers();
+
+ const types::particle_index local_next_particle_index = get_next_free_particle_index();
+ const types::particle_index particles_to_add_locally = positions.size();
+
+ // Determine the starting particle index of this process, which
+ // is the highest currently existing particle index plus the sum
+ // of the number of newly generated particles of all
+ // processes with a lower rank.
+ types::particle_index local_start_index = 0;
+ MPI_Scan(&particles_to_add_locally, &local_start_index, 1, PARTICLE_INDEX_MPI_TYPE, MPI_SUM, triangulation->get_communicator());
+
+ local_start_index -= particles_to_add_locally;
+ local_start_index += local_next_particle_index;
+
+ GridTools::Cache<dim,spacedim> cache(*triangulation, *mapping);
+ auto point_locations = GridTools::compute_point_locations(cache, positions);
+
+ auto &cells = std::get<0>(point_locations);
+ auto &local_positions = std::get<1>(point_locations);
+ auto &index_map = std::get<2>(point_locations);
+
+ 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());
+ 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])));
+ }
+ }
+
update_cached_numbers();
}
unsigned int
ParticleHandler<dim,spacedim>::n_particles_in_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell) const
{
- const types::LevelInd found_cell = std::make_pair<int, int> (cell->level(),cell->index());
+ const internal::LevelInd found_cell = std::make_pair<int, int> (cell->level(),cell->index());
if (cell->is_locally_owned())
return particles.count(found_cell);
// sorted_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<types::LevelInd, Particle<dim,spacedim> > > sorted_particles;
+ std::vector<std::pair<internal::LevelInd, Particle<dim,spacedim> > > sorted_particles;
std::map<types::subdomain_id, std::vector<particle_iterator> > moved_particles;
std::map<types::subdomain_id, std::vector<typename Triangulation<dim,spacedim>::active_cell_iterator> > moved_cells;
// Mark it for MPI transfer otherwise
if (current_cell->is_locally_owned())
{
- sorted_particles.push_back(std::make_pair(types::LevelInd(current_cell->level(),current_cell->index()),
+ sorted_particles.push_back(std::make_pair(internal::LevelInd(current_cell->level(),current_cell->index()),
(*it)->particle->second));
}
else
// Sort the updated particles. This pre-sort speeds up inserting
// them into particles to O(N) complexity.
- std::multimap<types::LevelInd,Particle <dim,spacedim> > sorted_particles_map;
+ std::multimap<internal::LevelInd,Particle <dim,spacedim> > sorted_particles_map;
// Exchange particles between processors if we have more than one process
if (dealii::Utilities::MPI::n_mpi_processes(triangulation->get_communicator()) > 1)
template <int dim, int spacedim>
void
ParticleHandler<dim,spacedim>::send_recv_particles(const std::map<types::subdomain_id, std::vector<particle_iterator> > &particles_to_send,
- std::multimap<types::LevelInd,Particle <dim,spacedim> > &received_particles,
+ std::multimap<internal::LevelInd,Particle <dim,spacedim> > &received_particles,
const std::map<types::subdomain_id, std::vector<typename Triangulation<dim,spacedim>::active_cell_iterator> > &send_cells)
{
// Determine the communication pattern
const typename Triangulation<dim,spacedim>::active_cell_iterator cell = id.to_cell(*triangulation);
- typename std::multimap<types::LevelInd,Particle <dim,spacedim> >::iterator recv_particle =
- received_particles.insert(std::make_pair(types::LevelInd(cell->level(),cell->index()),
+ 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)));
if (load_callback)
// particle map.
if (status == parallel::distributed::Triangulation<dim,spacedim>::CELL_PERSIST)
{
- typename std::multimap<types::LevelInd,Particle<dim,spacedim> >::iterator position_hint = particles.end();
+ typename std::multimap<internal::LevelInd,Particle<dim,spacedim> >::iterator position_hint = particles.end();
for (unsigned int i = 0; i < *n_particles_in_cell_ptr; ++i)
{
// Use std::multimap::emplace_hint to speed up insertion of
else if (status == parallel::distributed::Triangulation<dim,spacedim>::CELL_COARSEN)
{
- typename std::multimap<types::LevelInd,Particle<dim,spacedim> >::iterator position_hint = particles.end();
+ typename std::multimap<internal::LevelInd,Particle<dim,spacedim> >::iterator position_hint = particles.end();
for (unsigned int i = 0; i < *n_particles_in_cell_ptr; ++i)
{
// Use std::multimap::emplace_hint to speed up insertion of
}
else if (status == parallel::distributed::Triangulation<dim,spacedim>::CELL_REFINE)
{
- std::vector<typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator > position_hints(GeometryInfo<dim>::max_children_per_cell);
+ 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);
namespace Particles
{
template <int dim, int spacedim>
- ParticleIterator<dim,spacedim>::ParticleIterator (const std::multimap<types::LevelInd, Particle<dim,spacedim> > &map,
- const typename std::multimap<types::LevelInd, Particle<dim,spacedim> >::iterator &particle)
+ 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)
{}
-DEAL:0:2d/2d::Before sort particle id 0 is in cell 2.0 on process 0
-DEAL:0:2d/2d::Before sort particle id 1 is in cell 2.0 on process 0
-DEAL:0:2d/2d::After sort particle id 0 is in cell 2.0 on process 0
+DEAL:0:2d/2d::Particle id 0 is local particle on process 0
+DEAL:0:2d/2d::Particle id 1 is ghost particle on process 0
DEAL:0:2d/2d::OK
-DEAL:0:2d/3d::Before sort particle id 0 is in cell 2.0 on process 0
-DEAL:0:2d/3d::Before sort particle id 1 is in cell 2.0 on process 0
-DEAL:0:2d/3d::After sort particle id 0 is in cell 2.0 on process 0
+DEAL:0:2d/3d::Particle id 0 is local particle on process 0
+DEAL:0:2d/3d::Particle id 1 is ghost particle on process 0
DEAL:0:2d/3d::OK
-DEAL:0:3d/3d::Before sort particle id 0 is in cell 2.0 on process 0
-DEAL:0:3d/3d::Before sort particle id 1 is in cell 2.0 on process 0
-DEAL:0:3d/3d::After sort particle id 0 is in cell 2.0 on process 0
+DEAL:0:3d/3d::Particle id 0 is local particle on process 0
+DEAL:0:3d/3d::Particle id 1 is ghost particle on process 0
DEAL:0:3d/3d::OK
-DEAL:1:2d/2d::After sort particle id 1 is in cell 2.12 on process 1
-DEAL:1:2d/2d::After shift particle id 0 is in cell 2.8 on process 1
+DEAL:1:2d/2d::Particle id 1 is local particle on process 1
+DEAL:1:2d/2d::Particle id 0 is ghost particle on process 1
DEAL:1:2d/2d::OK
-DEAL:1:2d/3d::After sort particle id 1 is in cell 2.12 on process 1
-DEAL:1:2d/3d::After shift particle id 0 is in cell 2.8 on process 1
+DEAL:1:2d/3d::Particle id 1 is local particle on process 1
+DEAL:1:2d/3d::Particle id 0 is ghost particle on process 1
DEAL:1:2d/3d::OK
-DEAL:1:3d/3d::After sort particle id 1 is in cell 2.56 on process 1
-DEAL:1:3d/3d::After shift particle id 0 is in cell 2.32 on process 1
+DEAL:1:3d/3d::Particle id 1 is local particle on process 1
+DEAL:1:3d/3d::Particle id 0 is ghost particle on process 1
DEAL:1:3d/3d::OK
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2017 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 at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+// check the creation and destruction of particle within the particle handler class.
+
+#include "../tests.h"
+#include <deal.II/particles/particle_handler.h>
+
+#include <deal.II/distributed/tria.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_tools.h>
+
+#include <deal.II/fe/mapping_q.h>
+
+template <int dim, int spacedim>
+void test ()
+{
+ {
+ parallel::distributed::Triangulation<dim,spacedim> tr(MPI_COMM_WORLD);
+
+ GridGenerator::hyper_cube(tr);
+ tr.refine_global(2);
+
+ MappingQ<dim,spacedim> mapping(1);
+
+ Particles::ParticleHandler<dim,spacedim> particle_handler(tr,mapping);
+
+ std::vector<Point<spacedim> > points(10);
+ for (unsigned int i=0; i<10; ++i)
+ {
+ const double coordinate = static_cast<double> (i)/10.0;
+ for (unsigned int j=0; j<spacedim; ++j)
+ points[i][j] = coordinate;
+ }
+
+ particle_handler.insert_particles(points);
+ particle_handler.update_cached_numbers();
+
+ deallog << "Particle number: " << particle_handler.n_global_particles() << std::endl;
+
+ for (auto particle = particle_handler.begin(); particle != particle_handler.end(); ++particle)
+ deallog << "Particle id " << particle->get_id()
+ << " is in cell " << particle->get_surrounding_cell(tr) << std::endl
+ << " at location " << particle->get_location() << std::endl
+ << " at reference location " << particle->get_reference_location()
+ << std::endl;
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+
+
+int main (int argc, char *argv[])
+{
+ Utilities::MPI::MPI_InitFinalize mpi_initialization (argc, argv, 1);
+
+ initlog();
+
+ deallog.push("2d/2d");
+ test<2,2>();
+ deallog.pop();
+ deallog.push("2d/3d");
+ test<2,3>();
+ deallog.pop();
+ deallog.push("3d/3d");
+ test<3,3>();
+ deallog.pop();
+}
--- /dev/null
+
+DEAL::Particle location: 0.00000 0.00000
+DEAL::Particle location: 1.00000 0.00000
+DEAL::OK
+DEAL::Particle location: 0.00000 0.00000 0.00000
+DEAL::Particle location: 1.00000 0.00000 0.00000
+DEAL::OK
particle.set_property_pool(pool);
particle.set_properties(ArrayView<double>(&properties[0],properties.size()));
- std::multimap<types::LevelInd,Particles::Particle<dim> > map;
+ std::multimap<Particles::internal::LevelInd,Particles::Particle<dim> > map;
- types::LevelInd level_index = std::make_pair(0,0);
+ Particles::internal::LevelInd level_index = std::make_pair(0,0);
map.insert(std::make_pair(level_index,particle));
particle.get_properties()[0] = 0.05;