Reference documentation for deal.II version Git a0b41b6d0f 2020-02-26 20:08:13 -0600
\(\newcommand{\vcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\vcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
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 namespace Particles
28 {
29  namespace
30  {
31  template <int dim, int spacedim>
32  std::vector<char>
33  pack_particles(std::vector<Particle<dim, spacedim>> &particles)
34  {
35  std::vector<char> buffer;
36 
37  if (particles.size() == 0)
38  return buffer;
39 
40  buffer.resize(particles.size() *
41  particles.front().serialized_size_in_bytes());
42  void *current_data = buffer.data();
43 
44  for (const auto &particle : particles)
45  {
46  particle.write_data(current_data);
47  }
48 
49  return buffer;
50  }
51 
52  template <int dim, int spacedim>
53  std::vector<Particle<dim, spacedim>>
54  unpack_particles(
55  const boost::iterator_range<std::vector<char>::const_iterator>
56  & data_range,
57  PropertyPool &property_pool)
58  {
59  std::vector<Particle<dim, spacedim>> particles;
60 
61  if (data_range.empty())
62  return particles;
63 
64  Particle<dim, spacedim> particle;
65  particle.set_property_pool(property_pool);
66  const unsigned int particle_size = particle.serialized_size_in_bytes();
67 
68  particles.reserve(data_range.size() / particle_size);
69 
70  const void *data = static_cast<const void *>(&(*data_range.begin()));
71 
72  while (data < &(*data_range.end()))
73  {
74  particles.emplace_back(data, &property_pool);
75  }
76 
77  Assert(
78  data == &(*data_range.end()),
79  ExcMessage(
80  "The particle data could not be deserialized successfully. "
81  "Check that when deserializing the particles you expect the same "
82  "number of properties that were serialized."));
83 
84  return particles;
85  }
86  } // namespace
87 
88  template <int dim, int spacedim>
90  : triangulation()
91  , particles()
92  , ghost_particles()
93  , global_number_of_particles(0)
94  , global_max_particles_per_cell(0)
95  , next_free_particle_index(0)
96  , property_pool(new PropertyPool(0))
97  , size_callback()
98  , store_callback()
99  , load_callback()
100  , handle(numbers::invalid_unsigned_int)
101  {}
102 
103 
104 
105  template <int dim, int spacedim>
109  const unsigned int n_properties)
110  : triangulation(&triangulation, typeid(*this).name())
111  , mapping(&mapping, typeid(*this).name())
112  , particles()
113  , ghost_particles()
117  , property_pool(new PropertyPool(n_properties))
118  , size_callback()
119  , store_callback()
120  , load_callback()
121  , handle(numbers::invalid_unsigned_int)
122  {}
123 
124 
125 
126  template <int dim, int spacedim>
127  void
129  const Triangulation<dim, spacedim> &new_triangulation,
130  const Mapping<dim, spacedim> & new_mapping,
131  const unsigned int n_properties)
132  {
133  triangulation = &new_triangulation;
134  mapping = &new_mapping;
135 
136  // Create the memory pool that will store all particle properties
137  property_pool = std_cxx14::make_unique<PropertyPool>(n_properties);
138  }
139 
140 
141 
142  template <int dim, int spacedim>
143  void
145  {
146  clear_particles();
150  }
151 
152 
153 
154  template <int dim, int spacedim>
155  void
157  {
158  particles.clear();
159  }
160 
161 
162 
163  template <int dim, int spacedim>
164  void
166  {
167  types::particle_index locally_highest_index = 0;
168  unsigned int local_max_particles_per_cell = 0;
169  unsigned int current_particles_per_cell = 0;
171  triangulation->begin_active();
172 
173  for (particle_iterator particle = begin(); particle != end(); ++particle)
174  {
175  locally_highest_index =
176  std::max(locally_highest_index, particle->get_id());
177 
178  if (particle->get_surrounding_cell(*triangulation) != current_cell)
179  {
180  current_particles_per_cell = 0;
181  current_cell = particle->get_surrounding_cell(*triangulation);
182  }
183 
184  ++current_particles_per_cell;
185  local_max_particles_per_cell =
186  std::max(local_max_particles_per_cell, current_particles_per_cell);
187  }
188 
189  if (const auto parallel_triangulation =
190  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
191  &*triangulation))
192  {
193  global_number_of_particles = ::Utilities::MPI::sum(
194  particles.size(), parallel_triangulation->get_communicator());
197  0 :
198  ::Utilities::MPI::max(
199  locally_highest_index,
200  parallel_triangulation->get_communicator()) +
201  1;
202  global_max_particles_per_cell = ::Utilities::MPI::max(
203  local_max_particles_per_cell,
204  parallel_triangulation->get_communicator());
205  }
206  else
207  {
208  global_number_of_particles = particles.size();
210  global_number_of_particles == 0 ? 0 : locally_highest_index + 1;
211  global_max_particles_per_cell = local_max_particles_per_cell;
212  }
213  }
214 
215 
216 
217  template <int dim, int spacedim>
220  {
221  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin();
222  }
223 
224 
225 
226  template <int dim, int spacedim>
229  {
230  return particle_iterator(particles, particles.begin());
231  }
232 
233 
234 
235  template <int dim, int spacedim>
238  {
239  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end();
240  }
241 
242 
243 
244  template <int dim, int spacedim>
247  {
248  return particle_iterator(particles, particles.end());
249  }
250 
251 
252 
253  template <int dim, int spacedim>
256  {
257  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin_ghost();
258  }
259 
260 
261 
262  template <int dim, int spacedim>
265  {
267  }
268 
269 
270 
271  template <int dim, int spacedim>
274  {
275  return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end_ghost();
276  }
277 
278 
279 
280  template <int dim, int spacedim>
283  {
285  }
286 
287 
288 
289  template <int dim, int spacedim>
293  const
294  {
295  return (const_cast<ParticleHandler<dim, spacedim> *>(this))
296  ->particles_in_cell(cell);
297  }
298 
299 
300 
301  template <int dim, int spacedim>
305  {
306  const internal::LevelInd level_index =
307  std::make_pair(cell->level(), cell->index());
308 
309  if (cell->is_ghost())
310  {
311  const auto particles_in_cell = ghost_particles.equal_range(level_index);
312  return boost::make_iterator_range(
315  }
316 
317  const auto particles_in_cell = particles.equal_range(level_index);
318  return boost::make_iterator_range(
319  particle_iterator(particles, particles_in_cell.first),
320  particle_iterator(particles, particles_in_cell.second));
321  }
322 
323 
324 
325  template <int dim, int spacedim>
326  void
329  {
330  particles.erase(particle->particle);
331  }
332 
333 
334 
335  template <int dim, int spacedim>
338  const Particle<dim, spacedim> & particle,
340  {
341  typename std::multimap<internal::LevelInd,
342  Particle<dim, spacedim>>::iterator it =
343  particles.insert(
344  std::make_pair(internal::LevelInd(cell->level(), cell->index()),
345  particle));
346 
347  particle_iterator particle_it(particles, it);
348  particle_it->set_property_pool(*property_pool);
349 
350  if (particle.has_properties())
351  for (unsigned int n = 0; n < particle.get_properties().size(); ++n)
352  particle_it->get_properties()[n] = particle.get_properties()[n];
353 
354  return particle_it;
355  }
356 
357 
358 
359  template <int dim, int spacedim>
360  void
362  const std::multimap<
364  Particle<dim, spacedim>> &new_particles)
365  {
366  for (auto particle = new_particles.begin(); particle != new_particles.end();
367  ++particle)
368  particles.insert(
369  particles.end(),
370  std::make_pair(internal::LevelInd(particle->first->level(),
371  particle->first->index()),
372  particle->second));
373 
375  }
376 
377 
378 
379  template <int dim, int spacedim>
380  void
382  const std::vector<Point<spacedim>> &positions)
383  {
385 
386  // Determine the starting particle index of this process, which
387  // is the highest currently existing particle index plus the sum
388  // of the number of newly generated particles of all
389  // processes with a lower rank if in a parallel computation.
390  const types::particle_index local_next_particle_index =
392  types::particle_index local_start_index = 0;
393 
394 #ifdef DEAL_II_WITH_MPI
395  if (const auto parallel_triangulation =
396  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
397  &*triangulation))
398  {
399  types::particle_index particles_to_add_locally = positions.size();
400  const int ierr = MPI_Scan(&particles_to_add_locally,
401  &local_start_index,
402  1,
403  DEAL_II_PARTICLE_INDEX_MPI_TYPE,
404  MPI_SUM,
405  parallel_triangulation->get_communicator());
406  AssertThrowMPI(ierr);
407  local_start_index -= particles_to_add_locally;
408  }
409 #endif
410 
411  local_start_index += local_next_particle_index;
412 
414  auto point_locations = GridTools::compute_point_locations(cache, positions);
415 
416  auto &cells = std::get<0>(point_locations);
417  auto &local_positions = std::get<1>(point_locations);
418  auto &index_map = std::get<2>(point_locations);
419 
420  if (cells.size() == 0)
421  return;
422 
423  auto hint =
424  particles.find(std::make_pair(cells[0]->level(), cells[0]->index()));
425  for (unsigned int i = 0; i < cells.size(); ++i)
426  {
427  internal::LevelInd current_cell(cells[i]->level(), cells[i]->index());
428  for (unsigned int p = 0; p < local_positions[i].size(); ++p)
429  {
430  hint = particles.insert(
431  hint,
432  std::make_pair(current_cell,
433  Particle<dim, spacedim>(positions[index_map[i][p]],
434  local_positions[i][p],
435  local_start_index +
436  index_map[i][p])));
437  }
438  }
439 
441  }
442 
443 
444 
445  template <int dim, int spacedim>
446  std::map<unsigned int, IndexSet>
448  const std::vector<Point<spacedim>> &positions,
449  const std::vector<std::vector<BoundingBox<spacedim>>>
450  & global_bounding_boxes,
451  const std::vector<std::vector<double>> &properties)
452  {
453  if (!properties.empty())
454  {
455  AssertDimension(properties.size(), positions.size());
456 #ifdef DEBUG
457  for (const auto &p : properties)
459 #endif
460  }
461 
462  const auto tria =
464  &(*triangulation));
465  const auto comm =
466  (tria != nullptr ? tria->get_communicator() : MPI_COMM_WORLD);
467 
468  const auto n_mpi_processes = Utilities::MPI::n_mpi_processes(comm);
469 
471 
472  // Compute the global number of properties
473  const auto n_global_properties =
474  Utilities::MPI::sum(properties.size(), comm);
475 
476  // Gather the number of points per processor
477  const auto n_particles_per_proc =
478  Utilities::MPI::all_gather(comm, positions.size());
479 
480  // Calculate all starting points locally
481  std::vector<unsigned int> particle_start_indices(n_mpi_processes);
482 
483  unsigned int particle_start_index = 0;
484  for (unsigned int process = 0; process < particle_start_indices.size();
485  ++process)
486  {
487  particle_start_indices[process] = particle_start_index;
488  particle_start_index += n_particles_per_proc[process];
489  }
490 
491  // Get all local information
492  const auto cells_positions_and_index_maps =
494  positions,
495  global_bounding_boxes);
496 
497  // Unpack the information into several vectors:
498  // All cells that contain at least one particle
499  const auto &local_cells_containing_particles =
500  std::get<0>(cells_positions_and_index_maps);
501 
502  // The reference position of every particle in the local part of the
503  // triangulation.
504  const auto &local_reference_positions =
505  std::get<1>(cells_positions_and_index_maps);
506  // The original index in the positions vector for each particle in the
507  // local part of the triangulation
508  const auto &original_indices_of_local_particles =
509  std::get<2>(cells_positions_and_index_maps);
510  // The real spatial position of every particle in the local part of the
511  // triangulation.
512  const auto &local_positions = std::get<3>(cells_positions_and_index_maps);
513  // The MPI process that inserted each particle
514  const auto &calling_process_indices =
515  std::get<4>(cells_positions_and_index_maps);
516 
517  // Create the map of cpu to indices, indicating who sent us what
518  // particle.
519  std::map<unsigned int, IndexSet> original_process_to_local_particle_indices;
520  for (unsigned int i_cell = 0;
521  i_cell < local_cells_containing_particles.size();
522  ++i_cell)
523  {
524  for (unsigned int i_particle = 0;
525  i_particle < local_positions[i_cell].size();
526  ++i_particle)
527  {
528  const unsigned int local_id_on_calling_process =
529  original_indices_of_local_particles[i_cell][i_particle];
530  const unsigned int calling_process =
531  calling_process_indices[i_cell][i_particle];
532 
533  if (original_process_to_local_particle_indices.find(
534  calling_process) ==
535  original_process_to_local_particle_indices.end())
536  original_process_to_local_particle_indices.insert(
537  {calling_process,
538  IndexSet(n_particles_per_proc[calling_process])});
539 
540  original_process_to_local_particle_indices[calling_process]
541  .add_index(local_id_on_calling_process);
542  }
543  }
544  for (auto &process_and_particle_indices :
545  original_process_to_local_particle_indices)
546  process_and_particle_indices.second.compress();
547 
548 
549  // A map from mpi process to properties, ordered as in the IndexSet.
550  // Notice that this ordering maybe different from the ordering in the
551  // vectors above, since no local ordering is guaranteed by the
552  // distribute_compute_point_locations() call.
553  // This is only filled if n_global_properties is > 0
554  std::map<unsigned int, std::vector<std::vector<double>>>
555  locally_owned_properties_from_other_processes;
556 
557  if (n_global_properties > 0)
558  {
559  // Gather whom I sent my own particles to, to decide whom to send
560  // the particle properties
561  auto send_to_cpu = Utilities::MPI::some_to_some(
562  comm, original_process_to_local_particle_indices);
563 
564  std::map<unsigned int, std::vector<std::vector<double>>>
565  non_locally_owned_properties;
566 
567  // Prepare the vector of properties to send
568  for (const auto &it : send_to_cpu)
569  {
570  std::vector<std::vector<double>> properties_to_send(
571  it.second.n_elements(),
572  std::vector<double>(n_properties_per_particle()));
573  unsigned int index = 0;
574  for (const auto &el : it.second)
575  properties_to_send[index++] = properties[el];
576  non_locally_owned_properties.insert({it.first, properties_to_send});
577  }
578 
579  // Send the non locally owned properties to each mpi process
580  // that needs them
581  locally_owned_properties_from_other_processes =
582  Utilities::MPI::some_to_some(comm, non_locally_owned_properties);
583 
584  AssertDimension(locally_owned_properties_from_other_processes.size(),
585  original_process_to_local_particle_indices.size());
586  }
587 
588 
589  // Create the multimap of local particles
590  std::multimap<typename Triangulation<dim, spacedim>::active_cell_iterator,
592  particles;
593 
594  // Now fill up the actual particles
595  for (unsigned int i_cell = 0;
596  i_cell < local_cells_containing_particles.size();
597  ++i_cell)
598  {
599  for (unsigned int i_particle = 0;
600  i_particle < local_positions[i_cell].size();
601  ++i_particle)
602  {
603  const unsigned int local_id_on_calling_process =
604  original_indices_of_local_particles[i_cell][i_particle];
605  const unsigned int calling_process =
606  calling_process_indices[i_cell][i_particle];
607 
608  const unsigned int particle_id =
609  local_id_on_calling_process +
610  particle_start_indices[calling_process];
611 
612  Particle<dim, spacedim> particle(
613  local_positions[i_cell][i_particle],
614  local_reference_positions[i_cell][i_particle],
615  particle_id);
616 
617  if (n_global_properties > 0)
618  {
619  const unsigned int index_within_set =
620  original_process_to_local_particle_indices[calling_process]
621  .index_within_set(local_id_on_calling_process);
622 
623  const auto &this_particle_properties =
624  locally_owned_properties_from_other_processes
625  [calling_process][index_within_set];
626 
628  particle.set_properties(this_particle_properties);
629  }
630 
631  particles.emplace(local_cells_containing_particles[i_cell],
632  particle);
633  }
634  }
635 
636  this->insert_particles(particles);
637 
638  return original_process_to_local_particle_indices;
639  }
640 
641 
642 
643  template <int dim, int spacedim>
646  {
648  }
649 
650 
651 
652  template <int dim, int spacedim>
655  {
657  }
658 
659 
660 
661  template <int dim, int spacedim>
664  {
665  return particles.size();
666  }
667 
668 
669 
670  template <int dim, int spacedim>
671  unsigned int
673  {
674  return property_pool->n_properties_per_slot();
675  }
676 
677 
678 
679  template <int dim, int spacedim>
680  unsigned int
683  const
684  {
685  const internal::LevelInd found_cell =
686  std::make_pair(cell->level(), cell->index());
687 
688  if (cell->is_locally_owned())
689  return particles.count(found_cell);
690  else if (cell->is_ghost())
691  return ghost_particles.count(found_cell);
692  else if (cell->is_artificial())
693  AssertThrow(false, ExcInternalError());
694 
695  return 0;
696  }
697 
698 
699 
700  template <int dim, int spacedim>
703  {
705  }
706 
707 
708 
709  template <int dim, int spacedim>
710  IndexSet
712  {
714  for (const auto &p : *this)
715  set.add_index(p.get_id());
716  set.compress();
717  return set;
718  }
719 
720 
721 
722  template <int dim, int spacedim>
723  void
725  std::vector<Point<spacedim>> &positions,
726  const bool add_to_output_vector)
727  {
728  // There should be one point per particle to gather
729  AssertDimension(positions.size(), n_locally_owned_particles());
730 
731  unsigned int i = 0;
732  for (auto it = begin(); it != end(); ++it, ++i)
733  {
734  if (add_to_output_vector)
735  positions[i] = positions[i] + it->get_location();
736  else
737  positions[i] = it->get_location();
738  }
739  }
740 
741 
742 
743  template <int dim, int spacedim>
744  void
746  const std::vector<Point<spacedim>> &new_positions,
747  const bool displace_particles)
748  {
749  // There should be one point per particle to fix the new position
750  AssertDimension(new_positions.size(), n_locally_owned_particles());
751 
752  unsigned int i = 0;
753  for (auto it = begin(); it != end(); ++it, ++i)
754  {
755  if (displace_particles)
756  it->set_location(it->get_location() + new_positions[i]);
757  else
758  it->set_location(new_positions[i]);
759  }
761  }
762 
763 
764 
765  template <int dim, int spacedim>
766  void
768  const Function<spacedim> &function,
769  const bool displace_particles)
770  {
771  // The function should have sufficient components to displace the
772  // particles
773  AssertDimension(function.n_components, spacedim);
774 
775  Vector<double> new_position(spacedim);
776  for (auto &particle : *this)
777  {
778  Point<spacedim> particle_location = particle.get_location();
779  function.vector_value(particle_location, new_position);
780  if (displace_particles)
781  for (unsigned int d = 0; d < spacedim; ++d)
782  particle_location[d] += new_position[d];
783  else
784  for (unsigned int d = 0; d < spacedim; ++d)
785  particle_location[d] = new_position[d];
786  particle.set_location(particle_location);
787  }
789  }
790 
791 
792 
793  template <int dim, int spacedim>
794  PropertyPool &
796  {
797  return *property_pool;
798  }
799 
800 
801 
802  namespace
803  {
813  template <int dim>
814  bool
815  compare_particle_association(
816  const unsigned int a,
817  const unsigned int b,
818  const Tensor<1, dim> & particle_direction,
819  const std::vector<Tensor<1, dim>> &center_directions)
820  {
821  const double scalar_product_a = center_directions[a] * particle_direction;
822  const double scalar_product_b = center_directions[b] * particle_direction;
823 
824  // The function is supposed to return if a is before b. We are looking
825  // for the alignment of particle direction and center direction,
826  // therefore return if the scalar product of a is larger.
827  return (scalar_product_a > scalar_product_b);
828  }
829  } // namespace
830 
831 
832 
833  template <int dim, int spacedim>
834  void
836  {
837  // TODO: The current algorithm only works for particles that are in
838  // the local domain or in ghost cells, because it only knows the
839  // subdomain_id of ghost cells, but not of artificial cells. This
840  // can be extended using the distributed version of compute point
841  // locations.
842  // TODO: Extend this function to allow keeping particles on other
843  // processes around (with an invalid cell).
844 
845  std::vector<particle_iterator> particles_out_of_cell;
846  particles_out_of_cell.reserve(n_locally_owned_particles());
847 
848  // Now update the reference locations of the moved particles
849  for (particle_iterator it = begin(); it != end(); ++it)
850  {
851  const typename Triangulation<dim, spacedim>::cell_iterator cell =
852  it->get_surrounding_cell(*triangulation);
853 
854  try
855  {
856  const Point<dim> p_unit =
857  mapping->transform_real_to_unit_cell(cell, it->get_location());
859  {
860  it->set_reference_location(p_unit);
861  }
862  else
863  {
864  // The particle has left the cell
865  particles_out_of_cell.push_back(it);
866  }
867  }
868  catch (typename Mapping<dim>::ExcTransformationFailed &)
869  {
870  // The particle has left the cell
871  particles_out_of_cell.push_back(it);
872  }
873  }
874 
875  // There are three reasons why a particle is not in its old cell:
876  // It moved to another cell, to another subdomain or it left the mesh.
877  // Particles that moved to another cell are updated and stored inside the
878  // sorted_particles vector, particles that moved to another domain are
879  // collected in the moved_particles_domain vector. Particles that left
880  // the mesh completely are ignored and removed.
881  std::vector<std::pair<internal::LevelInd, Particle<dim, spacedim>>>
882  sorted_particles;
883  std::map<types::subdomain_id, std::vector<particle_iterator>>
884  moved_particles;
885  std::map<
887  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
888  moved_cells;
889 
890  // We do not know exactly how many particles are lost, exchanged between
891  // domains, or remain on this process. Therefore we pre-allocate
892  // approximate sizes for these vectors. If more space is needed an
893  // automatic and relatively fast (compared to other parts of this
894  // algorithm) re-allocation will happen.
895  using vector_size = typename std::vector<particle_iterator>::size_type;
896  sorted_particles.reserve(
897  static_cast<vector_size>(particles_out_of_cell.size() * 1.25));
898 
899  std::set<types::subdomain_id> ghost_owners;
900  if (const auto parallel_triangulation =
901  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
902  &*triangulation))
903  ghost_owners = parallel_triangulation->ghost_owners();
904 
905  for (const auto ghost_owner : ghost_owners)
906  moved_particles[ghost_owner].reserve(
907  static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
908  for (const auto ghost_owner : ghost_owners)
909  moved_cells[ghost_owner].reserve(
910  static_cast<vector_size>(particles_out_of_cell.size() * 0.25));
911 
912  {
913  // Create a map from vertices to adjacent cells
914  const std::vector<
915  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
916  vertex_to_cells(GridTools::vertex_to_cell_map(*triangulation));
917 
918  // Create a corresponding map of vectors from vertex to cell center
919  const std::vector<std::vector<Tensor<1, spacedim>>>
920  vertex_to_cell_centers(
922  vertex_to_cells));
923 
924  std::vector<unsigned int> neighbor_permutation;
925 
926  // Find the cells that the particles moved to.
927  typename std::vector<particle_iterator>::iterator
928  it = particles_out_of_cell.begin(),
929  end_particle = particles_out_of_cell.end();
930 
931  for (; it != end_particle; ++it)
932  {
933  // The cell the particle is in
934  Point<dim> current_reference_position;
935  bool found_cell = false;
936 
937  // Check if the particle is in one of the old cell's neighbors
938  // that are adjacent to the closest vertex
940  current_cell = (*it)->get_surrounding_cell(*triangulation);
941 
942  const unsigned int closest_vertex =
943  GridTools::find_closest_vertex_of_cell<dim, spacedim>(
944  current_cell, (*it)->get_location());
945  Tensor<1, spacedim> vertex_to_particle =
946  (*it)->get_location() - current_cell->vertex(closest_vertex);
947  vertex_to_particle /= vertex_to_particle.norm();
948 
949  const unsigned int closest_vertex_index =
950  current_cell->vertex_index(closest_vertex);
951  const unsigned int n_neighbor_cells =
952  vertex_to_cells[closest_vertex_index].size();
953 
954  neighbor_permutation.resize(n_neighbor_cells);
955  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
956  neighbor_permutation[i] = i;
957 
958  const auto cell_centers =
959  vertex_to_cell_centers[closest_vertex_index];
960  std::sort(neighbor_permutation.begin(),
961  neighbor_permutation.end(),
962  [&vertex_to_particle, &cell_centers](const unsigned int a,
963  const unsigned int b) {
964  return compare_particle_association(a,
965  b,
966  vertex_to_particle,
967  cell_centers);
968  });
969 
970  // Search all of the cells adjacent to the closest vertex of the
971  // previous cell Most likely we will find the particle in them.
972  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
973  {
974  try
975  {
976  typename std::set<typename Triangulation<dim, spacedim>::
977  active_cell_iterator>::const_iterator
978  cell = vertex_to_cells[closest_vertex_index].begin();
979 
980  std::advance(cell, neighbor_permutation[i]);
981  const Point<dim> p_unit =
982  mapping->transform_real_to_unit_cell(*cell,
983  (*it)->get_location());
985  {
986  current_cell = *cell;
987  current_reference_position = p_unit;
988  found_cell = true;
989  break;
990  }
991  }
992  catch (typename Mapping<dim>::ExcTransformationFailed &)
993  {}
994  }
995 
996  if (!found_cell)
997  {
998  // The particle is not in a neighbor of the old cell.
999  // Look for the new cell in the whole local domain.
1000  // This case is rare.
1001  try
1002  {
1003  const std::pair<const typename Triangulation<dim, spacedim>::
1004  active_cell_iterator,
1005  Point<dim>>
1006  current_cell_and_position =
1007  GridTools::find_active_cell_around_point<>(
1008  *mapping, *triangulation, (*it)->get_location());
1009  current_cell = current_cell_and_position.first;
1010  current_reference_position = current_cell_and_position.second;
1011  }
1012  catch (GridTools::ExcPointNotFound<spacedim> &)
1013  {
1014  // We can find no cell for this particle. It has left the
1015  // domain due to an integration error or an open boundary.
1016  continue;
1017  }
1018  }
1019 
1020  // If we are here, we found a cell and reference position for this
1021  // particle
1022  (*it)->set_reference_location(current_reference_position);
1023 
1024  // Reinsert the particle into our domain if we own its cell.
1025  // Mark it for MPI transfer otherwise
1026  if (current_cell->is_locally_owned())
1027  {
1028  sorted_particles.push_back(
1029  std::make_pair(internal::LevelInd(current_cell->level(),
1030  current_cell->index()),
1031  (*it)->particle->second));
1032  }
1033  else
1034  {
1035  moved_particles[current_cell->subdomain_id()].push_back(*it);
1036  moved_cells[current_cell->subdomain_id()].push_back(current_cell);
1037  }
1038  }
1039  }
1040 
1041  // Sort the updated particles. This pre-sort speeds up inserting
1042  // them into particles to O(N) complexity.
1043  std::multimap<internal::LevelInd, Particle<dim, spacedim>>
1044  sorted_particles_map;
1045 
1046  // Exchange particles between processors if we have more than one process
1047 #ifdef DEAL_II_WITH_MPI
1048  if (const auto parallel_triangulation =
1049  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
1050  &*triangulation))
1051  {
1053  parallel_triangulation->get_communicator()) > 1)
1054  send_recv_particles(moved_particles,
1055  sorted_particles_map,
1056  moved_cells);
1057  }
1058 #endif
1059 
1060  sorted_particles_map.insert(sorted_particles.begin(),
1061  sorted_particles.end());
1062 
1063  for (unsigned int i = 0; i < particles_out_of_cell.size(); ++i)
1064  remove_particle(particles_out_of_cell[i]);
1065 
1066  particles.insert(sorted_particles_map.begin(), sorted_particles_map.end());
1068  }
1069 
1070 
1071 
1072  template <int dim, int spacedim>
1073  void
1075  {
1076  // Nothing to do in serial computations
1077  const auto parallel_triangulation =
1078  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
1079  &*triangulation);
1080  if (parallel_triangulation != nullptr)
1081  {
1083  parallel_triangulation->get_communicator()) == 1)
1084  return;
1085  }
1086  else
1087  return;
1088 
1089 #ifdef DEAL_II_WITH_MPI
1090  // First clear the current ghost_particle information
1091  ghost_particles.clear();
1092 
1093  std::map<types::subdomain_id, std::vector<particle_iterator>>
1094  ghost_particles_by_domain;
1095 
1096  const std::set<types::subdomain_id> ghost_owners =
1097  parallel_triangulation->ghost_owners();
1098  for (const auto ghost_owner : ghost_owners)
1099  ghost_particles_by_domain[ghost_owner].reserve(
1100  static_cast<typename std::vector<particle_iterator>::size_type>(
1101  particles.size() * 0.25));
1102 
1103  std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain(
1104  triangulation->n_vertices());
1105 
1106  for (const auto &cell : triangulation->active_cell_iterators())
1107  {
1108  if (cell->is_ghost())
1109  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
1110  ++v)
1111  vertex_to_neighbor_subdomain[cell->vertex_index(v)].insert(
1112  cell->subdomain_id());
1113  }
1114 
1115  for (const auto &cell : triangulation->active_cell_iterators())
1116  {
1117  if (!cell->is_ghost())
1118  {
1119  std::set<unsigned int> cell_to_neighbor_subdomain;
1120  for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell;
1121  ++v)
1122  {
1123  cell_to_neighbor_subdomain.insert(
1124  vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
1125  vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
1126  }
1127 
1128  if (cell_to_neighbor_subdomain.size() > 0)
1129  {
1130  const particle_iterator_range particle_range =
1131  particles_in_cell(cell);
1132 
1133  for (const auto domain : cell_to_neighbor_subdomain)
1134  {
1135  for (typename particle_iterator_range::iterator particle =
1136  particle_range.begin();
1137  particle != particle_range.end();
1138  ++particle)
1139  ghost_particles_by_domain[domain].push_back(particle);
1140  }
1141  }
1142  }
1143  }
1144 
1145  send_recv_particles(ghost_particles_by_domain, ghost_particles);
1146 #endif
1147  }
1148 
1149 
1150 
1151 #ifdef DEAL_II_WITH_MPI
1152  template <int dim, int spacedim>
1153  void
1155  const std::map<types::subdomain_id, std::vector<particle_iterator>>
1156  &particles_to_send,
1157  std::multimap<internal::LevelInd, Particle<dim, spacedim>>
1158  &received_particles,
1159  const std::map<
1162  &send_cells)
1163  {
1164  const auto parallel_triangulation =
1165  dynamic_cast<const parallel::Triangulation<dim, spacedim> *>(
1166  &*triangulation);
1167  Assert(
1168  parallel_triangulation,
1169  ExcMessage(
1170  "This function is only implemented for parallel::Triangulation objects."));
1171 
1172  // Determine the communication pattern
1173  const std::set<types::subdomain_id> ghost_owners =
1174  parallel_triangulation->ghost_owners();
1175  const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
1176  ghost_owners.end());
1177  const unsigned int n_neighbors = neighbors.size();
1178 
1179  if (send_cells.size() != 0)
1180  Assert(particles_to_send.size() == send_cells.size(), ExcInternalError());
1181 
1182  // If we do not know the subdomain this particle needs to be send to,
1183  // throw an error
1184  Assert(particles_to_send.find(numbers::artificial_subdomain_id) ==
1185  particles_to_send.end(),
1186  ExcInternalError());
1187 
1188  // TODO: Implement the shipping of particles to processes that are not
1189  // ghost owners of the local domain
1190  for (auto send_particles = particles_to_send.begin();
1191  send_particles != particles_to_send.end();
1192  ++send_particles)
1193  Assert(ghost_owners.find(send_particles->first) != ghost_owners.end(),
1194  ExcNotImplemented());
1195 
1196  unsigned int n_send_particles = 0;
1197  for (auto send_particles = particles_to_send.begin();
1198  send_particles != particles_to_send.end();
1199  ++send_particles)
1200  n_send_particles += send_particles->second.size();
1201 
1202  const unsigned int cellid_size = sizeof(CellId::binary_type);
1203 
1204  // Containers for the amount and offsets of data we will send
1205  // to other processors and the data itself.
1206  std::vector<unsigned int> n_send_data(n_neighbors, 0);
1207  std::vector<unsigned int> send_offsets(n_neighbors, 0);
1208  std::vector<char> send_data;
1209 
1210  // Only serialize things if there are particles to be send.
1211  // We can not return early even if no particles
1212  // are send, because we might receive particles from other processes
1213  if (n_send_particles > 0)
1214  {
1215  // Allocate space for sending particle data
1216  const unsigned int particle_size =
1217  begin()->serialized_size_in_bytes() + cellid_size +
1218  (size_callback ? size_callback() : 0);
1219  send_data.resize(n_send_particles * particle_size);
1220  void *data = static_cast<void *>(&send_data.front());
1221 
1222  // Serialize the data sorted by receiving process
1223  for (unsigned int i = 0; i < n_neighbors; ++i)
1224  {
1225  send_offsets[i] = reinterpret_cast<std::size_t>(data) -
1226  reinterpret_cast<std::size_t>(&send_data.front());
1227 
1228  for (unsigned int j = 0;
1229  j < particles_to_send.at(neighbors[i]).size();
1230  ++j)
1231  {
1232  // If no target cells are given, use the iterator information
1234  cell;
1235  if (send_cells.size() == 0)
1236  cell =
1237  particles_to_send.at(neighbors[i])[j]->get_surrounding_cell(
1238  *triangulation);
1239  else
1240  cell = send_cells.at(neighbors[i])[j];
1241 
1242  const CellId::binary_type cellid =
1243  cell->id().template to_binary<dim>();
1244  memcpy(data, &cellid, cellid_size);
1245  data = static_cast<char *>(data) + cellid_size;
1246 
1247  particles_to_send.at(neighbors[i])[j]->write_data(data);
1248  if (store_callback)
1249  data =
1250  store_callback(particles_to_send.at(neighbors[i])[j], data);
1251  }
1252  n_send_data[i] = reinterpret_cast<std::size_t>(data) -
1253  send_offsets[i] -
1254  reinterpret_cast<std::size_t>(&send_data.front());
1255  }
1256  }
1257 
1258  // Containers for the data we will receive from other processors
1259  std::vector<unsigned int> n_recv_data(n_neighbors);
1260  std::vector<unsigned int> recv_offsets(n_neighbors);
1261 
1262  // Notify other processors how many particles we will send
1263  {
1264  const int mpi_tag = Utilities::MPI::internal::Tags::
1266 
1267  std::vector<MPI_Request> n_requests(2 * n_neighbors);
1268  for (unsigned int i = 0; i < n_neighbors; ++i)
1269  {
1270  const int ierr = MPI_Irecv(&(n_recv_data[i]),
1271  1,
1272  MPI_UNSIGNED,
1273  neighbors[i],
1274  mpi_tag,
1275  parallel_triangulation->get_communicator(),
1276  &(n_requests[2 * i]));
1277  AssertThrowMPI(ierr);
1278  }
1279  for (unsigned int i = 0; i < n_neighbors; ++i)
1280  {
1281  const int ierr = MPI_Isend(&(n_send_data[i]),
1282  1,
1283  MPI_UNSIGNED,
1284  neighbors[i],
1285  mpi_tag,
1286  parallel_triangulation->get_communicator(),
1287  &(n_requests[2 * i + 1]));
1288  AssertThrowMPI(ierr);
1289  }
1290  const int ierr =
1291  MPI_Waitall(2 * n_neighbors, n_requests.data(), MPI_STATUSES_IGNORE);
1292  AssertThrowMPI(ierr);
1293  }
1294 
1295  // Determine how many particles and data we will receive
1296  unsigned int total_recv_data = 0;
1297  for (unsigned int neighbor_id = 0; neighbor_id < n_neighbors; ++neighbor_id)
1298  {
1299  recv_offsets[neighbor_id] = total_recv_data;
1300  total_recv_data += n_recv_data[neighbor_id];
1301  }
1302 
1303  // Set up the space for the received particle data
1304  std::vector<char> recv_data(total_recv_data);
1305 
1306  // Exchange the particle data between domains
1307  {
1308  std::vector<MPI_Request> requests(2 * n_neighbors);
1309  unsigned int send_ops = 0;
1310  unsigned int recv_ops = 0;
1311 
1312  const int mpi_tag = Utilities::MPI::internal::Tags::
1314 
1315  for (unsigned int i = 0; i < n_neighbors; ++i)
1316  if (n_recv_data[i] > 0)
1317  {
1318  const int ierr =
1319  MPI_Irecv(&(recv_data[recv_offsets[i]]),
1320  n_recv_data[i],
1321  MPI_CHAR,
1322  neighbors[i],
1323  mpi_tag,
1324  parallel_triangulation->get_communicator(),
1325  &(requests[send_ops]));
1326  AssertThrowMPI(ierr);
1327  send_ops++;
1328  }
1329 
1330  for (unsigned int i = 0; i < n_neighbors; ++i)
1331  if (n_send_data[i] > 0)
1332  {
1333  const int ierr =
1334  MPI_Isend(&(send_data[send_offsets[i]]),
1335  n_send_data[i],
1336  MPI_CHAR,
1337  neighbors[i],
1338  mpi_tag,
1339  parallel_triangulation->get_communicator(),
1340  &(requests[send_ops + recv_ops]));
1341  AssertThrowMPI(ierr);
1342  recv_ops++;
1343  }
1344  const int ierr =
1345  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
1346  AssertThrowMPI(ierr);
1347  }
1348 
1349  // Put the received particles into the domain if they are in the
1350  // triangulation
1351  const void *recv_data_it = static_cast<const void *>(recv_data.data());
1352 
1353  while (reinterpret_cast<std::size_t>(recv_data_it) -
1354  reinterpret_cast<std::size_t>(recv_data.data()) <
1355  total_recv_data)
1356  {
1357  CellId::binary_type binary_cellid;
1358  memcpy(&binary_cellid, recv_data_it, cellid_size);
1359  const CellId id(binary_cellid);
1360  recv_data_it = static_cast<const char *>(recv_data_it) + cellid_size;
1361 
1363  id.to_cell(*triangulation);
1364 
1365  typename std::multimap<internal::LevelInd,
1366  Particle<dim, spacedim>>::iterator
1367  recv_particle = received_particles.insert(std::make_pair(
1368  internal::LevelInd(cell->level(), cell->index()),
1369  Particle<dim, spacedim>(recv_data_it, property_pool.get())));
1370 
1371  if (load_callback)
1372  recv_data_it =
1373  load_callback(particle_iterator(received_particles, recv_particle),
1374  recv_data_it);
1375  }
1376 
1377  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
1378  ExcMessage(
1379  "The amount of data that was read into new particles "
1380  "does not match the amount of data sent around."));
1381  }
1382 #endif
1383 
1384 
1385 
1386  template <int dim, int spacedim>
1387  void
1389  const std::function<std::size_t()> & size_callb,
1390  const std::function<void *(const particle_iterator &, void *)> &store_callb,
1391  const std::function<const void *(const particle_iterator &, const void *)>
1392  &load_callb)
1393  {
1394  size_callback = size_callb;
1395  store_callback = store_callb;
1396  load_callback = load_callb;
1397  }
1398 
1399 
1400 
1401  template <int dim, int spacedim>
1402  void
1404  {
1406  *non_const_triangulation =
1409  *>(&(*triangulation)));
1410  (void)non_const_triangulation;
1411 
1412  Assert(non_const_triangulation != nullptr, ::ExcNotImplemented());
1413 
1414 #ifdef DEAL_II_WITH_P4EST
1415  // Only save and load particles if there are any, we might get here for
1416  // example if somebody created a ParticleHandler but generated 0
1417  // particles.
1419 
1421  {
1422  const auto callback_function =
1423  [this](const typename Triangulation<dim, spacedim>::cell_iterator
1424  &cell_iterator,
1426  cell_status) {
1427  return this->store_particles(cell_iterator, cell_status);
1428  };
1429 
1430  handle = non_const_triangulation->register_data_attach(
1431  callback_function, /*returns_variable_size_data=*/true);
1432  }
1433 #endif
1434  }
1435 
1436 
1437 
1438  template <int dim, int spacedim>
1439  void
1441  const bool serialization)
1442  {
1443  // All particles have been stored, when we reach this point. Empty the
1444  // particle data.
1445  clear_particles();
1446 
1448  *non_const_triangulation =
1451  *>(&(*triangulation)));
1452  (void)non_const_triangulation;
1453 
1454  Assert(non_const_triangulation != nullptr, ::ExcNotImplemented());
1455 
1456 #ifdef DEAL_II_WITH_P4EST
1457  // If we are resuming from a checkpoint, we first have to register the
1458  // store function again, to set the triangulation in the same state as
1459  // before the serialization. Only by this it knows how to deserialize the
1460  // data correctly. Only do this if something was actually stored.
1461  if (serialization && (global_max_particles_per_cell > 0))
1462  {
1463  const auto callback_function =
1464  [this](const typename Triangulation<dim, spacedim>::cell_iterator
1465  &cell_iterator,
1467  cell_status) {
1468  return this->store_particles(cell_iterator, cell_status);
1469  };
1470 
1471  handle = non_const_triangulation->register_data_attach(
1472  callback_function, /*returns_variable_size_data=*/true);
1473  }
1474 
1475  // Check if something was stored and load it
1477  {
1478  const auto callback_function =
1479  [this](
1481  &cell_iterator,
1482  const typename Triangulation<dim, spacedim>::CellStatus cell_status,
1483  const boost::iterator_range<std::vector<char>::const_iterator>
1484  &range_iterator) {
1485  this->load_particles(cell_iterator, cell_status, range_iterator);
1486  };
1487 
1488  non_const_triangulation->notify_ready_to_unpack(handle,
1489  callback_function);
1490 
1491  // Reset handle and update global number of particles. The number
1492  // can change because of discarded or newly generated particles
1495  }
1496 #else
1497  (void)serialization;
1498 #endif
1499  }
1500 
1501 
1502 
1503  template <int dim, int spacedim>
1504  std::vector<char>
1506  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
1507  const typename Triangulation<dim, spacedim>::CellStatus status) const
1508  {
1509  std::vector<Particle<dim, spacedim>> stored_particles_on_cell;
1510 
1511  switch (status)
1512  {
1515  // If the cell persist or is refined store all particles of the
1516  // current cell.
1517  {
1518  unsigned int n_particles = 0;
1519 
1520  const internal::LevelInd level_index = {cell->level(),
1521  cell->index()};
1522  const auto particles_in_cell =
1523  (cell->is_ghost() ? ghost_particles.equal_range(level_index) :
1524  particles.equal_range(level_index));
1525 
1526  n_particles = n_particles_in_cell(cell);
1527  stored_particles_on_cell.reserve(n_particles);
1528 
1529  std::for_each(
1530  particles_in_cell.first,
1531  particles_in_cell.second,
1532  [&stored_particles_on_cell](
1533  const std::pair<internal::LevelInd, Particle<dim, spacedim>>
1534  &particle) {
1535  stored_particles_on_cell.push_back(particle.second);
1536  });
1537 
1538  AssertDimension(n_particles, stored_particles_on_cell.size());
1539  }
1540  break;
1541 
1543  // If this cell is the parent of children that will be coarsened,
1544  // collect the particles of all children.
1545  {
1546  unsigned int n_particles = 0;
1547 
1548  for (unsigned int child_index = 0;
1549  child_index < GeometryInfo<dim>::max_children_per_cell;
1550  ++child_index)
1551  {
1553  child = cell->child(child_index);
1554  n_particles += n_particles_in_cell(child);
1555  }
1556 
1557  stored_particles_on_cell.reserve(n_particles);
1558 
1559  for (unsigned int child_index = 0;
1560  child_index < GeometryInfo<dim>::max_children_per_cell;
1561  ++child_index)
1562  {
1564  child = cell->child(child_index);
1565  const internal::LevelInd level_index = {child->level(),
1566  child->index()};
1567  const auto particles_in_cell =
1568  (child->is_ghost() ?
1569  ghost_particles.equal_range(level_index) :
1570  particles.equal_range(level_index));
1571 
1572  std::for_each(
1573  particles_in_cell.first,
1574  particles_in_cell.second,
1575  [&stored_particles_on_cell](
1576  const std::pair<internal::LevelInd, Particle<dim, spacedim>>
1577  &particle) {
1578  stored_particles_on_cell.push_back(particle.second);
1579  });
1580  }
1581 
1582  AssertDimension(n_particles, stored_particles_on_cell.size());
1583  }
1584  break;
1585 
1586  default:
1587  Assert(false, ExcInternalError());
1588  break;
1589  }
1590 
1591  return pack_particles(stored_particles_on_cell);
1592  }
1593 
1594  template <int dim, int spacedim>
1595  void
1597  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
1598  const typename Triangulation<dim, spacedim>::CellStatus status,
1599  const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
1600  {
1601  // We leave this container non-const to be able to `std::move`
1602  // its contents directly into the particles multimap later.
1603  std::vector<Particle<dim, spacedim>> loaded_particles_on_cell =
1604  unpack_particles<dim, spacedim>(data_range, *property_pool);
1605 
1606  // Update the reference to the current property pool for all particles.
1607  // This was not stored, because they might be transported across process
1608  // domains.
1609  for (auto &particle : loaded_particles_on_cell)
1610  particle.set_property_pool(*property_pool);
1611 
1612  switch (status)
1613  {
1615  {
1616  auto position_hint = particles.end();
1617  for (const auto &particle : loaded_particles_on_cell)
1618  {
1619  // Use std::multimap::emplace_hint to speed up insertion of
1620  // particles. This is a C++11 function, but not all compilers
1621  // that report a -std=c++11 (like gcc 4.6) implement it, so
1622  // require C++14 instead.
1623 #ifdef DEAL_II_WITH_CXX14
1624  position_hint =
1625  particles.emplace_hint(position_hint,
1626  std::make_pair(cell->level(),
1627  cell->index()),
1628  std::move(particle));
1629 #else
1630  position_hint =
1631  particles.insert(position_hint,
1632  std::make_pair(std::make_pair(cell->level(),
1633  cell->index()),
1634  std::move(particle)));
1635 #endif
1636  // Move the hint position forward by one, i.e., for the next
1637  // particle. The 'hint' position will thus be right after the
1638  // one just inserted.
1639  ++position_hint;
1640  }
1641  }
1642  break;
1643 
1645  {
1646  typename std::multimap<internal::LevelInd,
1647  Particle<dim, spacedim>>::iterator
1648  position_hint = particles.end();
1649  for (auto &particle : loaded_particles_on_cell)
1650  {
1651  const Point<dim> p_unit =
1652  mapping->transform_real_to_unit_cell(cell,
1653  particle.get_location());
1654  particle.set_reference_location(p_unit);
1655  // Use std::multimap::emplace_hint to speed up insertion of
1656  // particles. This is a C++11 function, but not all compilers
1657  // that report a -std=c++11 (like gcc 4.6) implement it, so
1658  // require C++14 instead.
1659 #ifdef DEAL_II_WITH_CXX14
1660  position_hint =
1661  particles.emplace_hint(position_hint,
1662  std::make_pair(cell->level(),
1663  cell->index()),
1664  std::move(particle));
1665 #else
1666  position_hint =
1667  particles.insert(position_hint,
1668  std::make_pair(std::make_pair(cell->level(),
1669  cell->index()),
1670  std::move(particle)));
1671 #endif
1672  // Move the hint position forward by one, i.e., for the next
1673  // particle. The 'hint' position will thus be right after the
1674  // one just inserted.
1675  ++position_hint;
1676  }
1677  }
1678  break;
1679 
1681  {
1682  std::vector<
1683  typename std::multimap<internal::LevelInd,
1684  Particle<dim, spacedim>>::iterator>
1686  for (unsigned int child_index = 0;
1687  child_index < GeometryInfo<dim>::max_children_per_cell;
1688  ++child_index)
1689  {
1691  child = cell->child(child_index);
1692  position_hints[child_index] = particles.upper_bound(
1693  std::make_pair(child->level(), child->index()));
1694  }
1695 
1696  for (auto &particle : loaded_particles_on_cell)
1697  {
1698  for (unsigned int child_index = 0;
1699  child_index < GeometryInfo<dim>::max_children_per_cell;
1700  ++child_index)
1701  {
1703  child = cell->child(child_index);
1704 
1705  try
1706  {
1707  const Point<dim> p_unit =
1708  mapping->transform_real_to_unit_cell(
1709  child, particle.get_location());
1711  {
1712  particle.set_reference_location(p_unit);
1713  // Use std::multimap::emplace_hint to speed up
1714  // insertion of particles. This is a C++11
1715  // function, but not all compilers that report a
1716  // -std=c++11 (like gcc 4.6) implement it, so
1717  // require C++14 instead.
1718 #ifdef DEAL_II_WITH_CXX14
1719  position_hints[child_index] =
1720  particles.emplace_hint(
1721  position_hints[child_index],
1722  std::make_pair(child->level(), child->index()),
1723  std::move(particle));
1724 #else
1725  position_hints[child_index] = particles.insert(
1726  position_hints[child_index],
1727  std::make_pair(std::make_pair(child->level(),
1728  child->index()),
1729  std::move(particle)));
1730 #endif
1731  // Move the hint position forward by one, i.e.,
1732  // for the next particle. The 'hint' position will
1733  // thus be right after the one just inserted.
1734  ++position_hints[child_index];
1735  break;
1736  }
1737  }
1738  catch (typename Mapping<dim>::ExcTransformationFailed &)
1739  {}
1740  }
1741  }
1742  }
1743  break;
1744 
1745  default:
1746  Assert(false, ExcInternalError());
1747  break;
1748  }
1749  }
1750 } // namespace Particles
1751 
1752 #include "particle_handler.inst"
1753 
1754 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:4371
static const unsigned int invalid_unsigned_int
Definition: types.h:187
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1571
void set_property_pool(PropertyPool &property_pool)
Definition: particle.cc:286
void add_index(const size_type index)
Definition: index_set.h:1654
types::particle_index n_locally_owned_particles() const
void get_particle_positions(VectorType &output_vector, const bool add_to_output_vector=false)
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)
IndexSet locally_relevant_ids() const
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:2342
void set_properties(const ArrayView< const double > &new_properties)
Definition: particle.cc:295
types::particle_index next_free_particle_index
SmartPointer< const Triangulation< dim, spacedim >, ParticleHandler< dim, spacedim > > triangulation
return_type distributed_compute_point_locations(const GridTools::Cache< dim, spacedim > &cache, const std::vector< Point< spacedim >> &local_points, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bboxes)
Definition: grid_tools.cc:5256
const ArrayView< double > get_properties()
Definition: particle.cc:338
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
numbers::NumberTraits< Number >::real_type norm() const
Definition: tensor.h:1443
unsigned int global_max_particles_per_cell
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:78
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:110
void initialize(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const unsigned int n_properties=0)
unsigned int particle_index
Definition: particle.h:66
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
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_properties_per_particle() const
types::particle_index global_number_of_particles
#define Assert(cond, exc)
Definition: exceptions.h:1411
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:302
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:4401
particle_iterator begin() const
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:112
void set_particle_positions(const VectorType &input_vector, const bool displace_particles=true)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:100
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:1787
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:360
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1699
std::multimap< internal::LevelInd, Particle< dim, spacedim > > particles
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:4509
SmartPointer< const Mapping< dim, spacedim >, ParticleHandler< dim, spacedim > > mapping
std::vector< T > all_gather(const MPI_Comm &comm, const T &object_to_send)
std::map< unsigned int, IndexSet > insert_global_particles(const std::vector< Point< spacedim >> &positions, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, const std::vector< std::vector< double >> &properties={})
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)
std::map< unsigned int, T > some_to_some(const MPI_Comm &comm, const std::map< unsigned int, T > &objects_to_send)
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()