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