Reference documentation for deal.II version 9.1.1
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
particle_handler.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 #include <deal.II/base/std_cxx14/memory.h>
17 
18 #include <deal.II/grid/grid_tools.h>
19 #include <deal.II/grid/grid_tools_cache.h>
20 
21 #include <deal.II/particles/particle_handler.h>
22 
23 #include <utility>
24 
25 DEAL_II_NAMESPACE_OPEN
26 
27 #ifdef DEAL_II_WITH_P4EST
28 
29 namespace Particles
30 {
31  template <int dim, int spacedim>
33  : triangulation()
34  , particles()
35  , ghost_particles()
36  , global_number_of_particles(0)
37  , global_max_particles_per_cell(0)
38  , next_free_particle_index(0)
39  , property_pool(new PropertyPool(0))
40  , size_callback()
41  , store_callback()
42  , load_callback()
43  , handle(numbers::invalid_unsigned_int)
44  {}
45 
46 
47 
48  template <int dim, int spacedim>
51  const Mapping<dim, spacedim> & mapping,
52  const unsigned int n_properties)
53  : triangulation(&triangulation, typeid(*this).name())
54  , mapping(&mapping, typeid(*this).name())
55  , particles()
56  , ghost_particles()
57  , global_number_of_particles(0)
58  , global_max_particles_per_cell(0)
59  , next_free_particle_index(0)
60  , property_pool(new PropertyPool(n_properties))
61  , size_callback()
62  , store_callback()
63  , load_callback()
64  , handle(numbers::invalid_unsigned_int)
65  {}
66 
67 
68 
69  template <int dim, int spacedim>
70  void
73  & new_triangulation,
74  const Mapping<dim, spacedim> &new_mapping,
75  const unsigned int n_properties)
76  {
77  triangulation = &new_triangulation;
78  mapping = &new_mapping;
79 
80  // Create the memory pool that will store all particle properties
81  property_pool = std_cxx14::make_unique<PropertyPool>(n_properties);
82  }
83 
84 
85 
86  template <int dim, int spacedim>
87  void
89  {
90  clear_particles();
91  global_number_of_particles = 0;
92  next_free_particle_index = 0;
93  global_max_particles_per_cell = 0;
94  }
95 
96 
97 
98  template <int dim, int spacedim>
99  void
101  {
102  particles.clear();
103  }
104 
105 
106 
107  template <int dim, int spacedim>
108  void
110  {
111  types::particle_index locally_highest_index = 0;
112  unsigned int local_max_particles_per_cell = 0;
113  unsigned int current_particles_per_cell = 0;
115  triangulation->begin_active();
116 
117  for (particle_iterator particle = begin(); particle != end(); ++particle)
118  {
119  locally_highest_index =
120  std::max(locally_highest_index, particle->get_id());
121 
122  if (particle->get_surrounding_cell(*triangulation) != current_cell)
123  {
124  current_particles_per_cell = 0;
125  current_cell = particle->get_surrounding_cell(*triangulation);
126  }
127 
128  ++current_particles_per_cell;
129  local_max_particles_per_cell =
130  std::max(local_max_particles_per_cell, current_particles_per_cell);
131  }
132 
133  global_number_of_particles =
134  ::Utilities::MPI::sum(particles.size(),
135  triangulation->get_communicator());
136  next_free_particle_index =
137  ::Utilities::MPI::max(locally_highest_index,
138  triangulation->get_communicator()) +
139  1;
140  global_max_particles_per_cell =
141  ::Utilities::MPI::max(local_max_particles_per_cell,
142  triangulation->get_communicator());
143  }
144 
145 
146 
147  template <int dim, int spacedim>
150  {
151  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin();
152  }
153 
154 
155 
156  template <int dim, int spacedim>
159  {
160  return particle_iterator(particles, particles.begin());
161  }
162 
163 
164 
165  template <int dim, int spacedim>
168  {
169  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end();
170  }
171 
172 
173 
174  template <int dim, int spacedim>
177  {
178  return particle_iterator(particles, particles.end());
179  }
180 
181 
182 
183  template <int dim, int spacedim>
186  {
187  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin_ghost();
188  }
189 
190 
191 
192  template <int dim, int spacedim>
195  {
196  return particle_iterator(ghost_particles, ghost_particles.begin());
197  }
198 
199 
200 
201  template <int dim, int spacedim>
204  {
205  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end_ghost();
206  }
207 
208 
209 
210  template <int dim, int spacedim>
213  {
214  return particle_iterator(ghost_particles, ghost_particles.end());
215  }
216 
217 
218 
219  template <int dim, int spacedim>
223  const
224  {
225  return (const_cast<ParticleHandler<dim, spacedim> *>(this))
226  ->particles_in_cell(cell);
227  }
228 
229 
230 
231  template <int dim, int spacedim>
235  {
236  const internal::LevelInd level_index =
237  std::make_pair<int, int>(cell->level(), cell->index());
238 
239  if (cell->is_ghost())
240  {
241  const auto particles_in_cell = ghost_particles.equal_range(level_index);
242  return boost::make_iterator_range(
243  particle_iterator(ghost_particles, particles_in_cell.first),
244  particle_iterator(ghost_particles, particles_in_cell.second));
245  }
246 
247  const auto particles_in_cell = particles.equal_range(level_index);
248  return boost::make_iterator_range(
249  particle_iterator(particles, particles_in_cell.first),
250  particle_iterator(particles, particles_in_cell.second));
251  }
252 
253 
254 
255  template <int dim, int spacedim>
256  void
259  {
260  particles.erase(particle->particle);
261  }
262 
263 
264 
265  template <int dim, int spacedim>
268  const Particle<dim, spacedim> & particle,
270  {
271  typename std::multimap<internal::LevelInd,
272  Particle<dim, spacedim>>::iterator it =
273  particles.insert(
274  std::make_pair(internal::LevelInd(cell->level(), cell->index()),
275  particle));
276 
277  particle_iterator particle_it(particles, it);
278  particle_it->set_property_pool(*property_pool);
279 
280  if (particle.has_properties())
281  for (unsigned int n = 0; n < particle.get_properties().size(); ++n)
282  particle_it->get_properties()[n] = particle.get_properties()[n];
283 
284  return particle_it;
285  }
286 
287 
288 
289  template <int dim, int spacedim>
290  void
292  const std::multimap<
294  Particle<dim, spacedim>> &new_particles)
295  {
296  for (auto particle = new_particles.begin(); particle != new_particles.end();
297  ++particle)
298  particles.insert(
299  particles.end(),
300  std::make_pair(internal::LevelInd(particle->first->level(),
301  particle->first->index()),
302  particle->second));
303 
304  update_cached_numbers();
305  }
306 
307 
308 
309  template <int dim, int spacedim>
310  void
312  const std::vector<Point<spacedim>> &positions)
313  {
314  update_cached_numbers();
315 
316  // Determine the starting particle index of this process, which
317  // is the highest currently existing particle index plus the sum
318  // of the number of newly generated particles of all
319  // processes with a lower rank if in a parallel computation.
320  const types::particle_index local_next_particle_index =
321  get_next_free_particle_index();
322  types::particle_index local_start_index = 0;
323 
324 # ifdef DEAL_II_WITH_MPI
325  types::particle_index particles_to_add_locally = positions.size();
326  const int ierr = MPI_Scan(&particles_to_add_locally,
327  &local_start_index,
328  1,
329  DEAL_II_PARTICLE_INDEX_MPI_TYPE,
330  MPI_SUM,
331  triangulation->get_communicator());
332  AssertThrowMPI(ierr);
333  local_start_index -= particles_to_add_locally;
334 # endif
335 
336  local_start_index += local_next_particle_index;
337 
338  GridTools::Cache<dim, spacedim> cache(*triangulation, *mapping);
339  auto point_locations = GridTools::compute_point_locations(cache, positions);
340 
341  auto &cells = std::get<0>(point_locations);
342  auto &local_positions = std::get<1>(point_locations);
343  auto &index_map = std::get<2>(point_locations);
344 
345  if (cells.size() == 0)
346  return;
347 
348  auto hint =
349  particles.find(std::make_pair(cells[0]->level(), cells[0]->index()));
350  for (unsigned int i = 0; i < cells.size(); ++i)
351  {
352  internal::LevelInd current_cell(cells[i]->level(), cells[i]->index());
353  for (unsigned int p = 0; p < local_positions[i].size(); ++p)
354  {
355  hint = particles.insert(
356  hint,
357  std::make_pair(current_cell,
358  Particle<dim, spacedim>(positions[index_map[i][p]],
359  local_positions[i][p],
360  local_start_index +
361  index_map[i][p])));
362  }
363  }
364 
365  update_cached_numbers();
366  }
367 
368 
369 
370  template <int dim, int spacedim>
373  {
374  return global_number_of_particles;
375  }
376 
377 
378 
379  template <int dim, int spacedim>
382  {
383  return global_max_particles_per_cell;
384  }
385 
386 
387 
388  template <int dim, int spacedim>
391  {
392  return particles.size();
393  }
394 
395 
396 
397  template <int dim, int spacedim>
398  unsigned int
400  {
401  return property_pool->n_properties_per_slot();
402  }
403 
404 
405 
406  template <int dim, int spacedim>
407  unsigned int
410  const
411  {
412  const internal::LevelInd found_cell =
413  std::make_pair<int, int>(cell->level(), cell->index());
414 
415  if (cell->is_locally_owned())
416  return particles.count(found_cell);
417  else if (cell->is_ghost())
418  return ghost_particles.count(found_cell);
419  else if (cell->is_artificial())
420  AssertThrow(false, ExcInternalError());
421 
422  return 0;
423  }
424 
425 
426 
427  template <int dim, int spacedim>
430  {
431  return next_free_particle_index;
432  }
433 
434 
435 
436  template <int dim, int spacedim>
437  PropertyPool &
439  {
440  return *property_pool;
441  }
442 
443 
444 
445  namespace
446  {
456  template <int dim>
457  bool
458  compare_particle_association(
459  const unsigned int a,
460  const unsigned int b,
461  const Tensor<1, dim> & particle_direction,
462  const std::vector<Tensor<1, dim>> &center_directions)
463  {
464  const double scalar_product_a = center_directions[a] * particle_direction;
465  const double scalar_product_b = center_directions[b] * particle_direction;
466 
467  // The function is supposed to return if a is before b. We are looking
468  // for the alignment of particle direction and center direction,
469  // therefore return if the scalar product of a is larger.
470  return (scalar_product_a > scalar_product_b);
471  }
472  } // namespace
473 
474 
475 
476  template <int dim, int spacedim>
477  void
479  {
480  // TODO: The current algorithm only works for particles that are in
481  // the local domain or in ghost cells, because it only knows the
482  // subdomain_id of ghost cells, but not of artificial cells. This
483  // can be extended using the distributed version of compute point
484  // locations.
485  // TODO: Extend this function to allow keeping particles on other
486  // processes around (with an invalid cell).
487 
488  std::vector<particle_iterator> particles_out_of_cell;
489  particles_out_of_cell.reserve(n_locally_owned_particles());
490 
491  // Now update the reference locations of the moved particles
492  for (particle_iterator it = begin(); it != end(); ++it)
493  {
494  const typename Triangulation<dim, spacedim>::cell_iterator cell =
495  it->get_surrounding_cell(*triangulation);
496 
497  try
498  {
499  const Point<dim> p_unit =
500  mapping->transform_real_to_unit_cell(cell, it->get_location());
502  {
503  it->set_reference_location(p_unit);
504  }
505  else
506  {
507  // The particle has left the cell
508  particles_out_of_cell.push_back(it);
509  }
510  }
511  catch (typename Mapping<dim>::ExcTransformationFailed &)
512  {
513  // The particle has left the cell
514  particles_out_of_cell.push_back(it);
515  }
516  }
517 
518  // There are three reasons why a particle is not in its old cell:
519  // It moved to another cell, to another subdomain or it left the mesh.
520  // Particles that moved to another cell are updated and stored inside the
521  // sorted_particles vector, particles that moved to another domain are
522  // collected in the moved_particles_domain vector. Particles that left
523  // the mesh completely are ignored and removed.
524  std::vector<std::pair<internal::LevelInd, Particle<dim, spacedim>>>
525  sorted_particles;
526  std::map<types::subdomain_id, std::vector<particle_iterator>>
527  moved_particles;
528  std::map<
530  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
531  moved_cells;
532 
533  // We do not know exactly how many particles are lost, exchanged between
534  // domains, or remain on this process. Therefore we pre-allocate approximate
535  // sizes for these vectors. If more space is needed an automatic and
536  // relatively fast (compared to other parts of this algorithm)
537  // re-allocation will happen.
538  using vector_size = typename std::vector<particle_iterator>::size_type;
539  sorted_particles.reserve(
540  static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
541 
542  const std::set<types::subdomain_id> ghost_owners =
543  triangulation->ghost_owners();
544 
545  for (const auto ghost_owner : ghost_owners)
546  moved_particles[ghost_owner].reserve(
547  static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
548  for (const auto ghost_owner : ghost_owners)
549  moved_cells[ghost_owner].reserve(
550  static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
551 
552  {
553  // Create a map from vertices to adjacent cells
554  const std::vector<
555  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
556  vertex_to_cells(GridTools::vertex_to_cell_map(*triangulation));
557 
558  // Create a corresponding map of vectors from vertex to cell center
559  const std::vector<std::vector<Tensor<1, spacedim>>>
560  vertex_to_cell_centers(
562  vertex_to_cells));
563 
564  std::vector<unsigned int> neighbor_permutation;
565 
566  // Find the cells that the particles moved to.
567  typename std::vector<particle_iterator>::iterator
568  it = particles_out_of_cell.begin(),
569  end_particle = particles_out_of_cell.end();
570 
571  for (; it != end_particle; ++it)
572  {
573  // The cell the particle is in
574  Point<dim> current_reference_position;
575  bool found_cell = false;
576 
577  // Check if the particle is in one of the old cell's neighbors
578  // that are adjacent to the closest vertex
580  current_cell = (*it)->get_surrounding_cell(*triangulation);
581 
582  const unsigned int closest_vertex =
583  GridTools::find_closest_vertex_of_cell<dim, spacedim>(
584  current_cell, (*it)->get_location());
585  Tensor<1, spacedim> vertex_to_particle =
586  (*it)->get_location() - current_cell->vertex(closest_vertex);
587  vertex_to_particle /= vertex_to_particle.norm();
588 
589  const unsigned int closest_vertex_index =
590  current_cell->vertex_index(closest_vertex);
591  const unsigned int n_neighbor_cells =
592  vertex_to_cells[closest_vertex_index].size();
593 
594  neighbor_permutation.resize(n_neighbor_cells);
595  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
596  neighbor_permutation[i] = i;
597 
598  std::sort(neighbor_permutation.begin(),
599  neighbor_permutation.end(),
600  std::bind(&compare_particle_association<spacedim>,
601  std::placeholders::_1,
602  std::placeholders::_2,
603  std::cref(vertex_to_particle),
604  std::cref(
605  vertex_to_cell_centers[closest_vertex_index])));
606 
607  // Search all of the cells adjacent to the closest vertex of the
608  // previous cell Most likely we will find the particle in them.
609  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
610  {
611  try
612  {
613  typename std::set<typename Triangulation<dim, spacedim>::
614  active_cell_iterator>::const_iterator
615  cell = vertex_to_cells[closest_vertex_index].begin();
616 
617  std::advance(cell, neighbor_permutation[i]);
618  const Point<dim> p_unit =
619  mapping->transform_real_to_unit_cell(*cell,
620  (*it)->get_location());
622  {
623  current_cell = *cell;
624  current_reference_position = p_unit;
625  found_cell = true;
626  break;
627  }
628  }
629  catch (typename Mapping<dim>::ExcTransformationFailed &)
630  {}
631  }
632 
633  if (!found_cell)
634  {
635  // The particle is not in a neighbor of the old cell.
636  // Look for the new cell in the whole local domain.
637  // This case is rare.
638  try
639  {
640  const std::pair<const typename Triangulation<dim, spacedim>::
641  active_cell_iterator,
642  Point<dim>>
643  current_cell_and_position =
644  GridTools::find_active_cell_around_point<>(
645  *mapping, *triangulation, (*it)->get_location());
646  current_cell = current_cell_and_position.first;
647  current_reference_position = current_cell_and_position.second;
648  }
649  catch (GridTools::ExcPointNotFound<spacedim> &)
650  {
651  // We can find no cell for this particle. It has left the
652  // domain due to an integration error or an open boundary.
653  continue;
654  }
655  }
656 
657  // If we are here, we found a cell and reference position for this
658  // particle
659  (*it)->set_reference_location(current_reference_position);
660 
661  // Reinsert the particle into our domain if we own its cell.
662  // Mark it for MPI transfer otherwise
663  if (current_cell->is_locally_owned())
664  {
665  sorted_particles.push_back(
666  std::make_pair(internal::LevelInd(current_cell->level(),
667  current_cell->index()),
668  (*it)->particle->second));
669  }
670  else
671  {
672  moved_particles[current_cell->subdomain_id()].push_back(*it);
673  moved_cells[current_cell->subdomain_id()].push_back(current_cell);
674  }
675  }
676  }
677 
678  // Sort the updated particles. This pre-sort speeds up inserting
679  // them into particles to O(N) complexity.
680  std::multimap<internal::LevelInd, Particle<dim, spacedim>>
681  sorted_particles_map;
682 
683  // Exchange particles between processors if we have more than one process
684 # ifdef DEAL_II_WITH_MPI
686  triangulation->get_communicator()) > 1)
687  send_recv_particles(moved_particles, sorted_particles_map, moved_cells);
688 # endif
689 
690  sorted_particles_map.insert(sorted_particles.begin(),
691  sorted_particles.end());
692 
693  for (unsigned int i = 0; i < particles_out_of_cell.size(); ++i)
694  remove_particle(particles_out_of_cell[i]);
695 
696  particles.insert(sorted_particles_map.begin(), sorted_particles_map.end());
697  update_cached_numbers();
698  }
699 
700 
701 
702  template <int dim, int spacedim>
703  void
705  {
706  // Nothing to do in serial computations
708  triangulation->get_communicator()) == 1)
709  return;
710 
711 # ifdef DEAL_II_WITH_MPI
712  // First clear the current ghost_particle information
713  ghost_particles.clear();
714 
715  std::map<types::subdomain_id, std::vector<particle_iterator>>
716  ghost_particles_by_domain;
717 
718  const std::set<types::subdomain_id> ghost_owners =
719  triangulation->ghost_owners();
720  for (const auto ghost_owner : ghost_owners)
721  ghost_particles_by_domain[ghost_owner].reserve(
722  static_cast<typename std::vector<particle_iterator>::size_type>(
723  particles.size() * 0.25));
724 
725  std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain(
726  triangulation->n_vertices());
727 
728  for (const auto &cell : triangulation->active_cell_iterators())
729  {
730  if (cell->is_ghost())
731  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
732  ++v)
733  vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(
734  cell->subdomain_id());
735  }
736 
737  for (const auto &cell : triangulation->active_cell_iterators())
738  {
739  if (!cell->is_ghost())
740  {
741  std::set<unsigned int> cell_to_neighbor_subdomain;
742  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
743  ++v)
744  {
745  cell_to_neighbor_subdomain.insert(
746  vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
747  vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
748  }
749 
750  if (cell_to_neighbor_subdomain.size() > 0)
751  {
752  const particle_iterator_range particle_range =
753  particles_in_cell(cell);
754 
755  for (const auto domain : cell_to_neighbor_subdomain)
756  {
757  for (typename particle_iterator_range::iterator particle =
758  particle_range.begin();
759  particle != particle_range.end();
760  ++particle)
761  ghost_particles_by_domain[domain].push_back(particle);
762  }
763  }
764  }
765  }
766 
767  send_recv_particles(ghost_particles_by_domain, ghost_particles);
768 # endif
769  }
770 
771 
772 
773 # ifdef DEAL_II_WITH_MPI
774  template <int dim, int spacedim>
775  void
777  const std::map<types::subdomain_id, std::vector<particle_iterator>>
778  &particles_to_send,
779  std::multimap<internal::LevelInd, Particle<dim, spacedim>>
780  &received_particles,
781  const std::map<
784  &send_cells)
785  {
786  // Determine the communication pattern
787  const std::set<types::subdomain_id> ghost_owners =
788  triangulation->ghost_owners();
789  const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
790  ghost_owners.end());
791  const unsigned int n_neighbors = neighbors.size();
792 
793  if (send_cells.size() != 0)
794  Assert(particles_to_send.size() == send_cells.size(), ExcInternalError());
795 
796  // If we do not know the subdomain this particle needs to be send to, throw
797  // an error
798  Assert(particles_to_send.find(numbers::artificial_subdomain_id) ==
799  particles_to_send.end(),
800  ExcInternalError());
801 
802  // TODO: Implement the shipping of particles to processes that are not ghost
803  // owners of the local domain
804  for (auto send_particles = particles_to_send.begin();
805  send_particles != particles_to_send.end();
806  ++send_particles)
807  Assert(ghost_owners.find(send_particles->first) != ghost_owners.end(),
809 
810  unsigned int n_send_particles = 0;
811  for (auto send_particles = particles_to_send.begin();
812  send_particles != particles_to_send.end();
813  ++send_particles)
814  n_send_particles += send_particles->second.size();
815 
816  const unsigned int cellid_size = sizeof(CellId::binary_type);
817 
818  // Containers for the amount and offsets of data we will send
819  // to other processors and the data itself.
820  std::vector<unsigned int> n_send_data(n_neighbors, 0);
821  std::vector<unsigned int> send_offsets(n_neighbors, 0);
822  std::vector<char> send_data;
823 
824  // Only serialize things if there are particles to be send.
825  // We can not return early even if no particles
826  // are send, because we might receive particles from other processes
827  if (n_send_particles > 0)
828  {
829  // Allocate space for sending particle data
830  const unsigned int particle_size =
831  begin()->serialized_size_in_bytes() + cellid_size +
832  (size_callback ? size_callback() : 0);
833  send_data.resize(n_send_particles * particle_size);
834  void *data = static_cast<void *>(&send_data.front());
835 
836  // Serialize the data sorted by receiving process
837  for (unsigned int i = 0; i < n_neighbors; ++i)
838  {
839  send_offsets[i] = reinterpret_cast<std::size_t>(data) -
840  reinterpret_cast<std::size_t>(&send_data.front());
841 
842  for (unsigned int j = 0;
843  j < particles_to_send.at(neighbors[i]).size();
844  ++j)
845  {
846  // If no target cells are given, use the iterator information
848  cell;
849  if (send_cells.size() == 0)
850  cell =
851  particles_to_send.at(neighbors[i])[j]->get_surrounding_cell(
852  *triangulation);
853  else
854  cell = send_cells.at(neighbors[i])[j];
855 
856  const CellId::binary_type cellid =
857  cell->id().template to_binary<dim>();
858  memcpy(data, &cellid, cellid_size);
859  data = static_cast<char *>(data) + cellid_size;
860 
861  particles_to_send.at(neighbors[i])[j]->write_data(data);
862  if (store_callback)
863  data =
864  store_callback(particles_to_send.at(neighbors[i])[j], data);
865  }
866  n_send_data[i] = reinterpret_cast<std::size_t>(data) -
867  send_offsets[i] -
868  reinterpret_cast<std::size_t>(&send_data.front());
869  }
870  }
871 
872  // Containers for the data we will receive from other processors
873  std::vector<unsigned int> n_recv_data(n_neighbors);
874  std::vector<unsigned int> recv_offsets(n_neighbors);
875 
876  // Notify other processors how many particles we will send
877  {
878  std::vector<MPI_Request> n_requests(2 * n_neighbors);
879  for (unsigned int i = 0; i < n_neighbors; ++i)
880  {
881  const int ierr = MPI_Irecv(&(n_recv_data[i]),
882  1,
883  MPI_UNSIGNED,
884  neighbors[i],
885  0,
886  triangulation->get_communicator(),
887  &(n_requests[2 * i]));
888  AssertThrowMPI(ierr);
889  }
890  for (unsigned int i = 0; i < n_neighbors; ++i)
891  {
892  const int ierr = MPI_Isend(&(n_send_data[i]),
893  1,
894  MPI_UNSIGNED,
895  neighbors[i],
896  0,
897  triangulation->get_communicator(),
898  &(n_requests[2 * i + 1]));
899  AssertThrowMPI(ierr);
900  }
901  const int ierr =
902  MPI_Waitall(2 * n_neighbors, n_requests.data(), MPI_STATUSES_IGNORE);
903  AssertThrowMPI(ierr);
904  }
905 
906  // Determine how many particles and data we will receive
907  unsigned int total_recv_data = 0;
908  for (unsigned int neighbor_id = 0; neighbor_id < n_neighbors; ++neighbor_id)
909  {
910  recv_offsets[neighbor_id] = total_recv_data;
911  total_recv_data += n_recv_data[neighbor_id];
912  }
913 
914  // Set up the space for the received particle data
915  std::vector<char> recv_data(total_recv_data);
916 
917  // Exchange the particle data between domains
918  {
919  std::vector<MPI_Request> requests(2 * n_neighbors);
920  unsigned int send_ops = 0;
921  unsigned int recv_ops = 0;
922 
923  for (unsigned int i = 0; i < n_neighbors; ++i)
924  if (n_recv_data[i] > 0)
925  {
926  const int ierr = MPI_Irecv(&(recv_data[recv_offsets[i]]),
927  n_recv_data[i],
928  MPI_CHAR,
929  neighbors[i],
930  1,
931  triangulation->get_communicator(),
932  &(requests[send_ops]));
933  AssertThrowMPI(ierr);
934  send_ops++;
935  }
936 
937  for (unsigned int i = 0; i < n_neighbors; ++i)
938  if (n_send_data[i] > 0)
939  {
940  const int ierr = MPI_Isend(&(send_data[send_offsets[i]]),
941  n_send_data[i],
942  MPI_CHAR,
943  neighbors[i],
944  1,
945  triangulation->get_communicator(),
946  &(requests[send_ops + recv_ops]));
947  AssertThrowMPI(ierr);
948  recv_ops++;
949  }
950  const int ierr =
951  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
952  AssertThrowMPI(ierr);
953  }
954 
955  // Put the received particles into the domain if they are in the
956  // triangulation
957  const void *recv_data_it = static_cast<const void *>(recv_data.data());
958 
959  while (reinterpret_cast<std::size_t>(recv_data_it) -
960  reinterpret_cast<std::size_t>(recv_data.data()) <
961  total_recv_data)
962  {
963  CellId::binary_type binary_cellid;
964  memcpy(&binary_cellid, recv_data_it, cellid_size);
965  const CellId id(binary_cellid);
966  recv_data_it = static_cast<const char *>(recv_data_it) + cellid_size;
967 
969  id.to_cell(*triangulation);
970 
971  typename std::multimap<internal::LevelInd,
972  Particle<dim, spacedim>>::iterator
973  recv_particle = received_particles.insert(std::make_pair(
974  internal::LevelInd(cell->level(), cell->index()),
975  Particle<dim, spacedim>(recv_data_it, property_pool.get())));
976 
977  if (load_callback)
978  recv_data_it =
979  load_callback(particle_iterator(received_particles, recv_particle),
980  recv_data_it);
981  }
982 
983  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
984  ExcMessage(
985  "The amount of data that was read into new particles "
986  "does not match the amount of data sent around."));
987  }
988 # endif
989 
990 
991 
992  template <int dim, int spacedim>
993  void
995  const std::function<std::size_t()> & size_callb,
996  const std::function<void *(const particle_iterator &, void *)> &store_callb,
997  const std::function<const void *(const particle_iterator &, const void *)>
998  &load_callb)
999  {
1000  size_callback = size_callb;
1001  store_callback = store_callb;
1002  load_callback = load_callb;
1003  }
1004 
1005 
1006 
1007  template <int dim, int spacedim>
1008  void
1010  {
1012  *non_const_triangulation =
1014  &(*triangulation));
1015 
1016  // Only save and load particles if there are any, we might get here for
1017  // example if somebody created a ParticleHandler but generated 0 particles.
1018  update_cached_numbers();
1019 
1020  if (global_max_particles_per_cell > 0)
1021  {
1022  const std::function<std::vector<char>(
1025  callback_function =
1027  std::cref(*this),
1028  std::placeholders::_1,
1029  std::placeholders::_2);
1030 
1031  handle = non_const_triangulation->register_data_attach(
1032  callback_function, /*returns_variable_size_data=*/true);
1033  }
1034  }
1035 
1036 
1037 
1038  template <int dim, int spacedim>
1039  void
1041  const bool serialization)
1042  {
1043  // All particles have been stored, when we reach this point. Empty the
1044  // particle data.
1045  clear_particles();
1046 
1048  *non_const_triangulation =
1050  &(*triangulation));
1051 
1052  // If we are resuming from a checkpoint, we first have to register the
1053  // store function again, to set the triangulation in the same state as
1054  // before the serialization. Only by this it knows how to deserialize the
1055  // data correctly. Only do this if something was actually stored.
1056  if (serialization && (global_max_particles_per_cell > 0))
1057  {
1058  const std::function<std::vector<char>(
1061  callback_function =
1063  std::cref(*this),
1064  std::placeholders::_1,
1065  std::placeholders::_2);
1066 
1067  handle = non_const_triangulation->register_data_attach(
1068  callback_function, /*returns_variable_size_data=*/true);
1069  }
1070 
1071  // Check if something was stored and load it
1072  if (handle != numbers::invalid_unsigned_int)
1073  {
1074  const std::function<void(
1077  const boost::iterator_range<std::vector<char>::const_iterator> &)>
1078  callback_function =
1080  std::ref(*this),
1081  std::placeholders::_1,
1082  std::placeholders::_2,
1083  std::placeholders::_3);
1084 
1085  non_const_triangulation->notify_ready_to_unpack(handle,
1086  callback_function);
1087 
1088  // Reset handle and update global number of particles. The number
1089  // can change because of discarded or newly generated particles
1091  update_cached_numbers();
1092  }
1093  }
1094 
1095 
1096 
1097  template <int dim, int spacedim>
1098  std::vector<char>
1100  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1101  const typename Triangulation<dim, spacedim>::CellStatus status) const
1102  {
1103  std::vector<Particle<dim, spacedim>> stored_particles_on_cell;
1104 
1105  switch (status)
1106  {
1109  // If the cell persist or is refined store all particles of the
1110  // current cell.
1111  {
1112  unsigned int n_particles = 0;
1113 
1114  const internal::LevelInd level_index = {cell->level(),
1115  cell->index()};
1116  const auto particles_in_cell =
1117  (cell->is_ghost() ? ghost_particles.equal_range(level_index) :
1118  particles.equal_range(level_index));
1119 
1120  n_particles = n_particles_in_cell(cell);
1121  stored_particles_on_cell.reserve(n_particles);
1122 
1123  std::for_each(
1124  particles_in_cell.first,
1125  particles_in_cell.second,
1126  [&stored_particles_on_cell](
1127  const std::pair<internal::LevelInd, Particle<dim, spacedim>>
1128  &particle) {
1129  stored_particles_on_cell.push_back(particle.second);
1130  });
1131 
1132  AssertDimension(n_particles, stored_particles_on_cell.size());
1133  }
1134  break;
1135 
1137  // If this cell is the parent of children that will be coarsened,
1138  // collect the particles of all children.
1139  {
1140  unsigned int n_particles = 0;
1141 
1142  for (unsigned int child_index = 0;
1143  child_index < GeometryInfo<dim>::max_children_per_cell;
1144  ++child_index)
1145  {
1147  child = cell->child(child_index);
1148  n_particles += n_particles_in_cell(child);
1149  }
1150 
1151  stored_particles_on_cell.reserve(n_particles);
1152 
1153  for (unsigned int child_index = 0;
1154  child_index < GeometryInfo<dim>::max_children_per_cell;
1155  ++child_index)
1156  {
1158  child = cell->child(child_index);
1159  const internal::LevelInd level_index = {child->level(),
1160  child->index()};
1161  const auto particles_in_cell =
1162  (child->is_ghost() ?
1163  ghost_particles.equal_range(level_index) :
1164  particles.equal_range(level_index));
1165 
1166  std::for_each(
1167  particles_in_cell.first,
1168  particles_in_cell.second,
1169  [&stored_particles_on_cell](
1170  const std::pair<internal::LevelInd, Particle<dim, spacedim>>
1171  &particle) {
1172  stored_particles_on_cell.push_back(particle.second);
1173  });
1174  }
1175 
1176  AssertDimension(n_particles, stored_particles_on_cell.size());
1177  }
1178  break;
1179 
1180  default:
1181  Assert(false, ExcInternalError());
1182  break;
1183  }
1184 
1185  return Utilities::pack(stored_particles_on_cell,
1186  /*allow_compression=*/true);
1187  }
1188 
1189  template <int dim, int spacedim>
1190  void
1192  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
1193  const typename Triangulation<dim, spacedim>::CellStatus status,
1194  const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
1195  {
1196  // We leave this container non-const to be able to `std::move`
1197  // its contents directly into the particles multimap later.
1198  std::vector<Particle<dim, spacedim>> loaded_particles_on_cell =
1199  Utilities::unpack<std::vector<Particle<dim, spacedim>>>(
1200  data_range.begin(),
1201  data_range.end(),
1202  /*allow_compression=*/true);
1203 
1204  // Update the reference to the current property pool for all particles.
1205  // This was not stored, because they might be transported across process
1206  // domains.
1207  for (auto &particle : loaded_particles_on_cell)
1208  particle.set_property_pool(*property_pool);
1209 
1210  switch (status)
1211  {
1213  {
1214  auto position_hint = particles.end();
1215  for (const auto &particle : loaded_particles_on_cell)
1216  {
1217  // Use std::multimap::emplace_hint to speed up insertion of
1218  // particles. This is a C++11 function, but not all compilers
1219  // that report a -std=c++11 (like gcc 4.6) implement it, so
1220  // require C++14 instead.
1221 # ifdef DEAL_II_WITH_CXX14
1222  position_hint =
1223  particles.emplace_hint(position_hint,
1224  std::make_pair(cell->level(),
1225  cell->index()),
1226  std::move(particle));
1227 # else
1228  position_hint =
1229  particles.insert(position_hint,
1230  std::make_pair(std::make_pair(cell->level(),
1231  cell->index()),
1232  std::move(particle)));
1233 # endif
1234  // Move the hint position forward by one, i.e., for the next
1235  // particle. The 'hint' position will thus be right after the
1236  // one just inserted.
1237  ++position_hint;
1238  }
1239  }
1240  break;
1241 
1243  {
1244  typename std::multimap<internal::LevelInd,
1245  Particle<dim, spacedim>>::iterator
1246  position_hint = particles.end();
1247  for (auto &particle : loaded_particles_on_cell)
1248  {
1249  const Point<dim> p_unit =
1250  mapping->transform_real_to_unit_cell(cell,
1251  particle.get_location());
1252  particle.set_reference_location(p_unit);
1253  // Use std::multimap::emplace_hint to speed up insertion of
1254  // particles. This is a C++11 function, but not all compilers
1255  // that report a -std=c++11 (like gcc 4.6) implement it, so
1256  // require C++14 instead.
1257 # ifdef DEAL_II_WITH_CXX14
1258  position_hint =
1259  particles.emplace_hint(position_hint,
1260  std::make_pair(cell->level(),
1261  cell->index()),
1262  std::move(particle));
1263 # else
1264  position_hint =
1265  particles.insert(position_hint,
1266  std::make_pair(std::make_pair(cell->level(),
1267  cell->index()),
1268  std::move(particle)));
1269 # endif
1270  // Move the hint position forward by one, i.e., for the next
1271  // particle. The 'hint' position will thus be right after the
1272  // one just inserted.
1273  ++position_hint;
1274  }
1275  }
1276  break;
1277 
1279  {
1280  std::vector<
1281  typename std::multimap<internal::LevelInd,
1282  Particle<dim, spacedim>>::iterator>
1284  for (unsigned int child_index = 0;
1285  child_index < GeometryInfo<dim>::max_children_per_cell;
1286  ++child_index)
1287  {
1289  child = cell->child(child_index);
1290  position_hints[child_index] = particles.upper_bound(
1291  std::make_pair(child->level(), child->index()));
1292  }
1293 
1294  for (auto &particle : loaded_particles_on_cell)
1295  {
1296  for (unsigned int child_index = 0;
1297  child_index < GeometryInfo<dim>::max_children_per_cell;
1298  ++child_index)
1299  {
1301  child = cell->child(child_index);
1302 
1303  try
1304  {
1305  const Point<dim> p_unit =
1306  mapping->transform_real_to_unit_cell(
1307  child, particle.get_location());
1309  {
1310  particle.set_reference_location(p_unit);
1311  // Use std::multimap::emplace_hint to speed up
1312  // insertion of particles. This is a C++11 function,
1313  // but not all compilers that report a -std=c++11
1314  // (like gcc 4.6) implement it, so require C++14
1315  // instead.
1316 # ifdef DEAL_II_WITH_CXX14
1317  position_hints[child_index] =
1318  particles.emplace_hint(
1319  position_hints[child_index],
1320  std::make_pair(child->level(), child->index()),
1321  std::move(particle));
1322 # else
1323  position_hints[child_index] = particles.insert(
1324  position_hints[child_index],
1325  std::make_pair(std::make_pair(child->level(),
1326  child->index()),
1327  std::move(particle)));
1328 # endif
1329  // Move the hint position forward by one, i.e., for
1330  // the next particle. The 'hint' position will thus
1331  // be right after the one just inserted.
1332  ++position_hints[child_index];
1333  break;
1334  }
1335  }
1336  catch (typename Mapping<dim>::ExcTransformationFailed &)
1337  {}
1338  }
1339  }
1340  }
1341  break;
1342 
1343  default:
1344  Assert(false, ExcInternalError());
1345  break;
1346  }
1347  }
1348 } // namespace Particles
1349 
1350 #endif
1351 
1352 DEAL_II_NAMESPACE_CLOSE
1353 
1354 DEAL_II_NAMESPACE_OPEN
1355 
1356 #ifdef DEAL_II_WITH_P4EST
1357 # include "particle_handler.inst"
1358 #endif
1359 
1360 DEAL_II_NAMESPACE_CLOSE
types::particle_index n_global_particles() const
unsigned int register_data_attach(const std::function< std::vector< char >(const cell_iterator &, const CellStatus)> &pack_callback, const bool returns_variable_size_data)
Definition: tria.cc:4357
static const unsigned int invalid_unsigned_int
Definition: types.h:173
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1567
types::particle_index n_locally_owned_particles() const
std::vector< char > store_particles(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status) 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)
types::particle_index n_global_max_particles_per_cell() const
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > > vertex_to_cell_map(const Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2152
const ArrayView< double > get_properties()
Definition: particle.cc:327
#define AssertThrow(cond, exc)
Definition: exceptions.h:1519
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1318
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:73
cell_iterator end() const
Definition: tria.cc:11949
unsigned int particle_index
Definition: particle.h:64
void register_load_callback_function(const bool serialization)
static ::ExceptionBase & ExcMessage(std::string arg1)
unsigned int subdomain_id
Definition: types.h:43
void insert_particles(const std::multimap< typename Triangulation< dim, spacedim >::active_cell_iterator, Particle< dim, spacedim >> &particles)
particle_iterator begin_ghost() const
unsigned int n_properties_per_particle() const
#define Assert(cond, exc)
Definition: exceptions.h:1407
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
Abstract base class for mapping classes.
Definition: dof_tools.h:57
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)
Definition: tria.cc:4387
particle_iterator begin() const
size_t pack(const T &object, std::vector< char > &dest_buffer, const bool allow_compression=true)
Definition: utilities.h:1174
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:71
Definition: cell_id.h:63
std::vector< std::vector< Tensor< 1, spacedim > > > vertex_to_cell_centers_directions(const Triangulation< dim, spacedim > &mesh, const std::vector< std::set< typename Triangulation< dim, spacedim >::active_cell_iterator >> &vertex_to_cells)
Definition: grid_tools.cc:1597
const types::subdomain_id artificial_subdomain_id
Definition: types.h:275
types::particle_index get_next_free_particle_index() const
boost::iterator_range< particle_iterator > particle_iterator_range
bool has_properties() const
Definition: particle.cc:349
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1695
PropertyPool & get_property_pool() const
static ::ExceptionBase & ExcNotImplemented()
return_type compute_point_locations(const Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &points, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell_hint=typename Triangulation< dim, spacedim >::active_cell_iterator())
Definition: grid_tools.cc:4328
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)
particle_iterator end_ghost() const
particle_iterator end() 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)
unsigned int n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
static ::ExceptionBase & ExcInternalError()