Reference documentation for deal.II version GIT f9a1fd2718 2022-08-19 02:45:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
particle_handler.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2022 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 
18 
20 
22 
23 #include <memory>
24 #include <utility>
25 
27 
28 namespace Particles
29 {
30  namespace
31  {
32  template <int dim, int spacedim>
33  std::vector<char>
34  pack_particles(std::vector<ParticleIterator<dim, spacedim>> &particles)
35  {
36  std::vector<char> buffer;
37 
38  if (particles.size() == 0)
39  return buffer;
40 
41  buffer.resize(particles.size() *
42  particles.front()->serialized_size_in_bytes());
43  void *current_data = buffer.data();
44 
45  for (const auto &particle : particles)
46  {
47  current_data = particle->write_particle_data_to_memory(current_data);
48  }
49 
50  return buffer;
51  }
52  } // namespace
53 
54 
55 
56  template <int dim, int spacedim>
58  : triangulation()
59  , mapping()
60  , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(0))
61  , global_number_of_particles(0)
62  , number_of_locally_owned_particles(0)
63  , global_max_particles_per_cell(0)
64  , next_free_particle_index(0)
65  , size_callback()
66  , store_callback()
67  , load_callback()
68  , handle(numbers::invalid_unsigned_int)
69  , tria_listeners()
70  {
72  }
73 
74 
75 
76  template <int dim, int spacedim>
79  const Mapping<dim, spacedim> & mapping,
80  const unsigned int n_properties)
81  : triangulation(&triangulation, typeid(*this).name())
82  , mapping(&mapping, typeid(*this).name())
83  , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
84  , cells_to_particle_cache(triangulation.n_active_cells(), particles.end())
85  , global_number_of_particles(0)
86  , number_of_locally_owned_particles(0)
87  , global_max_particles_per_cell(0)
88  , next_free_particle_index(0)
89  , size_callback()
90  , store_callback()
91  , load_callback()
92  , handle(numbers::invalid_unsigned_int)
93  , triangulation_cache(
94  std::make_unique<GridTools::Cache<dim, spacedim>>(triangulation,
95  mapping))
96  , tria_listeners()
97  {
100  }
101 
102 
103 
104  template <int dim, int spacedim>
106  {
107  clear_particles();
108 
109  for (const auto &connection : tria_listeners)
110  connection.disconnect();
111  }
112 
113 
114 
115  template <int dim, int spacedim>
116  void
118  const Triangulation<dim, spacedim> &new_triangulation,
119  const Mapping<dim, spacedim> & new_mapping,
120  const unsigned int n_properties)
121  {
122  clear();
123 
124  triangulation = &new_triangulation;
125  mapping = &new_mapping;
126 
127  reset_particle_container(particles);
128 
129  // Create the memory pool that will store all particle properties
130  property_pool = std::make_unique<PropertyPool<dim, spacedim>>(n_properties);
131 
132  // Create the grid cache to cache the information about the triangulation
133  // that is used to locate the particles into subdomains and cells
134  triangulation_cache =
135  std::make_unique<GridTools::Cache<dim, spacedim>>(new_triangulation,
136  new_mapping);
137 
138  cells_to_particle_cache.resize(triangulation->n_active_cells(),
139  particles.end());
140 
141  connect_to_triangulation_signals();
142  }
143 
144 
145 
146  template <int dim, int spacedim>
147  void
149  const ParticleHandler<dim, spacedim> &particle_handler)
150  {
151  const unsigned int n_properties =
152  particle_handler.property_pool->n_properties_per_slot();
153  initialize(*particle_handler.triangulation,
154  *particle_handler.mapping,
155  n_properties);
156 
157  property_pool = std::make_unique<PropertyPool<dim, spacedim>>(
158  *(particle_handler.property_pool));
159 
160  // copy static members
161  global_number_of_particles = particle_handler.global_number_of_particles;
162  number_of_locally_owned_particles =
163  particle_handler.number_of_locally_owned_particles;
164 
165  global_max_particles_per_cell =
166  particle_handler.global_max_particles_per_cell;
167  next_free_particle_index = particle_handler.next_free_particle_index;
168 
169  // Manually copy over the particles because we do not want to touch the
170  // anchor iterators set by initialize()
171  particles.insert(particle_container_owned_end(),
172  particle_handler.particle_container_owned_begin(),
173  particle_handler.particle_container_owned_end());
174  particles.insert(particle_container_ghost_end(),
175  particle_handler.particle_container_ghost_begin(),
176  particle_handler.particle_container_ghost_end());
177 
178  for (auto it = particles.begin(); it != particles.end(); ++it)
179  if (!it->particles.empty())
180  cells_to_particle_cache[it->cell->active_cell_index()] = it;
181 
182  ghost_particles_cache.ghost_particles_by_domain =
183  particle_handler.ghost_particles_cache.ghost_particles_by_domain;
184  handle = particle_handler.handle;
185  }
186 
187 
188 
189  template <int dim, int spacedim>
190  void
192  {
193  clear_particles();
194  global_number_of_particles = 0;
195  number_of_locally_owned_particles = 0;
196  next_free_particle_index = 0;
197  global_max_particles_per_cell = 0;
198  }
199 
200 
201 
202  template <int dim, int spacedim>
203  void
205  {
206  for (auto &particles_in_cell : particles)
207  for (auto &particle : particles_in_cell.particles)
209  property_pool->deregister_particle(particle);
210 
211  cells_to_particle_cache.clear();
212  reset_particle_container(particles);
213  if (triangulation != nullptr)
214  cells_to_particle_cache.resize(triangulation->n_active_cells(),
215  particles.end());
216 
217  // the particle properties have already been deleted by their destructor,
218  // but the memory is still allocated. Return the memory as well.
219  property_pool->clear();
220  }
221 
222 
223 
224  template <int dim, int spacedim>
225  void
226  ParticleHandler<dim, spacedim>::reserve(const std::size_t n_particles)
227  {
228  property_pool->reserve(n_particles);
229  }
230 
231 
232 
233  template <int dim, int spacedim>
234  void
236  particle_container &given_particles)
237  {
238  // Make sure to set a valid past-the-end iterator also in case we have no
239  // triangulation
241  past_the_end_iterator =
242  triangulation != nullptr ?
243  triangulation->end() :
244  typename Triangulation<dim, spacedim>::cell_iterator(nullptr, -1, -1);
245 
246  given_particles.clear();
247  for (unsigned int i = 0; i < 3; ++i)
248  given_particles.emplace_back(
249  std::vector<typename PropertyPool<dim, spacedim>::Handle>(),
250  past_the_end_iterator);
251 
252  // Set the end of owned particles to the middle of the three elements
253  const_cast<typename particle_container::iterator &>(owned_particles_end) =
254  ++given_particles.begin();
255  }
256 
257 
258 
259  template <int dim, int spacedim>
260  void
262  {
263  // first sort the owned particles by the active cell index
264  bool sort_is_necessary = false;
265  {
266  auto previous = particle_container_owned_begin();
267  for (auto next = previous; next != particle_container_owned_end(); ++next)
268  {
269  if (previous->cell.state() == IteratorState::valid &&
270  next->cell.state() == IteratorState::valid &&
271  previous->cell > next->cell)
272  {
273  sort_is_necessary = true;
274  break;
275  }
276  previous = next;
277  }
278  }
279  if (sort_is_necessary)
280  {
281  // we could have tried to call std::list::sort with a custom
282  // comparator to get things sorted, but things are complicated by the
283  // three anchor entries that we do not want to move, and we would
284  // hence pay an O(N log(N)) algorithm with a large constant in front
285  // of it (on the order of 20+ instructions with many
286  // difficult-to-predict branches). Therefore, we simply copy the list
287  // into a new one (keeping alive the possibly large vectors with
288  // particles on cells) into a new container.
289  particle_container sorted_particles;
290 
291  // note that this call updates owned_particles_end, so that
292  // particle_container_owned_end() below already points to the
293  // new container
294  reset_particle_container(sorted_particles);
295 
296  // iterate over cells and insert the entries in the new order
297  for (const auto &cell : triangulation->active_cell_iterators())
298  if (!cell->is_artificial())
299  if (cells_to_particle_cache[cell->active_cell_index()] !=
300  particles.end())
301  {
302  // before we move the sorted_particles into particles
303  // particle_container_ghost_end() still points to the
304  // old particles container. Therefore this condition looks
305  // quirky.
306  typename particle_container::iterator insert_position =
307  cell->is_locally_owned() ? particle_container_owned_end() :
308  --sorted_particles.end();
309  typename particle_container::iterator new_entry =
310  sorted_particles.insert(
311  insert_position, typename particle_container::value_type());
312  new_entry->cell = cell;
313  new_entry->particles =
314  std::move(cells_to_particle_cache[cell->active_cell_index()]
315  ->particles);
316  }
317  particles = std::move(sorted_particles);
318 
319  // refresh cells_to_particle_cache
320  cells_to_particle_cache.clear();
321  cells_to_particle_cache.resize(triangulation->n_active_cells(),
322  particles.end());
323  for (auto it = particles.begin(); it != particles.end(); ++it)
324  if (!it->particles.empty())
325  cells_to_particle_cache[it->cell->active_cell_index()] = it;
326  }
327 
328  // Ensure that we did not accidentally modify the anchor entries with
329  // special purpose.
330  Assert(particles.front().cell.state() == IteratorState::past_the_end &&
331  particles.front().particles.empty() &&
332  particles.back().cell.state() == IteratorState::past_the_end &&
333  particles.back().particles.empty() &&
334  owned_particles_end->cell.state() == IteratorState::past_the_end &&
335  owned_particles_end->particles.empty(),
336  ExcInternalError());
337 
338 #ifdef DEBUG
339  // check that no cache element hits the three anchor states in the list of
340  // particles
341  for (const auto &it : cells_to_particle_cache)
342  Assert(it != particles.begin() && it != owned_particles_end &&
343  it != --(particles.end()),
344  ExcInternalError());
345 
346  // check that we only have locally owned particles in the first region of
347  // cells; note that we skip the very first anchor element
348  for (auto it = particle_container_owned_begin();
349  it != particle_container_owned_end();
350  ++it)
351  Assert(it->cell->is_locally_owned(), ExcInternalError());
352 
353  // check that the cache is consistent with the iterators
354  std::vector<typename particle_container::iterator> verify_cache(
355  triangulation->n_active_cells(), particles.end());
356  for (auto it = particles.begin(); it != particles.end(); ++it)
357  if (!it->particles.empty())
358  verify_cache[it->cell->active_cell_index()] = it;
359 
360  for (unsigned int i = 0; i < verify_cache.size(); ++i)
361  Assert(verify_cache[i] == cells_to_particle_cache[i], ExcInternalError());
362 #endif
363 
364  // now compute local result with the function above and then compute the
365  // collective results
366  number_of_locally_owned_particles = 0;
367 
368  types::particle_index result[2] = {};
369  for (const auto &particles_in_cell : particles)
370  {
371  const types::particle_index n_particles_in_cell =
372  particles_in_cell.particles.size();
373 
374  // local_max_particles_per_cell
375  result[0] = std::max(result[0], n_particles_in_cell);
376 
377  // number of locally owned particles
378  if (n_particles_in_cell > 0 &&
379  particles_in_cell.cell->is_locally_owned())
380  number_of_locally_owned_particles += n_particles_in_cell;
381 
382  // local_max_particle_index
383  for (const auto &particle : particles_in_cell.particles)
384  result[1] = std::max(result[1], property_pool->get_id(particle));
385  }
386 
387  global_number_of_particles =
388  ::Utilities::MPI::sum(number_of_locally_owned_particles,
389  triangulation->get_communicator());
390 
391  if (global_number_of_particles == 0)
392  {
393  next_free_particle_index = 0;
394  global_max_particles_per_cell = 0;
395  }
396  else
397  {
398  Utilities::MPI::max(result, triangulation->get_communicator(), result);
399 
400  next_free_particle_index = result[1] + 1;
401  global_max_particles_per_cell = result[0];
402  }
403  }
404 
405 
406 
407  template <int dim, int spacedim>
411  const
412  {
413  if (cells_to_particle_cache.size() == 0)
414  return 0;
415 
416  if (cell->is_artificial() == false)
417  {
418  return cells_to_particle_cache[cell->active_cell_index()] !=
419  particles.end() ?
420  cells_to_particle_cache[cell->active_cell_index()]
421  ->particles.size() :
422  0;
423  }
424  else
425  AssertThrow(false,
426  ExcMessage("You can't ask for the particles on an artificial "
427  "cell since we don't know what exists on these "
428  "kinds of cells."));
429 
431  }
432 
433 
434 
435  template <int dim, int spacedim>
439  const
440  {
441  return (const_cast<ParticleHandler<dim, spacedim> *>(this))
442  ->particles_in_cell(cell);
443  }
444 
445 
446 
447  template <int dim, int spacedim>
451  {
452  const unsigned int active_cell_index = cell->active_cell_index();
453 
454  if (cell->is_artificial() == false)
455  {
456  if (cells_to_particle_cache[active_cell_index] == particles.end())
457  {
459  particle_iterator(particles.begin(), *property_pool, 0),
460  particle_iterator(particles.begin(), *property_pool, 0));
461  }
462  else
463  {
464  const typename particle_container::iterator
465  particles_in_current_cell =
466  cells_to_particle_cache[active_cell_index];
467  typename particle_container::iterator particles_in_next_cell =
468  particles_in_current_cell;
469  ++particles_in_next_cell;
471  particle_iterator(particles_in_current_cell, *property_pool, 0),
472  particle_iterator(particles_in_next_cell, *property_pool, 0));
473  }
474  }
475  else
476  AssertThrow(false,
477  ExcMessage("You can't ask for the particles on an artificial "
478  "cell since we don't know what exists on these "
479  "kinds of cells."));
480 
481  return {};
482  }
483 
484 
485 
486  template <int dim, int spacedim>
487  void
490  {
491  auto &particles_on_cell = particle->particles_in_cell->particles;
492 
493  // if the particle has an invalid handle (e.g. because it has
494  // been duplicated before calling this function) do not try
495  // to deallocate its memory again
496  auto handle = particle->get_handle();
498  property_pool->deregister_particle(handle);
499 
500  // need to reduce the cached number before deleting, because the iterator
501  // may be invalid after removing the particle even if only
502  // accessing the cell
503  const auto cell = particle->get_surrounding_cell();
504  const bool owned_cell = cell->is_locally_owned();
505  if (owned_cell)
506  --number_of_locally_owned_particles;
507 
508  if (particles_on_cell.size() > 1)
509  {
510  particles_on_cell[particle->particle_index_within_cell] =
511  std::move(particles_on_cell.back());
512  particles_on_cell.resize(particles_on_cell.size() - 1);
513  }
514  else
515  {
516  particles.erase(particle->particles_in_cell);
517  cells_to_particle_cache[cell->active_cell_index()] = particles.end();
518  }
519  }
520 
521 
522 
523  template <int dim, int spacedim>
524  void
527  &particles_to_remove)
528  {
529  // We need to remove particles backwards on each cell to keep the particle
530  // iterators alive as we keep removing particles on the same cell. To
531  // ensure that this is safe, we either check if we already have sorted
532  // iterators or if we need to manually sort
533  const auto check_greater = [](const particle_iterator &a,
534  const particle_iterator &b) {
535  return a->particles_in_cell->cell > b->particles_in_cell->cell ||
536  (a->particles_in_cell->cell == b->particles_in_cell->cell &&
537  a->particle_index_within_cell > b->particle_index_within_cell);
538  };
539 
540  bool particles_are_sorted = true;
541  auto previous = particles_to_remove.begin();
542  for (auto next = previous; next != particles_to_remove.end(); ++next)
543  {
544  if (check_greater(*previous, *next))
545  {
546  particles_are_sorted = false;
547  break;
548  }
549  previous = next;
550  }
551  if (particles_are_sorted)
552  {
553  // pass along backwards in array
554  for (auto it = particles_to_remove.rbegin();
555  it != particles_to_remove.rend();
556  ++it)
557  remove_particle(*it);
558  }
559  else
560  {
561  std::vector<ParticleHandler<dim, spacedim>::particle_iterator>
562  sorted_particles(particles_to_remove);
563  std::sort(sorted_particles.begin(),
564  sorted_particles.end(),
565  check_greater);
566 
567  for (const auto &particle : sorted_particles)
568  remove_particle(particle);
569  }
570 
571  update_cached_numbers();
572  }
573 
574 
575 
576  template <int dim, int spacedim>
579  const Particle<dim, spacedim> & particle,
581  {
582  return insert_particle(particle.get_location(),
583  particle.get_reference_location(),
584  particle.get_id(),
585  cell,
586  particle.get_properties());
587  }
588 
589 
590 
591  template <int dim, int spacedim>
594  const typename PropertyPool<dim, spacedim>::Handle handle,
595  const typename Triangulation<dim, spacedim>::cell_iterator &cell)
596  {
597  const unsigned int active_cell_index = cell->active_cell_index();
598  typename particle_container::iterator &cache =
599  cells_to_particle_cache[active_cell_index];
600  if (cache == particles.end())
601  {
602  const typename particle_container::iterator insert_position =
603  cell->is_locally_owned() ? particle_container_owned_end() :
604  particle_container_ghost_end();
605  cache = particles.emplace(
606  insert_position,
607  std::vector<typename PropertyPool<dim, spacedim>::Handle>{handle},
608  cell);
609  }
610  else
611  {
612  cache->particles.push_back(handle);
613  Assert(cache->cell == cell, ExcInternalError());
614  }
615  return particle_iterator(cache,
616  *property_pool,
617  cache->particles.size() - 1);
618  }
619 
620 
621 
622  template <int dim, int spacedim>
625  const void *& data,
627  {
628  Assert(triangulation != nullptr, ExcInternalError());
629  Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(),
630  ExcInternalError());
631  Assert(cell->is_locally_owned(),
632  ExcMessage("You tried to insert particles into a cell that is not "
633  "locally owned. This is not supported."));
634 
635  particle_iterator particle_it =
636  insert_particle(property_pool->register_particle(), cell);
637 
638  data = particle_it->read_particle_data_from_memory(data);
639 
640  ++number_of_locally_owned_particles;
641 
642  return particle_it;
643  }
644 
645 
646 
647  template <int dim, int spacedim>
650  const Point<spacedim> & position,
651  const Point<dim> & reference_position,
654  const ArrayView<const double> &properties)
655  {
656  Assert(triangulation != nullptr, ExcInternalError());
657  Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(),
658  ExcInternalError());
659  Assert(cell.state() == IteratorState::valid, ExcInternalError());
660  Assert(cell->is_locally_owned(),
661  ExcMessage("You tried to insert particles into a cell that is not "
662  "locally owned. This is not supported."));
663 
664  particle_iterator particle_it =
665  insert_particle(property_pool->register_particle(), cell);
666 
667  particle_it->set_location(position);
668  particle_it->set_reference_location(reference_position);
669  particle_it->set_id(particle_index);
670 
671  if (properties.size() != 0)
672  particle_it->set_properties(properties);
673 
674  ++number_of_locally_owned_particles;
675 
676  return particle_it;
677  }
678 
679 
680 
681  template <int dim, int spacedim>
682  void
684  const std::multimap<
686  Particle<dim, spacedim>> &new_particles)
687  {
688  reserve(n_locally_owned_particles() + new_particles.size());
689  for (const auto &cell_and_particle : new_particles)
690  insert_particle(cell_and_particle.second, cell_and_particle.first);
691 
692  update_cached_numbers();
693  }
694 
695 
696 
697  template <int dim, int spacedim>
698  void
700  const std::vector<Point<spacedim>> &positions)
701  {
702  Assert(triangulation != nullptr, ExcInternalError());
703 
704  update_cached_numbers();
705  reserve(n_locally_owned_particles() + positions.size());
706 
707  // Determine the starting particle index of this process, which
708  // is the highest currently existing particle index plus the sum
709  // of the number of newly generated particles of all
710  // processes with a lower rank if in a parallel computation.
711  const types::particle_index local_next_particle_index =
712  get_next_free_particle_index();
713  types::particle_index local_start_index = 0;
714 
715 #ifdef DEAL_II_WITH_MPI
716  if (const auto parallel_triangulation =
717  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
718  &*triangulation))
719  {
720  types::particle_index particles_to_add_locally = positions.size();
721  const int ierr = MPI_Scan(&particles_to_add_locally,
722  &local_start_index,
723  1,
725  MPI_SUM,
726  parallel_triangulation->get_communicator());
727  AssertThrowMPI(ierr);
728  local_start_index -= particles_to_add_locally;
729  }
730 #endif
731 
732  local_start_index += local_next_particle_index;
733 
734  auto point_locations =
735  GridTools::compute_point_locations_try_all(*triangulation_cache,
736  positions);
737 
738  auto &cells = std::get<0>(point_locations);
739  auto &local_positions = std::get<1>(point_locations);
740  auto &index_map = std::get<2>(point_locations);
741  auto &missing_points = std::get<3>(point_locations);
742  // If a point was not found, throwing an error, as the old
743  // implementation of compute_point_locations would have done
744  AssertThrow(missing_points.size() == 0,
746 
747  (void)missing_points;
748 
749  for (unsigned int i = 0; i < cells.size(); ++i)
750  for (unsigned int p = 0; p < local_positions[i].size(); ++p)
751  insert_particle(positions[index_map[i][p]],
752  local_positions[i][p],
753  local_start_index + index_map[i][p],
754  cells[i]);
755 
756  update_cached_numbers();
757  }
758 
759 
760 
761  template <int dim, int spacedim>
762  std::map<unsigned int, IndexSet>
764  const std::vector<Point<spacedim>> &positions,
765  const std::vector<std::vector<BoundingBox<spacedim>>>
766  & global_bounding_boxes,
767  const std::vector<std::vector<double>> & properties,
768  const std::vector<types::particle_index> &ids)
769  {
770  if (!properties.empty())
771  {
772  AssertDimension(properties.size(), positions.size());
773 #ifdef DEBUG
774  for (const auto &p : properties)
775  AssertDimension(p.size(), n_properties_per_particle());
776 #endif
777  }
778 
779  if (!ids.empty())
780  AssertDimension(ids.size(), positions.size());
781 
782  const auto comm = triangulation->get_communicator();
783 
785 
786  // Compute the global number of properties
787  const auto n_global_properties =
788  Utilities::MPI::sum(properties.size(), comm);
789 
790  // Gather the number of points per processor
791  const auto n_particles_per_proc =
792  Utilities::MPI::all_gather(comm, positions.size());
793 
794  // Calculate all starting points locally
795  std::vector<unsigned int> particle_start_indices(n_mpi_processes);
796 
797  unsigned int particle_start_index = get_next_free_particle_index();
798  for (unsigned int process = 0; process < particle_start_indices.size();
799  ++process)
800  {
801  particle_start_indices[process] = particle_start_index;
802  particle_start_index += n_particles_per_proc[process];
803  }
804 
805  // Get all local information
806  const auto cells_positions_and_index_maps =
808  positions,
809  global_bounding_boxes);
810 
811  // Unpack the information into several vectors:
812  // All cells that contain at least one particle
813  const auto &local_cells_containing_particles =
814  std::get<0>(cells_positions_and_index_maps);
815 
816  // The reference position of every particle in the local part of the
817  // triangulation.
818  const auto &local_reference_positions =
819  std::get<1>(cells_positions_and_index_maps);
820  // The original index in the positions vector for each particle in the
821  // local part of the triangulation
822  const auto &original_indices_of_local_particles =
823  std::get<2>(cells_positions_and_index_maps);
824  // The real spatial position of every particle in the local part of the
825  // triangulation.
826  const auto &local_positions = std::get<3>(cells_positions_and_index_maps);
827  // The MPI process that inserted each particle
828  const auto &calling_process_indices =
829  std::get<4>(cells_positions_and_index_maps);
830 
831  // Create the map of cpu to indices, indicating who sent us what particle
832  std::map<unsigned int, std::vector<unsigned int>>
833  original_process_to_local_particle_indices_tmp;
834  for (unsigned int i_cell = 0;
835  i_cell < local_cells_containing_particles.size();
836  ++i_cell)
837  {
838  for (unsigned int i_particle = 0;
839  i_particle < local_positions[i_cell].size();
840  ++i_particle)
841  {
842  const unsigned int local_id_on_calling_process =
843  original_indices_of_local_particles[i_cell][i_particle];
844  const unsigned int calling_process =
845  calling_process_indices[i_cell][i_particle];
846 
847  original_process_to_local_particle_indices_tmp[calling_process]
848  .push_back(local_id_on_calling_process);
849  }
850  }
851  std::map<unsigned int, IndexSet> original_process_to_local_particle_indices;
852  for (auto &process_and_particle_indices :
853  original_process_to_local_particle_indices_tmp)
854  {
855  const unsigned int calling_process = process_and_particle_indices.first;
856  original_process_to_local_particle_indices.insert(
857  {calling_process, IndexSet(n_particles_per_proc[calling_process])});
858  std::sort(process_and_particle_indices.second.begin(),
859  process_and_particle_indices.second.end());
860  original_process_to_local_particle_indices[calling_process].add_indices(
861  process_and_particle_indices.second.begin(),
862  process_and_particle_indices.second.end());
863  original_process_to_local_particle_indices[calling_process].compress();
864  }
865 
866  // A map from mpi process to properties, ordered as in the IndexSet.
867  // Notice that this ordering may be different from the ordering in the
868  // vectors above, since no local ordering is guaranteed by the
869  // distribute_compute_point_locations() call.
870  // This is only filled if n_global_properties is > 0
871  std::map<unsigned int, std::vector<std::vector<double>>>
872  locally_owned_properties_from_other_processes;
873 
874  // A map from mpi process to ids, ordered as in the IndexSet.
875  // Notice that this ordering may be different from the ordering in the
876  // vectors above, since no local ordering is guaranteed by the
877  // distribute_compute_point_locations() call.
878  // This is only filled if ids.size() is > 0
879  std::map<unsigned int, std::vector<types::particle_index>>
880  locally_owned_ids_from_other_processes;
881 
882  if (n_global_properties > 0 || !ids.empty())
883  {
884  // Gather whom I sent my own particles to, to decide whom to send
885  // the particle properties or the ids
886  auto send_to_cpu = Utilities::MPI::some_to_some(
887  comm, original_process_to_local_particle_indices);
888 
889  // Prepare the vector of properties to send
890  if (n_global_properties > 0)
891  {
892  std::map<unsigned int, std::vector<std::vector<double>>>
893  non_locally_owned_properties;
894 
895  for (const auto &it : send_to_cpu)
896  {
897  std::vector<std::vector<double>> properties_to_send(
898  it.second.n_elements(),
899  std::vector<double>(n_properties_per_particle()));
900  unsigned int index = 0;
901  for (const auto el : it.second)
902  properties_to_send[index++] = properties[el];
903  non_locally_owned_properties.insert(
904  {it.first, properties_to_send});
905  }
906 
907  // Send the non locally owned properties to each mpi process
908  // that needs them
909  locally_owned_properties_from_other_processes =
910  Utilities::MPI::some_to_some(comm, non_locally_owned_properties);
911 
913  locally_owned_properties_from_other_processes.size(),
914  original_process_to_local_particle_indices.size());
915  }
916 
917  if (!ids.empty())
918  {
919  std::map<unsigned int, std::vector<types::particle_index>>
920  non_locally_owned_ids;
921  for (const auto &it : send_to_cpu)
922  {
923  std::vector<types::particle_index> ids_to_send(
924  it.second.n_elements());
925  unsigned int index = 0;
926  for (const auto el : it.second)
927  ids_to_send[index++] = ids[el];
928  non_locally_owned_ids.insert({it.first, ids_to_send});
929  }
930 
931  // Send the non locally owned ids to each mpi process
932  // that needs them
933  locally_owned_ids_from_other_processes =
934  Utilities::MPI::some_to_some(comm, non_locally_owned_ids);
935 
936  AssertDimension(locally_owned_ids_from_other_processes.size(),
937  original_process_to_local_particle_indices.size());
938  }
939  }
940 
941  // Now fill up the actual particles
942  for (unsigned int i_cell = 0;
943  i_cell < local_cells_containing_particles.size();
944  ++i_cell)
945  {
946  for (unsigned int i_particle = 0;
947  i_particle < local_positions[i_cell].size();
948  ++i_particle)
949  {
950  const unsigned int local_id_on_calling_process =
951  original_indices_of_local_particles[i_cell][i_particle];
952 
953  const unsigned int calling_process =
954  calling_process_indices[i_cell][i_particle];
955 
956  const unsigned int index_within_set =
957  original_process_to_local_particle_indices[calling_process]
958  .index_within_set(local_id_on_calling_process);
959 
960  const unsigned int particle_id =
961  ids.empty() ?
962  local_id_on_calling_process +
963  particle_start_indices[calling_process] :
964  locally_owned_ids_from_other_processes[calling_process]
965  [index_within_set];
966 
967  auto particle_it =
968  insert_particle(local_positions[i_cell][i_particle],
969  local_reference_positions[i_cell][i_particle],
970  particle_id,
971  local_cells_containing_particles[i_cell]);
972 
973  if (n_global_properties > 0)
974  {
975  particle_it->set_properties(
976  locally_owned_properties_from_other_processes
977  [calling_process][index_within_set]);
978  }
979  }
980  }
981 
982  update_cached_numbers();
983 
984  return original_process_to_local_particle_indices;
985  }
986 
987 
988 
989  template <int dim, int spacedim>
990  std::map<unsigned int, IndexSet>
992  const std::vector<Particle<dim, spacedim>> &particles,
993  const std::vector<std::vector<BoundingBox<spacedim>>>
994  &global_bounding_boxes)
995  {
996  // Store the positions in a vector of points, the ids in a vector of ids,
997  // and the properties, if any, in a vector of vector of properties.
998  std::vector<Point<spacedim>> positions;
999  std::vector<std::vector<double>> properties;
1000  std::vector<types::particle_index> ids;
1001  positions.resize(particles.size());
1002  ids.resize(particles.size());
1003  if (n_properties_per_particle() > 0)
1004  properties.resize(particles.size(),
1005  std::vector<double>(n_properties_per_particle()));
1006 
1007  unsigned int i = 0;
1008  for (const auto &p : particles)
1009  {
1010  positions[i] = p.get_location();
1011  ids[i] = p.get_id();
1012  if (p.has_properties())
1013  properties[i] = {p.get_properties().begin(),
1014  p.get_properties().end()};
1015  ++i;
1016  }
1017 
1018  return insert_global_particles(positions,
1019  global_bounding_boxes,
1020  properties,
1021  ids);
1022  }
1023 
1024 
1025 
1026  template <int dim, int spacedim>
1029  {
1030  return global_number_of_particles;
1031  }
1032 
1033 
1034 
1035  template <int dim, int spacedim>
1038  {
1039  return global_max_particles_per_cell;
1040  }
1041 
1042 
1043 
1044  template <int dim, int spacedim>
1047  {
1048  return number_of_locally_owned_particles;
1049  }
1050 
1051 
1052 
1053  template <int dim, int spacedim>
1054  unsigned int
1056  {
1057  return property_pool->n_properties_per_slot();
1058  }
1059 
1060 
1061 
1062  template <int dim, int spacedim>
1065  {
1066  return next_free_particle_index;
1067  }
1068 
1069 
1070 
1071  template <int dim, int spacedim>
1072  IndexSet
1074  {
1075  IndexSet set(get_next_free_particle_index());
1076  std::vector<types::particle_index> indices;
1077  indices.reserve(n_locally_owned_particles());
1078  for (const auto &p : *this)
1079  indices.push_back(p.get_id());
1080  set.add_indices(indices.begin(), indices.end());
1081  set.compress();
1082  return set;
1083  }
1084 
1085 
1086 
1087  template <int dim, int spacedim>
1090  {
1091  return property_pool->n_slots();
1092  }
1093 
1094 
1095 
1096  template <int dim, int spacedim>
1097  void
1099  std::vector<Point<spacedim>> &positions,
1100  const bool add_to_output_vector)
1101  {
1102  // There should be one point per particle to gather
1103  AssertDimension(positions.size(), n_locally_owned_particles());
1104 
1105  unsigned int i = 0;
1106  for (auto it = begin(); it != end(); ++it, ++i)
1107  {
1108  if (add_to_output_vector)
1109  positions[i] = positions[i] + it->get_location();
1110  else
1111  positions[i] = it->get_location();
1112  }
1113  }
1114 
1115 
1116 
1117  template <int dim, int spacedim>
1118  void
1120  const std::vector<Point<spacedim>> &new_positions,
1121  const bool displace_particles)
1122  {
1123  // There should be one point per particle to fix the new position
1124  AssertDimension(new_positions.size(), n_locally_owned_particles());
1125 
1126  unsigned int i = 0;
1127  for (auto it = begin(); it != end(); ++it, ++i)
1128  {
1129  if (displace_particles)
1130  it->set_location(it->get_location() + new_positions[i]);
1131  else
1132  it->set_location(new_positions[i]);
1133  }
1134  sort_particles_into_subdomains_and_cells();
1135  }
1136 
1137 
1138 
1139  template <int dim, int spacedim>
1140  void
1142  const Function<spacedim> &function,
1143  const bool displace_particles)
1144  {
1145  // The function should have sufficient components to displace the
1146  // particles
1147  AssertDimension(function.n_components, spacedim);
1148 
1149  Vector<double> new_position(spacedim);
1150  for (auto &particle : *this)
1151  {
1152  Point<spacedim> particle_location = particle.get_location();
1153  function.vector_value(particle_location, new_position);
1154  if (displace_particles)
1155  for (unsigned int d = 0; d < spacedim; ++d)
1156  particle_location[d] += new_position[d];
1157  else
1158  for (unsigned int d = 0; d < spacedim; ++d)
1159  particle_location[d] = new_position[d];
1160  particle.set_location(particle_location);
1161  }
1162  sort_particles_into_subdomains_and_cells();
1163  }
1164 
1165 
1166 
1167  template <int dim, int spacedim>
1170  {
1171  return *property_pool;
1172  }
1173 
1174 
1175 
1176  namespace
1177  {
1187  template <int dim>
1188  bool
1189  compare_particle_association(
1190  const unsigned int a,
1191  const unsigned int b,
1192  const Tensor<1, dim> & particle_direction,
1193  const std::vector<Tensor<1, dim>> &center_directions)
1194  {
1195  const double scalar_product_a = center_directions[a] * particle_direction;
1196  const double scalar_product_b = center_directions[b] * particle_direction;
1197 
1198  // The function is supposed to return if a is before b. We are looking
1199  // for the alignment of particle direction and center direction,
1200  // therefore return if the scalar product of a is larger.
1201  return (scalar_product_a > scalar_product_b);
1202  }
1203  } // namespace
1204 
1205 
1206 
1207  template <int dim, int spacedim>
1208  void
1210  {
1211  Assert(triangulation != nullptr, ExcInternalError());
1212  Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(),
1213  ExcInternalError());
1214 
1215  // TODO: The current algorithm only works for particles that are in
1216  // the local domain or in ghost cells, because it only knows the
1217  // subdomain_id of ghost cells, but not of artificial cells. This
1218  // can be extended using the distributed version of compute point
1219  // locations.
1220  // TODO: Extend this function to allow keeping particles on other
1221  // processes around (with an invalid cell).
1222 
1223  std::vector<particle_iterator> particles_out_of_cell;
1224 
1225  // Reserve some space for particles that need sorting to avoid frequent
1226  // re-allocation. Guess 25% of particles need sorting. Balance memory
1227  // overhead and performance.
1228  particles_out_of_cell.reserve(n_locally_owned_particles() / 4);
1229 
1230  // Now update the reference locations of the moved particles
1231  std::vector<Point<spacedim>> real_locations;
1232  std::vector<Point<dim>> reference_locations;
1233  real_locations.reserve(global_max_particles_per_cell);
1234  reference_locations.reserve(global_max_particles_per_cell);
1235 
1236  for (const auto &cell : triangulation->active_cell_iterators())
1237  {
1238  // Particles can be inserted into arbitrary cells, e.g. if their cell is
1239  // not known. However, for artificial cells we can not evaluate
1240  // the reference position of particles. Do not sort particles that are
1241  // not locally owned, because they will be sorted by the process that
1242  // owns them.
1243  if (cell->is_locally_owned() == false)
1244  {
1245  continue;
1246  }
1247 
1248  const unsigned int n_pic = n_particles_in_cell(cell);
1249  auto pic = particles_in_cell(cell);
1250 
1251  real_locations.clear();
1252  for (const auto &particle : pic)
1253  real_locations.push_back(particle.get_location());
1254 
1255  reference_locations.resize(n_pic);
1256  mapping->transform_points_real_to_unit_cell(cell,
1257  real_locations,
1258  reference_locations);
1259 
1260  auto particle = pic.begin();
1261  for (const auto &p_unit : reference_locations)
1262  {
1263  if (p_unit[0] == std::numeric_limits<double>::infinity() ||
1265  particles_out_of_cell.push_back(particle);
1266  else
1267  particle->set_reference_location(p_unit);
1268 
1269  ++particle;
1270  }
1271  }
1272 
1273  // There are three reasons why a particle is not in its old cell:
1274  // It moved to another cell, to another subdomain or it left the mesh.
1275  // Particles that moved to another cell are updated and moved inside the
1276  // particles vector, particles that moved to another domain are
1277  // collected in the moved_particles_domain vector. Particles that left
1278  // the mesh completely are ignored and removed.
1279  std::map<types::subdomain_id, std::vector<particle_iterator>>
1280  moved_particles;
1281  std::map<
1283  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1284  moved_cells;
1285 
1286  // We do not know exactly how many particles are lost, exchanged between
1287  // domains, or remain on this process. Therefore we pre-allocate
1288  // approximate sizes for these vectors. If more space is needed an
1289  // automatic and relatively fast (compared to other parts of this
1290  // algorithm) re-allocation will happen.
1291  using vector_size = typename std::vector<particle_iterator>::size_type;
1292 
1293  std::set<types::subdomain_id> ghost_owners;
1294  if (const auto parallel_triangulation =
1295  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1296  &*triangulation))
1297  ghost_owners = parallel_triangulation->ghost_owners();
1298 
1299  // Reserve some space for particles that need communication to avoid
1300  // frequent re-allocation. Guess 25% of particles out of their old cell need
1301  // communication. Balance memory overhead and performance.
1302  for (const auto &ghost_owner : ghost_owners)
1303  moved_particles[ghost_owner].reserve(particles_out_of_cell.size() / 4);
1304  for (const auto &ghost_owner : ghost_owners)
1305  moved_cells[ghost_owner].reserve(particles_out_of_cell.size() / 4);
1306 
1307  {
1308  // Create a map from vertices to adjacent cells using grid cache
1309  const std::vector<
1310  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1311  &vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
1312 
1313  // Create a corresponding map of vectors from vertex to cell center
1314  // using grid cache
1315  const std::vector<std::vector<Tensor<1, spacedim>>>
1316  &vertex_to_cell_centers =
1317  triangulation_cache->get_vertex_to_cell_centers_directions();
1318 
1319  std::vector<unsigned int> neighbor_permutation;
1320 
1321  // Reuse these vectors below, but only with a single element.
1322  // Avoid resizing for every particle.
1323  Point<dim> invalid_reference_point;
1324  Point<spacedim> invalid_point;
1325  invalid_reference_point[0] = std::numeric_limits<double>::infinity();
1326  invalid_point[0] = std::numeric_limits<double>::infinity();
1327  reference_locations.resize(1, invalid_reference_point);
1328  real_locations.resize(1, invalid_point);
1329 
1330  // Find the cells that the particles moved to.
1331  for (auto &out_particle : particles_out_of_cell)
1332  {
1333  // make a copy of the current cell, since we will modify the
1334  // variable current_cell below, but we need the original in
1335  // the case the particle is not found
1336  auto current_cell = out_particle->get_surrounding_cell();
1337 
1338  real_locations[0] = out_particle->get_location();
1339 
1340  // Record if the new cell was found
1341  bool found_cell = false;
1342 
1343  // Check if the particle is in one of the old cell's neighbors
1344  // that are adjacent to the closest vertex
1345  const unsigned int closest_vertex =
1346  GridTools::find_closest_vertex_of_cell<dim, spacedim>(
1347  current_cell, out_particle->get_location(), *mapping);
1348  Tensor<1, spacedim> vertex_to_particle =
1349  out_particle->get_location() - current_cell->vertex(closest_vertex);
1350  vertex_to_particle /= vertex_to_particle.norm();
1351 
1352  const unsigned int closest_vertex_index =
1353  current_cell->vertex_index(closest_vertex);
1354  const unsigned int n_neighbor_cells =
1355  vertex_to_cells[closest_vertex_index].size();
1356 
1357  neighbor_permutation.resize(n_neighbor_cells);
1358  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1359  neighbor_permutation[i] = i;
1360 
1361  const auto &cell_centers =
1362  vertex_to_cell_centers[closest_vertex_index];
1363  std::sort(neighbor_permutation.begin(),
1364  neighbor_permutation.end(),
1365  [&vertex_to_particle, &cell_centers](const unsigned int a,
1366  const unsigned int b) {
1367  return compare_particle_association(a,
1368  b,
1369  vertex_to_particle,
1370  cell_centers);
1371  });
1372 
1373  // Search all of the cells adjacent to the closest vertex of the
1374  // previous cell Most likely we will find the particle in them.
1375  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1376  {
1377  typename std::set<typename Triangulation<dim, spacedim>::
1378  active_cell_iterator>::const_iterator cell =
1379  vertex_to_cells[closest_vertex_index].begin();
1380 
1381  std::advance(cell, neighbor_permutation[i]);
1382  mapping->transform_points_real_to_unit_cell(*cell,
1383  real_locations,
1384  reference_locations);
1385 
1387  reference_locations[0]))
1388  {
1389  current_cell = *cell;
1390  found_cell = true;
1391  break;
1392  }
1393  }
1394 
1395  if (!found_cell)
1396  {
1397  // The particle is not in a neighbor of the old cell.
1398  // Look for the new cell in the whole local domain.
1399  // This case is rare.
1400  std::vector<std::pair<Point<spacedim>, unsigned int>>
1401  closest_vertex_in_domain;
1402  triangulation_cache->get_used_vertices_rtree().query(
1403  boost::geometry::index::nearest(out_particle->get_location(),
1404  1),
1405  std::back_inserter(closest_vertex_in_domain));
1406 
1407  // We should have one and only one result
1408  AssertDimension(closest_vertex_in_domain.size(), 1);
1409  const unsigned int closest_vertex_index_in_domain =
1410  closest_vertex_in_domain[0].second;
1411 
1412  // Search all of the cells adjacent to the closest vertex of the
1413  // domain. Most likely we will find the particle in them.
1414  for (const auto &cell :
1415  vertex_to_cells[closest_vertex_index_in_domain])
1416  {
1417  mapping->transform_points_real_to_unit_cell(
1418  cell, real_locations, reference_locations);
1419 
1421  reference_locations[0]))
1422  {
1423  current_cell = cell;
1424  found_cell = true;
1425  break;
1426  }
1427  }
1428  }
1429 
1430  if (!found_cell)
1431  {
1432  // We can find no cell for this particle. It has left the
1433  // domain due to an integration error or an open boundary.
1434  // Signal the loss and move on.
1435  signals.particle_lost(out_particle,
1436  out_particle->get_surrounding_cell());
1437  continue;
1438  }
1439 
1440  // If we are here, we found a cell and reference position for this
1441  // particle
1442  out_particle->set_reference_location(reference_locations[0]);
1443 
1444  // Reinsert the particle into our domain if we own its cell.
1445  // Mark it for MPI transfer otherwise
1446  if (current_cell->is_locally_owned())
1447  {
1448  typename PropertyPool<dim, spacedim>::Handle &old =
1449  out_particle->particles_in_cell
1450  ->particles[out_particle->particle_index_within_cell];
1451 
1452  // Avoid deallocating the memory of this particle
1453  const auto old_value = old;
1455 
1456  // Allocate particle with the old handle
1457  insert_particle(old_value, current_cell);
1458  }
1459  else
1460  {
1461  moved_particles[current_cell->subdomain_id()].push_back(
1462  out_particle);
1463  moved_cells[current_cell->subdomain_id()].push_back(current_cell);
1464  }
1465  }
1466  }
1467 
1468  // Exchange particles between processors if we have more than one process
1469 #ifdef DEAL_II_WITH_MPI
1470  if (const auto parallel_triangulation =
1471  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1472  &*triangulation))
1473  {
1475  parallel_triangulation->get_communicator()) > 1)
1476  send_recv_particles(moved_particles, moved_cells);
1477  }
1478 #endif
1479 
1480  // remove_particles also calls update_cached_numbers()
1481  remove_particles(particles_out_of_cell);
1482 
1483  // now make sure particle data is sorted in order of iteration
1484  std::vector<typename PropertyPool<dim, spacedim>::Handle> unsorted_handles;
1485  unsorted_handles.reserve(property_pool->n_registered_slots());
1486 
1487  typename PropertyPool<dim, spacedim>::Handle sorted_handle = 0;
1488  for (auto &particles_in_cell : particles)
1489  for (auto &particle : particles_in_cell.particles)
1490  {
1491  unsorted_handles.push_back(particle);
1492  particle = sorted_handle++;
1493  }
1494 
1495  property_pool->sort_memory_slots(unsorted_handles);
1496 
1497  } // namespace Particles
1498 
1499 
1500 
1501  template <int dim, int spacedim>
1502  void
1504  const bool enable_cache)
1505  {
1506  // Nothing to do in serial computations
1507  const auto parallel_triangulation =
1508  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1509  &*triangulation);
1510  if (parallel_triangulation != nullptr)
1511  {
1513  parallel_triangulation->get_communicator()) == 1)
1514  return;
1515  }
1516  else
1517  return;
1518 
1519 #ifndef DEAL_II_WITH_MPI
1520  (void)enable_cache;
1521 #else
1522  // Clear ghost particles and their properties
1523  for (const auto &cell : triangulation->active_cell_iterators())
1524  if (cell->is_ghost() &&
1525  cells_to_particle_cache[cell->active_cell_index()] != particles.end())
1526  {
1527  Assert(cells_to_particle_cache[cell->active_cell_index()]->cell ==
1528  cell,
1529  ExcInternalError());
1530  // Clear particle properties
1531  for (auto &ghost_particle :
1532  cells_to_particle_cache[cell->active_cell_index()]->particles)
1533  property_pool->deregister_particle(ghost_particle);
1534 
1535  // Clear particles themselves
1536  particles.erase(cells_to_particle_cache[cell->active_cell_index()]);
1537  cells_to_particle_cache[cell->active_cell_index()] = particles.end();
1538  }
1539 
1540  // Clear ghost particles cache and invalidate it
1541  ghost_particles_cache.ghost_particles_by_domain.clear();
1542  ghost_particles_cache.valid = false;
1543 
1544  const std::set<types::subdomain_id> ghost_owners =
1545  parallel_triangulation->ghost_owners();
1546  for (const auto ghost_owner : ghost_owners)
1547  ghost_particles_cache.ghost_particles_by_domain[ghost_owner].reserve(
1548  n_locally_owned_particles() / 4);
1549 
1550  const std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain =
1551  triangulation_cache->get_vertex_to_neighbor_subdomain();
1552 
1553  for (const auto &cell : triangulation->active_cell_iterators())
1554  {
1555  if (cell->is_locally_owned())
1556  {
1557  std::set<unsigned int> cell_to_neighbor_subdomain;
1558  for (const unsigned int v : cell->vertex_indices())
1559  {
1560  cell_to_neighbor_subdomain.insert(
1561  vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
1562  vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
1563  }
1564 
1565  if (cell_to_neighbor_subdomain.size() > 0)
1566  {
1567  const particle_iterator_range particle_range =
1568  particles_in_cell(cell);
1569 
1570  for (const auto domain : cell_to_neighbor_subdomain)
1571  {
1572  for (typename particle_iterator_range::iterator particle =
1573  particle_range.begin();
1574  particle != particle_range.end();
1575  ++particle)
1576  ghost_particles_cache.ghost_particles_by_domain[domain]
1577  .push_back(particle);
1578  }
1579  }
1580  }
1581  }
1582 
1583  send_recv_particles(
1584  ghost_particles_cache.ghost_particles_by_domain,
1585  std::map<
1587  std::vector<
1589  enable_cache);
1590 #endif
1591  }
1592 
1593 
1594 
1595  template <int dim, int spacedim>
1596  void
1598  {
1599  // Nothing to do in serial computations
1600  const auto parallel_triangulation =
1601  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1602  &*triangulation);
1603  if (parallel_triangulation == nullptr ||
1605  parallel_triangulation->get_communicator()) == 1)
1606  {
1607  return;
1608  }
1609 
1610 
1611 #ifdef DEAL_II_WITH_MPI
1612  // First clear the current ghost_particle information
1613  // ghost_particles.clear();
1614  Assert(ghost_particles_cache.valid,
1615  ExcMessage(
1616  "Ghost particles cannot be updated if they first have not been "
1617  "exchanged at least once with the cache enabled"));
1618 
1619 
1620  send_recv_particles_properties_and_location(
1621  ghost_particles_cache.ghost_particles_by_domain);
1622 #endif
1623  }
1624 
1625 
1626 
1627 #ifdef DEAL_II_WITH_MPI
1628  template <int dim, int spacedim>
1629  void
1631  const std::map<types::subdomain_id, std::vector<particle_iterator>>
1632  &particles_to_send,
1633  const std::map<
1636  & send_cells,
1637  const bool build_cache)
1638  {
1639  Assert(triangulation != nullptr, ExcInternalError());
1640  Assert(cells_to_particle_cache.size() == triangulation->n_active_cells(),
1641  ExcInternalError());
1642 
1643  ghost_particles_cache.valid = build_cache;
1644 
1645  const auto parallel_triangulation =
1646  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1647  &*triangulation);
1648  Assert(parallel_triangulation,
1649  ExcMessage("This function is only implemented for "
1650  "parallel::TriangulationBase objects."));
1651 
1652  // Determine the communication pattern
1653  const std::set<types::subdomain_id> ghost_owners =
1654  parallel_triangulation->ghost_owners();
1655  const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
1656  ghost_owners.end());
1657  const unsigned int n_neighbors = neighbors.size();
1658 
1659  if (send_cells.size() != 0)
1660  Assert(particles_to_send.size() == send_cells.size(), ExcInternalError());
1661 
1662  // If we do not know the subdomain this particle needs to be send to,
1663  // throw an error
1664  Assert(particles_to_send.find(numbers::artificial_subdomain_id) ==
1665  particles_to_send.end(),
1666  ExcInternalError());
1667 
1668  // TODO: Implement the shipping of particles to processes that are not
1669  // ghost owners of the local domain
1670  for (auto send_particles = particles_to_send.begin();
1671  send_particles != particles_to_send.end();
1672  ++send_particles)
1673  Assert(ghost_owners.find(send_particles->first) != ghost_owners.end(),
1674  ExcNotImplemented());
1675 
1676  std::size_t n_send_particles = 0;
1677  for (auto send_particles = particles_to_send.begin();
1678  send_particles != particles_to_send.end();
1679  ++send_particles)
1680  n_send_particles += send_particles->second.size();
1681 
1682  const unsigned int cellid_size = sizeof(CellId::binary_type);
1683 
1684  // Containers for the amount and offsets of data we will send
1685  // to other processors and the data itself.
1686  std::vector<unsigned int> n_send_data(n_neighbors, 0);
1687  std::vector<unsigned int> send_offsets(n_neighbors, 0);
1688  std::vector<char> send_data;
1689 
1690  Particle<dim, spacedim> test_particle;
1691  test_particle.set_property_pool(*property_pool);
1692 
1693  const unsigned int individual_particle_data_size =
1694  test_particle.serialized_size_in_bytes() +
1695  (size_callback ? size_callback() : 0);
1696 
1697  const unsigned int individual_total_particle_data_size =
1698  individual_particle_data_size + cellid_size;
1699 
1700  // Only serialize things if there are particles to be send.
1701  // We can not return early even if no particles
1702  // are send, because we might receive particles from other processes
1703  if (n_send_particles > 0)
1704  {
1705  // Allocate space for sending particle data
1706  send_data.resize(n_send_particles *
1707  individual_total_particle_data_size);
1708 
1709  void *data = static_cast<void *>(&send_data.front());
1710 
1711  // Serialize the data sorted by receiving process
1712  for (unsigned int i = 0; i < n_neighbors; ++i)
1713  {
1714  send_offsets[i] = reinterpret_cast<std::size_t>(data) -
1715  reinterpret_cast<std::size_t>(&send_data.front());
1716 
1717  const unsigned int n_particles_to_send =
1718  particles_to_send.at(neighbors[i]).size();
1719 
1720  Assert(static_cast<std::size_t>(n_particles_to_send) *
1721  individual_total_particle_data_size ==
1722  static_cast<std::size_t>(
1723  n_particles_to_send *
1724  individual_total_particle_data_size),
1725  ExcMessage("Overflow when trying to send particle "
1726  "data"));
1727 
1728  for (unsigned int j = 0; j < n_particles_to_send; ++j)
1729  {
1730  // If no target cells are given, use the iterator
1731  // information
1733  cell;
1734  if (send_cells.size() == 0)
1735  cell = particles_to_send.at(neighbors[i])[j]
1736  ->get_surrounding_cell();
1737  else
1738  cell = send_cells.at(neighbors[i])[j];
1739 
1740  const CellId::binary_type cellid =
1741  cell->id().template to_binary<dim>();
1742  memcpy(data, &cellid, cellid_size);
1743  data = static_cast<char *>(data) + cellid_size;
1744 
1745  data = particles_to_send.at(neighbors[i])[j]
1746  ->write_particle_data_to_memory(data);
1747  if (store_callback)
1748  data =
1749  store_callback(particles_to_send.at(neighbors[i])[j], data);
1750  }
1751  n_send_data[i] = n_particles_to_send;
1752  }
1753  }
1754 
1755  // Containers for the data we will receive from other processors
1756  std::vector<unsigned int> n_recv_data(n_neighbors);
1757  std::vector<unsigned int> recv_offsets(n_neighbors);
1758 
1759  {
1760  const int mpi_tag = Utilities::MPI::internal::Tags::
1762 
1763  std::vector<MPI_Request> n_requests(2 * n_neighbors);
1764  for (unsigned int i = 0; i < n_neighbors; ++i)
1765  {
1766  const int ierr = MPI_Irecv(&(n_recv_data[i]),
1767  1,
1768  MPI_UNSIGNED,
1769  neighbors[i],
1770  mpi_tag,
1771  parallel_triangulation->get_communicator(),
1772  &(n_requests[2 * i]));
1773  AssertThrowMPI(ierr);
1774  }
1775  for (unsigned int i = 0; i < n_neighbors; ++i)
1776  {
1777  const int ierr = MPI_Isend(&(n_send_data[i]),
1778  1,
1779  MPI_UNSIGNED,
1780  neighbors[i],
1781  mpi_tag,
1782  parallel_triangulation->get_communicator(),
1783  &(n_requests[2 * i + 1]));
1784  AssertThrowMPI(ierr);
1785  }
1786  const int ierr =
1787  MPI_Waitall(2 * n_neighbors, n_requests.data(), MPI_STATUSES_IGNORE);
1788  AssertThrowMPI(ierr);
1789  }
1790 
1791  // Determine how many particles and data we will receive
1792  unsigned int total_recv_data = 0;
1793  for (unsigned int neighbor_id = 0; neighbor_id < n_neighbors; ++neighbor_id)
1794  {
1795  recv_offsets[neighbor_id] = total_recv_data;
1796  total_recv_data +=
1797  n_recv_data[neighbor_id] * individual_total_particle_data_size;
1798  }
1799 
1800  // Set up the space for the received particle data
1801  std::vector<char> recv_data(total_recv_data);
1802 
1803  // Exchange the particle data between domains
1804  {
1805  std::vector<MPI_Request> requests(2 * n_neighbors);
1806  unsigned int send_ops = 0;
1807  unsigned int recv_ops = 0;
1808 
1809  const int mpi_tag = Utilities::MPI::internal::Tags::
1811 
1812  for (unsigned int i = 0; i < n_neighbors; ++i)
1813  if (n_recv_data[i] > 0)
1814  {
1815  const int ierr =
1816  MPI_Irecv(&(recv_data[recv_offsets[i]]),
1817  n_recv_data[i] * individual_total_particle_data_size,
1818  MPI_CHAR,
1819  neighbors[i],
1820  mpi_tag,
1821  parallel_triangulation->get_communicator(),
1822  &(requests[send_ops]));
1823  AssertThrowMPI(ierr);
1824  send_ops++;
1825  }
1826 
1827  for (unsigned int i = 0; i < n_neighbors; ++i)
1828  if (n_send_data[i] > 0)
1829  {
1830  const int ierr =
1831  MPI_Isend(&(send_data[send_offsets[i]]),
1832  n_send_data[i] * individual_total_particle_data_size,
1833  MPI_CHAR,
1834  neighbors[i],
1835  mpi_tag,
1836  parallel_triangulation->get_communicator(),
1837  &(requests[send_ops + recv_ops]));
1838  AssertThrowMPI(ierr);
1839  recv_ops++;
1840  }
1841  const int ierr =
1842  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
1843  AssertThrowMPI(ierr);
1844  }
1845 
1846  // Put the received particles into the domain if they are in the
1847  // triangulation
1848  const void *recv_data_it = static_cast<const void *>(recv_data.data());
1849 
1850  // Store the particle iterators in the cache
1851  auto &ghost_particles_iterators =
1852  ghost_particles_cache.ghost_particles_iterators;
1853 
1854  if (build_cache)
1855  {
1856  ghost_particles_iterators.clear();
1857 
1858  auto &send_pointers_particles = ghost_particles_cache.send_pointers;
1859  send_pointers_particles.assign(n_neighbors + 1, 0);
1860 
1861  for (unsigned int i = 0; i < n_neighbors; ++i)
1862  send_pointers_particles[i + 1] =
1863  send_pointers_particles[i] +
1864  n_send_data[i] * individual_particle_data_size;
1865 
1866  auto &recv_pointers_particles = ghost_particles_cache.recv_pointers;
1867  recv_pointers_particles.assign(n_neighbors + 1, 0);
1868 
1869  for (unsigned int i = 0; i < n_neighbors; ++i)
1870  recv_pointers_particles[i + 1] =
1871  recv_pointers_particles[i] +
1872  n_recv_data[i] * individual_particle_data_size;
1873 
1874  ghost_particles_cache.neighbors = neighbors;
1875 
1876  ghost_particles_cache.send_data.resize(
1877  ghost_particles_cache.send_pointers.back());
1878  ghost_particles_cache.recv_data.resize(
1879  ghost_particles_cache.recv_pointers.back());
1880  }
1881 
1882  while (reinterpret_cast<std::size_t>(recv_data_it) -
1883  reinterpret_cast<std::size_t>(recv_data.data()) <
1884  total_recv_data)
1885  {
1886  CellId::binary_type binary_cellid;
1887  memcpy(&binary_cellid, recv_data_it, cellid_size);
1888  const CellId id(binary_cellid);
1889  recv_data_it = static_cast<const char *>(recv_data_it) + cellid_size;
1890 
1892  triangulation->create_cell_iterator(id);
1893 
1894  insert_particle(property_pool->register_particle(), cell);
1895  const typename particle_container::iterator &cache =
1896  cells_to_particle_cache[cell->active_cell_index()];
1897  Assert(cache->cell == cell, ExcInternalError());
1898 
1899  particle_iterator particle_it(cache,
1900  *property_pool,
1901  cache->particles.size() - 1);
1902 
1903  recv_data_it =
1904  particle_it->read_particle_data_from_memory(recv_data_it);
1905 
1906  if (load_callback)
1907  recv_data_it = load_callback(particle_it, recv_data_it);
1908 
1909  if (build_cache) // TODO: is this safe?
1910  ghost_particles_iterators.push_back(particle_it);
1911  }
1912 
1913  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
1914  ExcMessage(
1915  "The amount of data that was read into new particles "
1916  "does not match the amount of data sent around."));
1917  }
1918 #endif
1919 
1920 
1921 
1922 #ifdef DEAL_II_WITH_MPI
1923  template <int dim, int spacedim>
1924  void
1926  const std::map<types::subdomain_id, std::vector<particle_iterator>>
1927  &particles_to_send)
1928  {
1929  const auto parallel_triangulation =
1930  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1931  &*triangulation);
1932  Assert(
1933  parallel_triangulation,
1934  ExcMessage(
1935  "This function is only implemented for parallel::TriangulationBase "
1936  "objects."));
1937 
1938  const auto &neighbors = ghost_particles_cache.neighbors;
1939  const auto &send_pointers = ghost_particles_cache.send_pointers;
1940  const auto &recv_pointers = ghost_particles_cache.recv_pointers;
1941 
1942  std::vector<char> &send_data = ghost_particles_cache.send_data;
1943 
1944  // Fill data to send
1945  if (send_pointers.back() > 0)
1946  {
1947  void *data = static_cast<void *>(&send_data.front());
1948 
1949  // Serialize the data sorted by receiving process
1950  for (const auto i : neighbors)
1951  for (const auto &p : particles_to_send.at(i))
1952  {
1953  data = p->write_particle_data_to_memory(data);
1954  if (store_callback)
1955  data = store_callback(p, data);
1956  }
1957  }
1958 
1959  std::vector<char> &recv_data = ghost_particles_cache.recv_data;
1960 
1961  // Exchange the particle data between domains
1962  {
1963  std::vector<MPI_Request> requests(2 * neighbors.size());
1964  unsigned int send_ops = 0;
1965  unsigned int recv_ops = 0;
1966 
1967  const int mpi_tag = Utilities::MPI::internal::Tags::
1969 
1970  for (unsigned int i = 0; i < neighbors.size(); ++i)
1971  if ((recv_pointers[i + 1] - recv_pointers[i]) > 0)
1972  {
1973  const int ierr =
1974  MPI_Irecv(recv_data.data() + recv_pointers[i],
1975  recv_pointers[i + 1] - recv_pointers[i],
1976  MPI_CHAR,
1977  neighbors[i],
1978  mpi_tag,
1979  parallel_triangulation->get_communicator(),
1980  &(requests[send_ops]));
1981  AssertThrowMPI(ierr);
1982  send_ops++;
1983  }
1984 
1985  for (unsigned int i = 0; i < neighbors.size(); ++i)
1986  if ((send_pointers[i + 1] - send_pointers[i]) > 0)
1987  {
1988  const int ierr =
1989  MPI_Isend(send_data.data() + send_pointers[i],
1990  send_pointers[i + 1] - send_pointers[i],
1991  MPI_CHAR,
1992  neighbors[i],
1993  mpi_tag,
1994  parallel_triangulation->get_communicator(),
1995  &(requests[send_ops + recv_ops]));
1996  AssertThrowMPI(ierr);
1997  recv_ops++;
1998  }
1999  const int ierr =
2000  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
2001  AssertThrowMPI(ierr);
2002  }
2003 
2004  // Put the received particles into the domain if they are in the
2005  // triangulation
2006  const void *recv_data_it = static_cast<const void *>(recv_data.data());
2007 
2008  // Gather ghost particle iterators from the cache
2009  auto &ghost_particles_iterators =
2010  ghost_particles_cache.ghost_particles_iterators;
2011 
2012  for (auto &recv_particle : ghost_particles_iterators)
2013  {
2014  // Update particle data using previously allocated memory space
2015  // for efficiency reasons
2016  recv_data_it =
2017  recv_particle->read_particle_data_from_memory(recv_data_it);
2018 
2019  Assert(recv_particle->particles_in_cell->cell->is_ghost(),
2020  ExcInternalError());
2021 
2022  if (load_callback)
2023  recv_data_it = load_callback(
2024  particle_iterator(recv_particle->particles_in_cell,
2025  *property_pool,
2026  recv_particle->particle_index_within_cell),
2027  recv_data_it);
2028  }
2029 
2030  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
2031  ExcMessage(
2032  "The amount of data that was read into new particles "
2033  "does not match the amount of data sent around."));
2034  }
2035 #endif
2036 
2037  template <int dim, int spacedim>
2038  void
2040  const std::function<std::size_t()> & size_callb,
2041  const std::function<void *(const particle_iterator &, void *)> &store_callb,
2042  const std::function<const void *(const particle_iterator &, const void *)>
2043  &load_callb)
2044  {
2045  size_callback = size_callb;
2046  store_callback = store_callb;
2047  load_callback = load_callb;
2048  }
2049 
2050 
2051  template <int dim, int spacedim>
2052  void
2054  {
2055  // First disconnect existing connections
2056  for (const auto &connection : tria_listeners)
2057  connection.disconnect();
2058 
2059  tria_listeners.clear();
2060 
2061  tria_listeners.push_back(triangulation->signals.create.connect([&]() {
2062  this->initialize(*(this->triangulation),
2063  *(this->mapping),
2064  this->property_pool->n_properties_per_slot());
2065  }));
2066 
2067  this->tria_listeners.push_back(
2068  this->triangulation->signals.clear.connect([&]() { this->clear(); }));
2069 
2070  // for distributed triangulations, connect to distributed signals
2072  *>(&(*triangulation)) != nullptr)
2073  {
2074  tria_listeners.push_back(
2075  triangulation->signals.post_distributed_refinement.connect(
2076  [&]() { this->post_mesh_change_action(); }));
2077  tria_listeners.push_back(
2078  triangulation->signals.post_distributed_repartition.connect(
2079  [&]() { this->post_mesh_change_action(); }));
2080  tria_listeners.push_back(
2081  triangulation->signals.post_distributed_load.connect(
2082  [&]() { this->post_mesh_change_action(); }));
2083  }
2084  else
2085  {
2086  tria_listeners.push_back(triangulation->signals.post_refinement.connect(
2087  [&]() { this->post_mesh_change_action(); }));
2088  }
2089  }
2090 
2091 
2092 
2093  template <int dim, int spacedim>
2094  void
2096  {
2097  Assert(triangulation != nullptr, ExcInternalError());
2098 
2099  const bool distributed_triangulation =
2100  dynamic_cast<
2102  &(*triangulation)) != nullptr;
2103  (void)distributed_triangulation;
2104 
2105  Assert(
2106  distributed_triangulation || number_of_locally_owned_particles == 0,
2107  ExcMessage(
2108  "Mesh refinement in a non-distributed triangulation is not supported "
2109  "by the ParticleHandler class. Either insert particles after mesh "
2110  "creation, or use a distributed triangulation."));
2111 
2112  // Resize the container if it is possible without
2113  // transferring particles
2114  if (number_of_locally_owned_particles == 0)
2115  cells_to_particle_cache.resize(triangulation->n_active_cells(),
2116  particles.end());
2117  }
2118 
2119 
2120 
2121  template <int dim, int spacedim>
2122  void
2124  {
2125  register_data_attach();
2126  }
2127 
2128 
2129 
2130  template <int dim, int spacedim>
2131  void
2133  {
2134  register_data_attach();
2135  }
2136 
2137 
2138 
2139  template <int dim, int spacedim>
2140  void
2142  {
2143  register_data_attach();
2144  }
2145 
2146 
2147 
2148  template <int dim, int spacedim>
2149  void
2151  {
2153  *distributed_triangulation =
2156  *>(&(*triangulation)));
2157  (void)distributed_triangulation;
2158 
2159  Assert(
2160  distributed_triangulation != nullptr,
2161  ExcMessage(
2162  "Mesh refinement in a non-distributed triangulation is not supported "
2163  "by the ParticleHandler class. Either insert particles after mesh "
2164  "creation and do not refine afterwards, or use a distributed triangulation."));
2165 
2166 #ifdef DEAL_II_WITH_P4EST
2167  const auto callback_function =
2168  [this](
2170  & cell_iterator,
2171  const typename Triangulation<dim, spacedim>::CellStatus cell_status) {
2172  return this->pack_callback(cell_iterator, cell_status);
2173  };
2174 
2175  handle = distributed_triangulation->register_data_attach(
2176  callback_function, /*returns_variable_size_data=*/true);
2177 #endif
2178  }
2179 
2180 
2181 
2182  template <int dim, int spacedim>
2183  void
2185  {
2186  const bool serialization = false;
2187  notify_ready_to_unpack(serialization);
2188  }
2189 
2190 
2191 
2192  template <int dim, int spacedim>
2193  void
2195  {
2196  const bool serialization = true;
2197  notify_ready_to_unpack(serialization);
2198  }
2199 
2200 
2201  template <int dim, int spacedim>
2202  void
2204  const bool serialization)
2205  {
2206  notify_ready_to_unpack(serialization);
2207  }
2208 
2209 
2210 
2211  template <int dim, int spacedim>
2212  void
2214  const bool serialization)
2215  {
2217  *distributed_triangulation =
2220  *>(&(*triangulation)));
2221  (void)distributed_triangulation;
2222 
2223  Assert(
2224  distributed_triangulation != nullptr,
2225  ExcMessage(
2226  "Mesh refinement in a non-distributed triangulation is not supported "
2227  "by the ParticleHandler class. Either insert particles after mesh "
2228  "creation and do not refine afterwards, or use a distributed triangulation."));
2229 
2230  // First prepare container for insertion
2231  clear();
2232 
2233 #ifdef DEAL_II_WITH_P4EST
2234  // If we are resuming from a checkpoint, we first have to register the
2235  // store function again, to set the triangulation to the same state as
2236  // before the serialization. Only afterwards we know how to deserialize the
2237  // data correctly.
2238  if (serialization)
2239  register_data_attach();
2240 
2241  // Check if something was stored and load it
2242  if (handle != numbers::invalid_unsigned_int)
2243  {
2244  const auto callback_function =
2245  [this](
2247  &cell_iterator,
2248  const typename Triangulation<dim, spacedim>::CellStatus cell_status,
2249  const boost::iterator_range<std::vector<char>::const_iterator>
2250  &range_iterator) {
2251  this->unpack_callback(cell_iterator, cell_status, range_iterator);
2252  };
2253 
2254  distributed_triangulation->notify_ready_to_unpack(handle,
2255  callback_function);
2256 
2257  // Reset handle and update global numbers.
2259  update_cached_numbers();
2260  }
2261 #else
2262  (void)serialization;
2263 #endif
2264  }
2265 
2266 
2267 
2268  template <int dim, int spacedim>
2269  std::vector<char>
2271  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2272  const typename Triangulation<dim, spacedim>::CellStatus status) const
2273  {
2274  std::vector<particle_iterator> stored_particles_on_cell;
2275 
2276  switch (status)
2277  {
2280  // If the cell persist or is refined store all particles of the
2281  // current cell.
2282  {
2283  const unsigned int n_particles = n_particles_in_cell(cell);
2284  stored_particles_on_cell.reserve(n_particles);
2285 
2286  for (unsigned int i = 0; i < n_particles; ++i)
2287  stored_particles_on_cell.push_back(particle_iterator(
2288  cells_to_particle_cache[cell->active_cell_index()],
2289  *property_pool,
2290  i));
2291  }
2292  break;
2293 
2295  // If this cell is the parent of children that will be coarsened,
2296  // collect the particles of all children.
2297  {
2298  for (const auto &child : cell->child_iterators())
2299  {
2300  const unsigned int n_particles = n_particles_in_cell(child);
2301 
2302  stored_particles_on_cell.reserve(
2303  stored_particles_on_cell.size() + n_particles);
2304 
2305  const typename particle_container::iterator &cache =
2306  cells_to_particle_cache[child->active_cell_index()];
2307  for (unsigned int i = 0; i < n_particles; ++i)
2308  stored_particles_on_cell.push_back(
2309  particle_iterator(cache, *property_pool, i));
2310  }
2311  }
2312  break;
2313 
2314  default:
2315  Assert(false, ExcInternalError());
2316  break;
2317  }
2318 
2319  return pack_particles(stored_particles_on_cell);
2320  }
2321 
2322 
2323 
2324  template <int dim, int spacedim>
2325  void
2327  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
2328  const typename Triangulation<dim, spacedim>::CellStatus status,
2329  const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
2330  {
2331  if (data_range.begin() == data_range.end())
2332  return;
2333 
2334  const auto cell_to_store_particles =
2336  cell :
2337  cell->child(0);
2338 
2339  // deserialize particles and insert into local storage
2340  if (data_range.begin() != data_range.end())
2341  {
2342  const void *data = static_cast<const void *>(&(*data_range.begin()));
2343  const void *end = static_cast<const void *>(
2344  &(*data_range.begin()) + (data_range.end() - data_range.begin()));
2345 
2346  while (data < end)
2347  insert_particle(data, cell_to_store_particles);
2348 
2349  Assert(data == end,
2350  ExcMessage(
2351  "The particle data could not be deserialized successfully. "
2352  "Check that when deserializing the particles you expect "
2353  "the same number of properties that were serialized."));
2354  }
2355 
2356  auto loaded_particles_on_cell = particles_in_cell(cell_to_store_particles);
2357 
2358  // now update particle storage location and properties if necessary
2359  switch (status)
2360  {
2362  {
2363  // all particles are correctly inserted
2364  }
2365  break;
2366 
2368  {
2369  // all particles are in correct cell, but their reference location
2370  // has changed
2371  for (auto &particle : loaded_particles_on_cell)
2372  {
2373  const Point<dim> p_unit =
2374  mapping->transform_real_to_unit_cell(cell_to_store_particles,
2375  particle.get_location());
2376  particle.set_reference_location(p_unit);
2377  }
2378  }
2379  break;
2380 
2382  {
2383  // we need to find the correct child to store the particles and
2384  // their reference location has changed
2385  typename particle_container::iterator &cache =
2386  cells_to_particle_cache[cell_to_store_particles
2387  ->active_cell_index()];
2388 
2389  // make sure that the call above has inserted an entry
2390  Assert(cache != particles.end(), ExcInternalError());
2391 
2392  // Cannot use range-based loop, because number of particles in cell
2393  // is going to change
2394  auto particle = loaded_particles_on_cell.begin();
2395  for (unsigned int i = 0; i < cache->particles.size();)
2396  {
2397  bool found_new_cell = false;
2398 
2399  for (unsigned int child_index = 0;
2400  child_index < GeometryInfo<dim>::max_children_per_cell;
2401  ++child_index)
2402  {
2404  child = cell->child(child_index);
2405  Assert(child->is_locally_owned(), ExcInternalError());
2406 
2407  try
2408  {
2409  const Point<dim> p_unit =
2410  mapping->transform_real_to_unit_cell(
2411  child, particle->get_location());
2413  1e-12))
2414  {
2415  found_new_cell = true;
2416  particle->set_reference_location(p_unit);
2417 
2418  // if the particle is not in child 0, we stored the
2419  // handle in the wrong place; move the handle and
2420  // redo the loop; otherwise move on to next particle
2421  if (child_index != 0)
2422  {
2423  insert_particle(cache->particles[i], child);
2424 
2425  cache->particles[i] = cache->particles.back();
2426  cache->particles.resize(
2427  cache->particles.size() - 1);
2428  }
2429  else
2430  {
2431  ++i;
2432  ++particle;
2433  }
2434  break;
2435  }
2436  }
2437  catch (typename Mapping<dim>::ExcTransformationFailed &)
2438  {}
2439  }
2440 
2441  if (found_new_cell == false)
2442  {
2443  // If we get here, we did not find the particle in any
2444  // child. This case may happen for particles that are at the
2445  // boundary for strongly curved cells. We apply a tolerance
2446  // in the call to GeometryInfo<dim>::is_inside_unit_cell to
2447  // account for this, but if that is not enough, we still
2448  // need to prevent an endless loop here. Delete the particle
2449  // and move on.
2450  signals.particle_lost(particle,
2451  particle->get_surrounding_cell());
2452  cache->particles[i] = cache->particles.back();
2453  cache->particles.resize(cache->particles.size() - 1);
2454  }
2455  }
2456  // clean up in case child 0 has no particle left
2457  if (cache->particles.empty())
2458  {
2459  particles.erase(cache);
2460  cache = particles.end();
2461  }
2462  }
2463  break;
2464 
2465  default:
2466  Assert(false, ExcInternalError());
2467  break;
2468  }
2469  }
2470 } // namespace Particles
2471 
2472 #include "particle_handler.inst"
2473 
std::size_t size() const
Definition: array_view.h:576
Definition: cell_id.h:71
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
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)
void exchange_ghost_particles(const bool enable_ghost_cache=false)
types::particle_index n_global_particles() const
unsigned int global_max_particles_per_cell
internal::GhostParticlePartitioner< dim, spacedim > ghost_particles_cache
particle_container::iterator particle_container_ghost_begin() const
unsigned int n_properties_per_particle() const
types::particle_index global_number_of_particles
void get_particle_positions(VectorType &output_vector, const bool add_to_output_vector=false)
particle_container::iterator particle_container_owned_end() const
particle_container::iterator particle_container_ghost_end() const
boost::iterator_range< particle_iterator > particle_iterator_range
types::particle_index number_of_locally_owned_particles
PropertyPool< dim, spacedim > & get_property_pool() const
void initialize(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const unsigned int n_properties=0)
std::unique_ptr< PropertyPool< dim, spacedim > > property_pool
types::particle_index get_max_local_particle_index() const
void reserve(const std::size_t n_particles)
void notify_ready_to_unpack(const bool serialization)
void register_load_callback_function(const bool serialization)
std::enable_if_t< std::is_convertible< VectorType *, Function< spacedim > * >::value==false > set_particle_positions(const VectorType &input_vector, const bool displace_particles=true)
types::particle_index next_free_particle_index
void send_recv_particles_properties_and_location(const std::map< types::subdomain_id, std::vector< particle_iterator >> &particles_to_send)
void unpack_callback(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)
void remove_particle(const particle_iterator &particle)
types::particle_index n_locally_owned_particles() const
void reset_particle_container(particle_container &particles)
types::particle_index n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
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={}, const std::vector< types::particle_index > &ids={})
typename ParticleAccessor< dim, spacedim >::particle_container particle_container
particle_iterator_range particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
SmartPointer< const Triangulation< dim, spacedim >, ParticleHandler< dim, spacedim > > triangulation
particle_iterator insert_particle(const Particle< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
void remove_particles(const std::vector< particle_iterator > &particles)
types::particle_index n_global_max_particles_per_cell() const
particle_container::iterator particle_container_owned_begin() const
void insert_particles(const std::multimap< typename Triangulation< dim, spacedim >::active_cell_iterator, Particle< dim, spacedim >> &particles)
void send_recv_particles(const std::map< types::subdomain_id, std::vector< particle_iterator >> &particles_to_send, 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 >>(), const bool enable_cache=false)
void copy_from(const ParticleHandler< dim, spacedim > &particle_handler)
SmartPointer< const Mapping< dim, spacedim >, ParticleHandler< dim, spacedim > > mapping
types::particle_index get_next_free_particle_index() const
particle_container particles
std::vector< char > pack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status) const
IndexSet locally_owned_particle_ids() const
void set_property_pool(PropertyPool< dim, spacedim > &property_pool)
Definition: particle.h:578
const Point< dim > & get_reference_location() const
Definition: particle.h:551
const Point< spacedim > & get_location() const
Definition: particle.h:533
std::size_t serialized_size_in_bytes() const
Definition: particle.cc:287
const ArrayView< double > get_properties()
Definition: particle.cc:330
types::particle_index get_id() const
Definition: particle.h:560
void sort_memory_slots(const std::vector< Handle > &handles_to_sort)
numbers::NumberTraits< Number >::real_type norm() const
Definition: vector.h:110
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_base.cc:798
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_base.cc:768
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
__global__ void set(Number *val, const Number s, const size_type N)
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1790
static ::ExceptionBase & ExcPointNotAvailableHere()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
TriaIterator< CellAccessor< dim, spacedim > > cell_iterator
Definition: tria.h:1355
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)
return_type compute_point_locations_try_all(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:5615
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, const double tolerance=1e-10)
Definition: grid_tools.cc:5783
@ past_the_end
Iterator reached end of container.
@ valid
Iterator points to a valid object.
types::global_dof_index size_type
Definition: cuda_kernels.h:45
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
VectorType::value_type * begin(VectorType &V)
VectorType::value_type * end(VectorType &V)
@ particle_handler_send_recv_particles_send
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:114
@ particle_handler_send_recv_particles_setup
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:110
std::map< unsigned int, T > some_to_some(const MPI_Comm &comm, const std::map< unsigned int, T > &objects_to_send)
T sum(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:145
std::vector< T > all_gather(const MPI_Comm &comm, const T &object_to_send)
T max(const T &t, const MPI_Comm &mpi_communicator)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13794
const types::subdomain_id artificial_subdomain_id
Definition: types.h:298
static const unsigned int invalid_unsigned_int
Definition: types.h:201
unsigned int subdomain_id
Definition: types.h:43
unsigned int particle_index
Definition: property_pool.h:65
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
Definition: property_pool.h:72
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
const MPI_Comm & comm