16 #include <deal.II/base/std_cxx14/memory.h> 18 #include <deal.II/grid/grid_tools.h> 19 #include <deal.II/grid/grid_tools_cache.h> 21 #include <deal.II/particles/particle_handler.h> 25 DEAL_II_NAMESPACE_OPEN
27 #ifdef DEAL_II_WITH_P4EST 31 template <
int dim,
int spacedim>
36 , global_number_of_particles(0)
37 , global_max_particles_per_cell(0)
38 , next_free_particle_index(0)
43 , handle(
numbers::invalid_unsigned_int)
48 template <
int dim,
int spacedim>
52 const unsigned int n_properties)
53 : triangulation(&triangulation, typeid(*this).name())
54 , mapping(&mapping, typeid(*this).name())
69 template <
int dim,
int spacedim>
75 template <
int dim,
int spacedim>
80 const unsigned int n_properties)
86 property_pool = std_cxx14::make_unique<PropertyPool>(n_properties);
91 template <
int dim,
int spacedim>
103 template <
int dim,
int spacedim>
112 template <
int dim,
int spacedim>
117 unsigned int local_max_particles_per_cell = 0;
118 unsigned int current_particles_per_cell = 0;
124 locally_highest_index =
125 std::max(locally_highest_index, particle->get_id());
127 if (particle->get_surrounding_cell(*
triangulation) != current_cell)
129 current_particles_per_cell = 0;
130 current_cell = particle->get_surrounding_cell(*
triangulation);
133 ++current_particles_per_cell;
134 local_max_particles_per_cell =
135 std::max(local_max_particles_per_cell, current_particles_per_cell);
142 ::Utilities::MPI::max(locally_highest_index,
146 ::Utilities::MPI::max(local_max_particles_per_cell,
152 template <
int dim,
int spacedim>
161 template <
int dim,
int spacedim>
170 template <
int dim,
int spacedim>
179 template <
int dim,
int spacedim>
188 template <
int dim,
int spacedim>
197 template <
int dim,
int spacedim>
206 template <
int dim,
int spacedim>
215 template <
int dim,
int spacedim>
224 template <
int dim,
int spacedim>
231 ->particles_in_cell(cell);
236 template <
int dim,
int spacedim>
241 const internal::LevelInd level_index =
242 std::make_pair<int, int>(cell->level(), cell->index());
244 if (cell->is_ghost())
247 return boost::make_iterator_range(
253 return boost::make_iterator_range(
260 template <
int dim,
int spacedim>
270 template <
int dim,
int spacedim>
276 typename std::multimap<internal::LevelInd,
279 std::make_pair(internal::LevelInd(cell->level(), cell->index()),
286 for (
unsigned int n = 0; n < particle.
get_properties().size(); ++n)
294 template <
int dim,
int spacedim>
301 for (
auto particle = new_particles.begin(); particle != new_particles.end();
305 std::make_pair(internal::LevelInd(particle->first->level(),
306 particle->first->index()),
314 template <
int dim,
int spacedim>
329 # ifdef DEAL_II_WITH_MPI 331 const int ierr = MPI_Scan(&particles_to_add_locally,
334 PARTICLE_INDEX_MPI_TYPE,
338 local_start_index -= particles_to_add_locally;
341 local_start_index += local_next_particle_index;
346 auto &cells = std::get<0>(point_locations);
347 auto &local_positions = std::get<1>(point_locations);
348 auto &index_map = std::get<2>(point_locations);
350 if (cells.size() == 0)
354 particles.find(std::make_pair(cells[0]->level(), cells[0]->index()));
355 for (
unsigned int i = 0; i < cells.size(); ++i)
357 internal::LevelInd current_cell(cells[i]->level(), cells[i]->index());
358 for (
unsigned int p = 0; p < local_positions[i].size(); ++p)
362 std::make_pair(current_cell,
364 local_positions[i][p],
375 template <
int dim,
int spacedim>
384 template <
int dim,
int spacedim>
393 template <
int dim,
int spacedim>
402 template <
int dim,
int spacedim>
411 template <
int dim,
int spacedim>
417 const internal::LevelInd found_cell =
418 std::make_pair<int, int>(cell->level(), cell->index());
420 if (cell->is_locally_owned())
422 else if (cell->is_ghost())
424 else if (cell->is_artificial())
432 template <
int dim,
int spacedim>
441 template <
int dim,
int spacedim>
463 compare_particle_association(
464 const unsigned int a,
465 const unsigned int b,
469 const double scalar_product_a = center_directions[a] * particle_direction;
470 const double scalar_product_b = center_directions[b] * particle_direction;
475 return (scalar_product_a > scalar_product_b);
481 template <
int dim,
int spacedim>
493 std::vector<particle_iterator> particles_out_of_cell;
505 mapping->transform_real_to_unit_cell(cell, it->get_location());
508 it->set_reference_location(p_unit);
513 particles_out_of_cell.push_back(it);
519 particles_out_of_cell.push_back(it);
529 std::vector<std::pair<internal::LevelInd, Particle<dim, spacedim>>>
531 std::map<types::subdomain_id, std::vector<particle_iterator>>
535 std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
543 using vector_size =
typename std::vector<particle_iterator>::size_type;
544 sorted_particles.reserve(
545 static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
547 const std::set<types::subdomain_id> ghost_owners =
550 for (
auto ghost_domain_id = ghost_owners.begin();
551 ghost_domain_id != ghost_owners.end();
553 moved_particles[*ghost_domain_id].reserve(
554 static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
555 for (
auto ghost_domain_id = ghost_owners.begin();
556 ghost_domain_id != ghost_owners.end();
558 moved_cells[*ghost_domain_id].reserve(
559 static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
564 std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
568 const std::vector<std::vector<Tensor<1, spacedim>>>
569 vertex_to_cell_centers(
573 std::vector<unsigned int> neighbor_permutation;
576 typename std::vector<particle_iterator>::iterator
577 it = particles_out_of_cell.begin(),
578 end_particle = particles_out_of_cell.end();
580 for (; it != end_particle; ++it)
584 bool found_cell =
false;
591 const unsigned int closest_vertex =
592 GridTools::find_closest_vertex_of_cell<dim, spacedim>(
593 current_cell, (*it)->get_location());
595 (*it)->get_location() - current_cell->vertex(closest_vertex);
596 vertex_to_particle /= vertex_to_particle.
norm();
598 const unsigned int closest_vertex_index =
599 current_cell->vertex_index(closest_vertex);
600 const unsigned int n_neighbor_cells =
601 vertex_to_cells[closest_vertex_index].size();
603 neighbor_permutation.resize(n_neighbor_cells);
604 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
605 neighbor_permutation[i] = i;
607 std::sort(neighbor_permutation.begin(),
608 neighbor_permutation.end(),
609 std::bind(&compare_particle_association<spacedim>,
610 std::placeholders::_1,
611 std::placeholders::_2,
612 std::cref(vertex_to_particle),
614 vertex_to_cell_centers[closest_vertex_index])));
618 for (
unsigned int i = 0; i < n_neighbor_cells; ++i)
622 typename std::set<typename Triangulation<dim, spacedim>::
623 active_cell_iterator>::const_iterator
624 cell = vertex_to_cells[closest_vertex_index].begin();
626 std::advance(cell, neighbor_permutation[i]);
628 mapping->transform_real_to_unit_cell(*cell,
629 (*it)->get_location());
632 current_cell = *cell;
633 current_reference_position = p_unit;
649 const std::pair<const typename Triangulation<dim, spacedim>::
650 active_cell_iterator,
652 current_cell_and_position =
653 GridTools::find_active_cell_around_point<>(
655 current_cell = current_cell_and_position.first;
656 current_reference_position = current_cell_and_position.second;
658 catch (GridTools::ExcPointNotFound<spacedim> &)
668 (*it)->set_reference_location(current_reference_position);
672 if (current_cell->is_locally_owned())
674 sorted_particles.push_back(
675 std::make_pair(internal::LevelInd(current_cell->level(),
676 current_cell->index()),
677 (*it)->particle->second));
681 moved_particles[current_cell->subdomain_id()].push_back(*it);
682 moved_cells[current_cell->subdomain_id()].push_back(current_cell);
689 std::multimap<internal::LevelInd, Particle<dim, spacedim>>
690 sorted_particles_map;
693 # ifdef DEAL_II_WITH_MPI 699 sorted_particles_map.insert(sorted_particles.begin(),
700 sorted_particles.end());
702 for (
unsigned int i = 0; i < particles_out_of_cell.size(); ++i)
705 particles.insert(sorted_particles_map.begin(), sorted_particles_map.end());
711 template <
int dim,
int spacedim>
720 # ifdef DEAL_II_WITH_MPI 724 std::map<types::subdomain_id, std::vector<particle_iterator>>
725 ghost_particles_by_domain;
727 const std::set<types::subdomain_id> ghost_owners =
729 for (
auto ghost_domain_id = ghost_owners.begin();
730 ghost_domain_id != ghost_owners.end();
732 ghost_particles_by_domain[*ghost_domain_id].reserve(
733 static_cast<typename std::vector<particle_iterator>::size_type
>(
736 std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain(
742 for (; cell != endc; ++cell)
744 if (cell->is_ghost())
745 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
747 vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(
748 cell->subdomain_id());
752 for (; cell != endc; ++cell)
754 if (!cell->is_ghost())
756 std::set<unsigned int> cell_to_neighbor_subdomain;
757 for (
unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
760 cell_to_neighbor_subdomain.insert(
761 vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
762 vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
765 if (cell_to_neighbor_subdomain.size() > 0)
770 for (std::set<types::subdomain_id>::const_iterator domain =
771 cell_to_neighbor_subdomain.begin();
772 domain != cell_to_neighbor_subdomain.end();
775 for (
typename particle_iterator_range::iterator particle =
776 particle_range.begin();
777 particle != particle_range.end();
779 ghost_particles_by_domain[*domain].push_back(particle);
791 # ifdef DEAL_II_WITH_MPI 792 template <
int dim,
int spacedim>
805 const std::set<types::subdomain_id> ghost_owners =
807 const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
809 const unsigned int n_neighbors = neighbors.size();
811 if (send_cells.size() != 0)
817 particles_to_send.end(),
822 for (
auto send_particles = particles_to_send.begin();
823 send_particles != particles_to_send.end();
825 Assert(ghost_owners.find(send_particles->first) != ghost_owners.end(),
828 unsigned int n_send_particles = 0;
829 for (
auto send_particles = particles_to_send.begin();
830 send_particles != particles_to_send.end();
832 n_send_particles += send_particles->second.size();
838 std::vector<unsigned int> n_send_data(n_neighbors, 0);
839 std::vector<unsigned int> send_offsets(n_neighbors, 0);
840 std::vector<char> send_data;
845 if (n_send_particles > 0)
848 const unsigned int particle_size =
849 begin()->serialized_size_in_bytes() + cellid_size +
851 send_data.resize(n_send_particles * particle_size);
852 void *data =
static_cast<void *
>(&send_data.front());
855 for (
unsigned int i = 0; i < n_neighbors; ++i)
857 send_offsets[i] =
reinterpret_cast<std::size_t
>(data) -
858 reinterpret_cast<std::size_t>(&send_data.front());
860 for (
unsigned int j = 0;
861 j < particles_to_send.at(neighbors[i]).size();
867 if (send_cells.size() == 0)
869 particles_to_send.at(neighbors[i])[j]->get_surrounding_cell(
872 cell = send_cells.at(neighbors[i])[j];
875 cell->id().template to_binary<dim>();
876 memcpy(data, &cellid, cellid_size);
877 data =
static_cast<char *
>(data) + cellid_size;
879 particles_to_send.at(neighbors[i])[j]->write_data(data);
884 n_send_data[i] =
reinterpret_cast<std::size_t
>(data) -
886 reinterpret_cast<std::size_t>(&send_data.front());
891 std::vector<unsigned int> n_recv_data(n_neighbors);
892 std::vector<unsigned int> recv_offsets(n_neighbors);
896 std::vector<MPI_Request> n_requests(2 * n_neighbors);
897 for (
unsigned int i = 0; i < n_neighbors; ++i)
899 const int ierr = MPI_Irecv(&(n_recv_data[i]),
905 &(n_requests[2 * i]));
908 for (
unsigned int i = 0; i < n_neighbors; ++i)
910 const int ierr = MPI_Isend(&(n_send_data[i]),
916 &(n_requests[2 * i + 1]));
920 MPI_Waitall(2 * n_neighbors, &n_requests[0], MPI_STATUSES_IGNORE);
925 unsigned int total_recv_data = 0;
926 for (
unsigned int neighbor_id = 0; neighbor_id < n_neighbors; ++neighbor_id)
928 recv_offsets[neighbor_id] = total_recv_data;
929 total_recv_data += n_recv_data[neighbor_id];
933 std::vector<char> recv_data(total_recv_data);
937 std::vector<MPI_Request> requests(2 * n_neighbors);
938 unsigned int send_ops = 0;
939 unsigned int recv_ops = 0;
941 for (
unsigned int i = 0; i < n_neighbors; ++i)
942 if (n_recv_data[i] > 0)
944 const int ierr = MPI_Irecv(&(recv_data[recv_offsets[i]]),
950 &(requests[send_ops]));
955 for (
unsigned int i = 0; i < n_neighbors; ++i)
956 if (n_send_data[i] > 0)
958 const int ierr = MPI_Isend(&(send_data[send_offsets[i]]),
964 &(requests[send_ops + recv_ops]));
969 MPI_Waitall(send_ops + recv_ops, &requests[0], MPI_STATUSES_IGNORE);
975 const void *recv_data_it =
static_cast<const void *
>(recv_data.data());
977 while (reinterpret_cast<std::size_t>(recv_data_it) -
978 reinterpret_cast<std::size_t
>(recv_data.data()) <
982 memcpy(&binary_cellid, recv_data_it, cellid_size);
983 const CellId id(binary_cellid);
984 recv_data_it =
static_cast<const char *
>(recv_data_it) + cellid_size;
989 typename std::multimap<internal::LevelInd,
991 recv_particle = received_particles.insert(std::make_pair(
992 internal::LevelInd(cell->level(), cell->index()),
1001 AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
1003 "The amount of data that was read into new particles " 1004 "does not match the amount of data sent around."));
1010 template <
int dim,
int spacedim>
1013 const std::function<std::size_t()> & size_callb,
1025 template <
int dim,
int spacedim>
1030 *non_const_triangulation =
1040 const std::function<std::vector<char>(
1046 std::placeholders::_1,
1047 std::placeholders::_2);
1050 callback_function,
false);
1056 template <
int dim,
int spacedim>
1059 const bool serialization)
1066 *non_const_triangulation =
1076 const std::function<std::vector<char>(
1082 std::placeholders::_1,
1083 std::placeholders::_2);
1086 callback_function,
false);
1092 const std::function<void(
1095 const boost::iterator_range<std::vector<char>::const_iterator> &)>
1099 std::placeholders::_1,
1100 std::placeholders::_2,
1101 std::placeholders::_3);
1115 template <
int dim,
int spacedim>
1129 const std::size_t size_per_particle =
1131 begin()->serialized_size_in_bytes() :
1141 const bool serialization =
false;
1142 const std::size_t transfer_size_per_cell =
1143 sizeof(
unsigned int) +
1145 (serialization ? 1 : std::pow(2, dim));
1147 std::vector<char> data_vector(transfer_size_per_cell);
1148 void * data =
static_cast<void *
>(data_vector.data());
1151 unsigned int n_particles(0);
1160 const boost::iterator_range<particle_iterator> particle_range =
1163 std::distance(particle_range.begin(), particle_range.end());
1165 unsigned int *ndata =
static_cast<unsigned int *
>(data);
1166 *ndata = n_particles;
1167 data =
static_cast<void *
>(ndata + 1);
1170 particle != particle_range.end();
1173 particle->write_data(data);
1181 for (
unsigned int child_index = 0;
1182 child_index < GeometryInfo<dim>::max_children_per_cell;
1186 cell->child(child_index);
1190 unsigned int *ndata =
static_cast<unsigned int *
>(data);
1191 *ndata = n_particles;
1193 data =
static_cast<void *
>(ndata + 1);
1195 for (
unsigned int child_index = 0;
1196 child_index < GeometryInfo<dim>::max_children_per_cell;
1200 cell->child(child_index);
1201 const boost::iterator_range<particle_iterator> particle_range =
1205 particle != particle_range.end();
1208 particle->write_data(data);
1218 template <
int dim,
int spacedim>
1223 const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
1229 const void *data =
static_cast<const void *
>(&*(data_range.begin()));
1232 const unsigned int *n_particles_in_cell_ptr =
1233 static_cast<const unsigned int *
>(data);
1235 reinterpret_cast<const void *
>(n_particles_in_cell_ptr + 1);
1237 if (*n_particles_in_cell_ptr == 0)
1245 typename std::multimap<internal::LevelInd,
1248 for (
unsigned int i = 0; i < *n_particles_in_cell_ptr; ++i)
1254 # ifdef DEAL_II_WITH_CXX14 1257 std::make_pair(cell->level(), cell->index()),
1262 std::make_pair(std::make_pair(cell->level(), cell->index()),
1273 typename std::multimap<internal::LevelInd,
1276 for (
unsigned int i = 0; i < *n_particles_in_cell_ptr; ++i)
1282 # ifdef DEAL_II_WITH_CXX14 1285 std::make_pair(cell->level(), cell->index()),
1290 std::make_pair(std::make_pair(cell->level(), cell->index()),
1295 cell, position_hint->second.get_location());
1296 position_hint->second.set_reference_location(p_unit);
1303 std::vector<
typename std::multimap<internal::LevelInd,
1306 for (
unsigned int child_index = 0;
1307 child_index < GeometryInfo<dim>::max_children_per_cell;
1311 cell->child(child_index);
1312 position_hints[child_index] =
particles.upper_bound(
1313 std::make_pair(child->level(), child->index()));
1316 for (
unsigned int i = 0; i < *n_particles_in_cell_ptr; ++i)
1320 for (
unsigned int child_index = 0;
1321 child_index < GeometryInfo<dim>::max_children_per_cell;
1325 child = cell->child(child_index);
1330 mapping->transform_real_to_unit_cell(child,
1334 p.set_reference_location(p_unit);
1339 # ifdef DEAL_II_WITH_CXX14 1340 position_hints[child_index] =
1341 particles.emplace_hint(position_hints[child_index],
1342 std::make_pair(child->level(),
1346 position_hints[child_index] =
particles.insert(
1347 position_hints[child_index],
1349 std::make_pair(child->level(), child->index()), p));
1351 ++position_hints[child_index];
1365 DEAL_II_NAMESPACE_CLOSE
1367 DEAL_II_NAMESPACE_OPEN
1369 #ifdef DEAL_II_WITH_P4EST 1370 # include "particle_handler.inst" 1373 DEAL_II_NAMESPACE_CLOSE
std::multimap< internal::LevelInd, Particle< dim, spacedim > > ghost_particles
void sort_particles_into_subdomains_and_cells()
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const CellStatus)> &pack_callback, const bool returns_variable_size_data)
static const unsigned int invalid_unsigned_int
types::particle_index n_global_max_particles_per_cell() const
void remove_particle(const particle_iterator &particle)
void initialize(const parallel::distributed::Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const unsigned int n_properties=0)
particle_iterator begin_ghost() const
std::function< const void *(const particle_iterator &, const void *)> load_callback
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
unsigned long long int particle_index
types::particle_index next_free_particle_index
SmartPointer< const parallel::distributed::Triangulation< dim, spacedim >, ParticleHandler< dim, spacedim > > triangulation
void update_cached_numbers()
const ArrayView< double > get_properties()
numbers::NumberTraits< Number >::real_type norm() const
#define AssertThrow(cond, exc)
unsigned int global_max_particles_per_cell
std::array< unsigned int, 4 > binary_type
void register_load_callback_function(const bool serialization)
static::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
void insert_particles(const std::multimap< typename Triangulation< dim, spacedim >::active_cell_iterator, Particle< dim, spacedim >> &particles)
particle_iterator end_ghost() const
types::particle_index global_number_of_particles
#define Assert(cond, exc)
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
std::function< std::size_t()> size_callback
Abstract base class for mapping classes.
particle_iterator begin() const
void register_store_callback_function()
void notify_ready_to_unpack(const unsigned int handle, const std::function< void(const cell_iterator &, const CellStatus, const boost::iterator_range< std::vector< char >::const_iterator > &)> &unpack_callback)
particle_iterator end() const
void exchange_ghost_particles()
PropertyPool & get_property_pool() const
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
std::unique_ptr< PropertyPool > property_pool
const types::subdomain_id artificial_subdomain_id
boost::iterator_range< particle_iterator > particle_iterator_range
ParticleIterator< dim, spacedim > particle_iterator
#define AssertThrowMPI(error_code)
std::multimap< internal::LevelInd, Particle< dim, spacedim > > particles
std::function< void *(const particle_iterator &, void *)> store_callback
unsigned int n_properties_per_particle() const
~ParticleHandler() override
types::particle_index n_locally_owned_particles() const
static::ExceptionBase & ExcNotImplemented()
types::particle_index n_global_particles() const
std::vector< char > store_particles(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status) const
SmartPointer< const Mapping< dim, spacedim >, ParticleHandler< dim, spacedim > > mapping
types::particle_index get_next_free_particle_index() const
unsigned int n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
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, 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 >>())
void load_particles(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status, const boost::iterator_range< std::vector< char >::const_iterator > &data_range)
bool has_properties() const
void register_additional_store_load_functions(const std::function< std::size_t()> &size_callback, const std::function< void *(const particle_iterator &, void *)> &store_callback, const std::function< const void *(const particle_iterator &, const void *)> &load_callback)
static::ExceptionBase & ExcInternalError()