From c66e7a60d161e04b11ba15d6b7f2ebde3c34e94e Mon Sep 17 00:00:00 2001 From: Rene Gassmoeller Date: Fri, 10 Nov 2017 17:17:07 -0700 Subject: [PATCH] Hide LevelInd, and add insert_particles(points) --- include/deal.II/particles/particle.h | 14 +-- include/deal.II/particles/particle_accessor.h | 9 +- include/deal.II/particles/particle_handler.h | 33 ++++--- include/deal.II/particles/particle_iterator.h | 9 +- source/particles/particle_accessor.cc | 6 +- source/particles/particle_handler.cc | 89 +++++++++++++++---- source/particles/particle_iterator.cc | 4 +- .../particle_handler_06.mpirun=2.output | 27 +++--- tests/particles/particle_handler_07.cc | 83 +++++++++++++++++ tests/particles/particle_handler_07.output | 7 ++ tests/particles/particle_iterator_01.cc | 4 +- 11 files changed, 223 insertions(+), 62 deletions(-) create mode 100644 tests/particles/particle_handler_07.cc create mode 100644 tests/particles/particle_handler_07.output diff --git a/include/deal.II/particles/particle.h b/include/deal.II/particles/particle.h index b492ecc5a5..0e9b80641e 100644 --- a/include/deal.II/particles/particle.h +++ b/include/deal.II/particles/particle.h @@ -29,12 +29,6 @@ DEAL_II_NAMESPACE_OPEN */ namespace types { - /** - * Typedef of cell level/index pair. TODO: replace this by the - * active_cell_index. - */ - typedef std::pair LevelInd; - /* Type definitions */ #ifdef DEAL_II_WITH_64BIT_INDICES @@ -84,6 +78,14 @@ namespace types */ namespace Particles { + namespace internal + { + /** + * Internal typedef of cell level/index pair. + */ + typedef std::pair LevelInd; + } + /** * Base class of particles - represents a particle with position, * an ID number and a variable number of properties. This class diff --git a/include/deal.II/particles/particle_accessor.h b/include/deal.II/particles/particle_accessor.h index 754a638a98..958846358a 100644 --- a/include/deal.II/particles/particle_accessor.h +++ b/include/deal.II/particles/particle_accessor.h @@ -181,27 +181,28 @@ namespace Particles * This constructor is protected so that it can only be accessed by friend * classes. */ - ParticleAccessor (const std::multimap > &map, - const typename std::multimap >::iterator &particle); + ParticleAccessor (const std::multimap > &map, + const typename std::multimap >::iterator &particle); private: /** * A pointer to the container that stores the particles. Obviously, * this accessor is invalidated if the container changes. */ - std::multimap > *map; + std::multimap > *map; /** * An iterator into the container of particles. Obviously, * this accessor is invalidated if the container changes. */ - typename std::multimap >::iterator particle; + typename std::multimap >::iterator particle; /** * Make ParticleIterator a friend to allow it constructing ParticleAccessors. */ template friend class ParticleIterator; template friend class ParticleHandler; + }; } diff --git a/include/deal.II/particles/particle_handler.h b/include/deal.II/particles/particle_handler.h index 87256ea959..03b3ce2da7 100644 --- a/include/deal.II/particles/particle_handler.h +++ b/include/deal.II/particles/particle_handler.h @@ -185,12 +185,24 @@ namespace Particles const typename Triangulation::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 > &particles); + insert_particles(const std::multimap::active_cell_iterator, + Particle > &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 > &positions); /** * This function allows to register three additional functions that are @@ -227,7 +239,7 @@ namespace Particles /** * 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. * @@ -237,7 +249,7 @@ namespace Particles /** * 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. */ @@ -250,8 +262,9 @@ namespace Particles 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; @@ -333,18 +346,18 @@ 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; + std::multimap > 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 > ghost_particles; + std::multimap > 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; @@ -439,7 +452,7 @@ namespace Particles */ void send_recv_particles(const std::map > &particles_to_send, - std::multimap > &received_particles, + std::multimap > &received_particles, const std::map::active_cell_iterator> > &new_cells_for_particles = std::map::active_cell_iterator> > ()); diff --git a/include/deal.II/particles/particle_iterator.h b/include/deal.II/particles/particle_iterator.h index c8ecf24446..d9aedcb2a5 100644 --- a/include/deal.II/particles/particle_iterator.h +++ b/include/deal.II/particles/particle_iterator.h @@ -16,12 +16,15 @@ #ifndef dealii_particles_particle_iterator_h #define dealii_particles_particle_iterator_h +#include #include DEAL_II_NAMESPACE_OPEN namespace Particles { + template class ParticleHandler; + /** * A class that is used to iterate over particles. Together with the * ParticleAccessor class this is used to hide the internal implementation @@ -38,10 +41,10 @@ namespace Particles /** * 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 > &map, - const typename std::multimap >::iterator &particle); + ParticleIterator (const std::multimap > &map, + const typename std::multimap >::iterator &particle); /** * Dereferencing operator, returns a reference to an accessor. Usage is thus diff --git a/source/particles/particle_accessor.cc b/source/particles/particle_accessor.cc index 92c0026e6b..f7f3eb62b6 100644 --- a/source/particles/particle_accessor.cc +++ b/source/particles/particle_accessor.cc @@ -29,10 +29,10 @@ namespace Particles template - ParticleAccessor::ParticleAccessor (const std::multimap > &map, - const typename std::multimap >::iterator &particle) + ParticleAccessor::ParticleAccessor (const std::multimap > &map, + const typename std::multimap >::iterator &particle) : - map (const_cast > *> (&map)), + map (const_cast > *> (&map)), particle (particle) {} diff --git a/source/particles/particle_handler.cc b/source/particles/particle_handler.cc index 785ed3a380..ae6d927ba4 100644 --- a/source/particles/particle_handler.cc +++ b/source/particles/particle_handler.cc @@ -16,6 +16,8 @@ #include #include +#include + #include DEAL_II_NAMESPACE_OPEN @@ -221,10 +223,10 @@ namespace Particles typename ParticleHandler::particle_iterator_range ParticleHandler::particles_in_cell(const typename Triangulation::active_cell_iterator &cell) { - const types::LevelInd level_index = std::make_pair (cell->level(),cell->index()); + const internal::LevelInd level_index = std::make_pair (cell->level(),cell->index()); - std::pair >::iterator, - typename std::multimap >::iterator> particles_in_cell; + std::pair >::iterator, + typename std::multimap >::iterator> particles_in_cell; if (!cell->is_ghost()) particles_in_cell = particles.equal_range(level_index); @@ -251,8 +253,8 @@ namespace Particles ParticleHandler::insert_particle(const Particle &particle, const typename Triangulation::active_cell_iterator &cell) { - typename std::multimap >::iterator it = - particles.insert(std::make_pair(types::LevelInd(cell->level(),cell->index()),particle)); + typename std::multimap >::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); @@ -268,9 +270,62 @@ namespace Particles template void - ParticleHandler::insert_particles(const std::multimap > &new_particles) + ParticleHandler::insert_particles(const std::multimap::active_cell_iterator, + Particle > &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 + void + ParticleHandler::insert_particles(const std::vector > &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 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; ilevel(),cells[i]->index()); + for (unsigned int p=0; p(positions[index_map[i][p]], + local_positions[i][p], + local_start_index+index_map[i][p]))); + } + } + update_cached_numbers(); } @@ -316,7 +371,7 @@ namespace Particles unsigned int ParticleHandler::n_particles_in_cell(const typename Triangulation::active_cell_iterator &cell) const { - const types::LevelInd found_cell = std::make_pair (cell->level(),cell->index()); + const internal::LevelInd found_cell = std::make_pair (cell->level(),cell->index()); if (cell->is_locally_owned()) return particles.count(found_cell); @@ -420,7 +475,7 @@ namespace Particles // 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 > > sorted_particles; + std::vector > > sorted_particles; std::map > moved_particles; std::map::active_cell_iterator> > moved_cells; @@ -537,7 +592,7 @@ namespace Particles // 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 @@ -550,7 +605,7 @@ namespace Particles // Sort the updated particles. This pre-sort speeds up inserting // them into particles to O(N) complexity. - std::multimap > sorted_particles_map; + std::multimap > 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) @@ -631,7 +686,7 @@ namespace Particles template void ParticleHandler::send_recv_particles(const std::map > &particles_to_send, - std::multimap > &received_particles, + std::multimap > &received_particles, const std::map::active_cell_iterator> > &send_cells) { // Determine the communication pattern @@ -759,8 +814,8 @@ namespace Particles const typename Triangulation::active_cell_iterator cell = id.to_cell(*triangulation); - typename std::multimap >::iterator recv_particle = - received_particles.insert(std::make_pair(types::LevelInd(cell->level(),cell->index()), + typename std::multimap >::iterator recv_particle = + received_particles.insert(std::make_pair(internal::LevelInd(cell->level(),cell->index()), Particle(recv_data_it,*property_pool))); if (load_callback) @@ -979,7 +1034,7 @@ namespace Particles // particle map. if (status == parallel::distributed::Triangulation::CELL_PERSIST) { - typename std::multimap >::iterator position_hint = particles.end(); + typename std::multimap >::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 @@ -1001,7 +1056,7 @@ namespace Particles else if (status == parallel::distributed::Triangulation::CELL_COARSEN) { - typename std::multimap >::iterator position_hint = particles.end(); + typename std::multimap >::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 @@ -1024,7 +1079,7 @@ namespace Particles } else if (status == parallel::distributed::Triangulation::CELL_REFINE) { - std::vector >::iterator > position_hints(GeometryInfo::max_children_per_cell); + std::vector >::iterator > position_hints(GeometryInfo::max_children_per_cell); for (unsigned int child_index=0; child_index::max_children_per_cell; ++child_index) { const typename Triangulation::cell_iterator child = cell->child(child_index); diff --git a/source/particles/particle_iterator.cc b/source/particles/particle_iterator.cc index 67c645f913..90e3071b84 100644 --- a/source/particles/particle_iterator.cc +++ b/source/particles/particle_iterator.cc @@ -20,8 +20,8 @@ DEAL_II_NAMESPACE_OPEN namespace Particles { template - ParticleIterator::ParticleIterator (const std::multimap > &map, - const typename std::multimap >::iterator &particle) + ParticleIterator::ParticleIterator (const std::multimap > &map, + const typename std::multimap >::iterator &particle) : accessor (map, particle) {} diff --git a/tests/particles/particle_handler_06.mpirun=2.output b/tests/particles/particle_handler_06.mpirun=2.output index 01df526a8d..af7cf97de2 100644 --- a/tests/particles/particle_handler_06.mpirun=2.output +++ b/tests/particles/particle_handler_06.mpirun=2.output @@ -1,24 +1,21 @@ -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 diff --git a/tests/particles/particle_handler_07.cc b/tests/particles/particle_handler_07.cc new file mode 100644 index 0000000000..336b4f964a --- /dev/null +++ b/tests/particles/particle_handler_07.cc @@ -0,0 +1,83 @@ +// --------------------------------------------------------------------- +// +// 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 + +#include +#include +#include + +#include + +template +void test () +{ + { + parallel::distributed::Triangulation tr(MPI_COMM_WORLD); + + GridGenerator::hyper_cube(tr); + tr.refine_global(2); + + MappingQ mapping(1); + + Particles::ParticleHandler particle_handler(tr,mapping); + + std::vector > points(10); + for (unsigned int i=0; i<10; ++i) + { + const double coordinate = static_cast (i)/10.0; + for (unsigned int j=0; jget_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(); +} diff --git a/tests/particles/particle_handler_07.output b/tests/particles/particle_handler_07.output new file mode 100644 index 0000000000..cbfa261dee --- /dev/null +++ b/tests/particles/particle_handler_07.output @@ -0,0 +1,7 @@ + +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 diff --git a/tests/particles/particle_iterator_01.cc b/tests/particles/particle_iterator_01.cc index 2eac94b1c6..517d2b2e3f 100644 --- a/tests/particles/particle_iterator_01.cc +++ b/tests/particles/particle_iterator_01.cc @@ -52,9 +52,9 @@ void test () particle.set_property_pool(pool); particle.set_properties(ArrayView(&properties[0],properties.size())); - std::multimap > map; + std::multimap > 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; -- 2.39.5