Reference documentation for deal.II version Git e67cfc96e3 2021-10-16 16:07:33 +0200
\(\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 - 2021 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  template <int dim, int spacedim>
56  : triangulation()
57  , mapping()
58  , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(0))
59  , particles()
60  , global_number_of_particles(0)
61  , number_of_locally_owned_particles(0)
62  , global_max_particles_per_cell(0)
63  , next_free_particle_index(0)
64  , size_callback()
65  , store_callback()
66  , load_callback()
67  , handle(numbers::invalid_unsigned_int)
68  , tria_listeners()
69  {}
70 
71 
72 
73  template <int dim, int spacedim>
77  const unsigned int n_properties)
78  : triangulation(&triangulation, typeid(*this).name())
79  , mapping(&mapping, typeid(*this).name())
80  , property_pool(std::make_unique<PropertyPool<dim, spacedim>>(n_properties))
81  , particles(triangulation.n_active_cells())
86  , size_callback()
87  , store_callback()
88  , load_callback()
91  std::make_unique<GridTools::Cache<dim, spacedim>>(triangulation,
92  mapping))
93  , tria_listeners()
94  {
96  }
97 
98 
99 
100  template <int dim, int spacedim>
102  {
103  clear_particles();
104 
105  for (const auto &connection : tria_listeners)
106  connection.disconnect();
107  }
108 
109 
110 
111  template <int dim, int spacedim>
112  void
114  const Triangulation<dim, spacedim> &new_triangulation,
115  const Mapping<dim, spacedim> & new_mapping,
116  const unsigned int n_properties)
117  {
118  clear();
119 
120  triangulation = &new_triangulation;
121  mapping = &new_mapping;
122 
123  // Create the memory pool that will store all particle properties
124  property_pool = std::make_unique<PropertyPool<dim, spacedim>>(n_properties);
125 
126  // Create the grid cache to cache the information about the triangulation
127  // that is used to locate the particles into subdomains and cells
129  std::make_unique<GridTools::Cache<dim, spacedim>>(new_triangulation,
130  new_mapping);
131 
132  particles.resize(triangulation->n_active_cells());
133 
135  }
136 
137 
138 
139  template <int dim, int spacedim>
140  void
142  const ParticleHandler<dim, spacedim> &particle_handler)
143  {
144  // clear this object before copying particles
145  clear();
146 
147  const unsigned int n_properties =
148  particle_handler.property_pool->n_properties_per_slot();
149  initialize(*particle_handler.triangulation,
150  *particle_handler.mapping,
151  n_properties);
152 
153  property_pool = std::make_unique<PropertyPool<dim, spacedim>>(
154  *(particle_handler.property_pool));
155 
156  // copy static members
159  particle_handler.number_of_locally_owned_particles;
160 
162  particle_handler.global_max_particles_per_cell;
164  particles = particle_handler.particles;
165 
166  ghost_particles_cache.ghost_particles_by_domain =
167  particle_handler.ghost_particles_cache.ghost_particles_by_domain;
168  handle = particle_handler.handle;
169  }
170 
171 
172 
173  template <int dim, int spacedim>
174  void
176  {
177  clear_particles();
182  }
183 
184 
185 
186  template <int dim, int spacedim>
187  void
189  {
190  for (auto &particles_in_cell : particles)
191  for (auto &particle : particles_in_cell)
193  property_pool->deregister_particle(particle);
194 
195  particles.clear();
196  if (triangulation != nullptr)
197  particles.resize(triangulation->n_active_cells());
198 
199  // the particle properties have already been deleted by their destructor,
200  // but the memory is still allocated. Return the memory as well.
201  property_pool->clear();
202  }
203 
204 
205 
206  template <int dim, int spacedim>
207  void
209  {
211  types::particle_index local_max_particle_index = 0;
212  types::particle_index local_max_particles_per_cell = 0;
213 
214  for (const auto &cell : triangulation->active_cell_iterators())
215  {
216  const auto &particles_in_cell = particles[cell->active_cell_index()];
217 
219  particles_in_cell.size();
220 
221  local_max_particles_per_cell =
222  std::max(local_max_particles_per_cell, n_particles_in_cell);
223 
224  if (cell->is_locally_owned())
225  number_of_locally_owned_particles += n_particles_in_cell;
226 
227  for (const auto &particle : particles_in_cell)
228  {
229  local_max_particle_index =
230  std::max(local_max_particle_index,
231  property_pool->get_id(particle));
232  }
233  }
234 
235  if (const auto parallel_triangulation =
236  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
237  &*triangulation))
238  {
241  parallel_triangulation->get_communicator());
242 
244  {
247  }
248  else
249  {
250  types::particle_index local_max_values[2] = {
251  local_max_particle_index, local_max_particles_per_cell};
252  types::particle_index global_max_values[2];
253 
254  Utilities::MPI::max(local_max_values,
255  parallel_triangulation->get_communicator(),
256  global_max_values);
257 
258  next_free_particle_index = global_max_values[0] + 1;
259  global_max_particles_per_cell = global_max_values[1];
260  }
261  }
262  else
263  {
266  global_number_of_particles == 0 ? 0 : local_max_particle_index + 1;
267  global_max_particles_per_cell = local_max_particles_per_cell;
268  }
269  }
270 
271 
272 
273  template <int dim, int spacedim>
277  const
278  {
279  if (particles.size() == 0)
280  return 0;
281 
282  if (cell->is_artificial() == false)
283  {
284  return particles[cell->active_cell_index()].size();
285  }
286  else
287  AssertThrow(false,
288  ExcMessage("You can't ask for the particles on an artificial "
289  "cell since we don't know what exists on these "
290  "kinds of cells."));
291 
293  }
294 
295 
296 
297  template <int dim, int spacedim>
301  const
302  {
303  return (const_cast<ParticleHandler<dim, spacedim> *>(this))
304  ->particles_in_cell(cell);
305  }
306 
307 
308 
309  template <int dim, int spacedim>
313  {
314  const unsigned int active_cell_index = cell->active_cell_index();
315 
316  if (cell->is_artificial() == false)
317  {
318  if (particles[active_cell_index].size() == 0)
319  {
323  }
324  else
325  {
328  *property_pool,
329  cell,
330  particles[active_cell_index].size() - 1);
331  // end needs to point to the particle after the last one in the
332  // cell.
333  ++end;
334 
335  return boost::make_iterator_range(begin, end);
336  }
337  }
338  else
339  AssertThrow(false,
340  ExcMessage("You can't ask for the particles on an artificial "
341  "cell since we don't know what exists on these "
342  "kinds of cells."));
343 
344  return {};
345  }
346 
347 
348 
349  template <int dim, int spacedim>
350  void
353  {
354  auto &particles_on_cell = *(particle->particles_on_cell);
355 
356  // if the particle has an invalid handle (e.g. because it has
357  // been duplicated before calling this function) do not try
358  // to deallocate its memory again
359  auto handle = particle->get_handle();
361  {
362  property_pool->deregister_particle(handle);
363  }
364 
365  // need to reduce the cached number before deleting, because the iterator
366  // may be invalid after removing the particle even if only
367  // accessing the cell
368  if (particle->get_surrounding_cell()->is_locally_owned())
370 
371  if (particles_on_cell.size() > 1)
372  {
373  particles_on_cell[particle->particle_index_within_cell] =
374  std::move(particles_on_cell.back());
375  particles_on_cell.resize(particles_on_cell.size() - 1);
376  }
377  else
378  {
379  particles_on_cell.clear();
380  }
381  }
382 
383 
384 
385  template <int dim, int spacedim>
386  void
389  &particles_to_remove)
390  {
391  unsigned int n_particles_removed = 0;
392 
393  for (auto particles_on_cell = particles.begin();
394  particles_on_cell != particles.end();
395  ++particles_on_cell)
396  {
397  // If there is nothing left to remove, break and return
398  if (n_particles_removed == particles_to_remove.size())
399  break;
400 
401  // Skip cells where there is nothing to remove
402  if (particles_to_remove[n_particles_removed]->particles_on_cell !=
403  particles_on_cell)
404  continue;
405 
406  const unsigned int n_particles_in_cell = particles_on_cell->size();
407  unsigned int move_to = 0;
408  for (unsigned int move_from = 0; move_from < n_particles_in_cell;
409  ++move_from)
410  {
411  if (n_particles_removed != particles_to_remove.size() &&
412  particles_to_remove[n_particles_removed]->particles_on_cell ==
413  particles_on_cell &&
414  particles_to_remove[n_particles_removed]
415  ->particle_index_within_cell == move_from)
416  {
417  auto handle =
418  particles_to_remove[n_particles_removed]->get_handle();
419 
420  // if some particles have invalid handles (e.g. because they are
421  // multiple times in the vector, or have been moved from
422  // before), do not try to deallocate their memory again
424  property_pool->deregister_particle(handle);
425  ++n_particles_removed;
426  continue;
427  }
428  else
429  {
430  (*particles_on_cell)[move_to] =
431  std::move((*particles_on_cell)[move_from]);
432  ++move_to;
433  }
434  }
435  particles_on_cell->resize(move_to);
436  }
437 
439  }
440 
441 
442 
443  template <int dim, int spacedim>
446  const Particle<dim, spacedim> & particle,
448  {
449  return insert_particle(particle.get_location(),
450  particle.get_reference_location(),
451  particle.get_id(),
452  cell,
453  particle.get_properties());
454  }
455 
456 
457 
458  template <int dim, int spacedim>
461  const void *& data,
463  {
464  Assert(triangulation != nullptr, ExcInternalError());
465  Assert(particles.size() == triangulation->n_active_cells(),
466  ExcInternalError());
467  Assert(
468  cell->is_locally_owned(),
469  ExcMessage(
470  "You can't insert particles in a cell that is not locally owned."));
471 
472  const unsigned int active_cell_index = cell->active_cell_index();
473  particles[active_cell_index].push_back(property_pool->register_particle());
474  particle_iterator particle_it(particles,
475  *property_pool,
476  cell,
477  particles[active_cell_index].size() - 1);
478 
479  data = particle_it->read_particle_data_from_memory(data);
480 
482 
483  return particle_it;
484  }
485 
486 
487 
488  template <int dim, int spacedim>
491  const Point<spacedim> & position,
492  const Point<dim> & reference_position,
495  const ArrayView<const double> &properties)
496  {
497  Assert(triangulation != nullptr, ExcInternalError());
498  Assert(particles.size() == triangulation->n_active_cells(),
499  ExcInternalError());
501  Assert(
502  cell->is_locally_owned(),
503  ExcMessage(
504  "You tried to insert a particle into a cell that is not locally owned. This is not supported."));
505 
506  const unsigned int active_cell_index = cell->active_cell_index();
507  particles[active_cell_index].push_back(property_pool->register_particle());
508  particle_iterator particle_it(particles,
509  *property_pool,
510  cell,
511  particles[active_cell_index].size() - 1);
512 
513  particle_it->set_location(position);
514  particle_it->set_reference_location(reference_position);
515  particle_it->set_id(particle_index);
516 
517  if (properties.size() != 0)
518  particle_it->set_properties(properties);
519 
521 
522  return particle_it;
523  }
524 
525 
526 
527  template <int dim, int spacedim>
528  void
530  const std::multimap<
532  Particle<dim, spacedim>> &new_particles)
533  {
534  for (const auto &cell_and_particle : new_particles)
535  insert_particle(cell_and_particle.second, cell_and_particle.first);
536 
538  }
539 
540 
541 
542  template <int dim, int spacedim>
543  void
545  const std::vector<Point<spacedim>> &positions)
546  {
547  Assert(triangulation != nullptr, ExcInternalError());
548 
550 
551  // Determine the starting particle index of this process, which
552  // is the highest currently existing particle index plus the sum
553  // of the number of newly generated particles of all
554  // processes with a lower rank if in a parallel computation.
555  const types::particle_index local_next_particle_index =
557  types::particle_index local_start_index = 0;
558 
559 #ifdef DEAL_II_WITH_MPI
560  if (const auto parallel_triangulation =
561  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
562  &*triangulation))
563  {
564  types::particle_index particles_to_add_locally = positions.size();
565  const int ierr = MPI_Scan(&particles_to_add_locally,
566  &local_start_index,
567  1,
569  MPI_SUM,
570  parallel_triangulation->get_communicator());
571  AssertThrowMPI(ierr);
572  local_start_index -= particles_to_add_locally;
573  }
574 #endif
575 
576  local_start_index += local_next_particle_index;
577 
578  auto point_locations =
580  positions);
581 
582  auto &cells = std::get<0>(point_locations);
583  auto &local_positions = std::get<1>(point_locations);
584  auto &index_map = std::get<2>(point_locations);
585  auto &missing_points = std::get<3>(point_locations);
586  // If a point was not found, throwing an error, as the old
587  // implementation of compute_point_locations would have done
588  AssertThrow(missing_points.size() == 0,
590 
591  (void)missing_points;
592 
593  for (unsigned int i = 0; i < cells.size(); ++i)
594  for (unsigned int p = 0; p < local_positions[i].size(); ++p)
595  insert_particle(positions[index_map[i][p]],
596  local_positions[i][p],
597  local_start_index + index_map[i][p],
598  cells[i]);
599 
601  }
602 
603 
604 
605  template <int dim, int spacedim>
606  std::map<unsigned int, IndexSet>
608  const std::vector<Point<spacedim>> &positions,
609  const std::vector<std::vector<BoundingBox<spacedim>>>
610  & global_bounding_boxes,
611  const std::vector<std::vector<double>> & properties,
612  const std::vector<types::particle_index> &ids)
613  {
614  if (!properties.empty())
615  {
616  AssertDimension(properties.size(), positions.size());
617 #ifdef DEBUG
618  for (const auto &p : properties)
620 #endif
621  }
622 
623  if (!ids.empty())
624  AssertDimension(ids.size(), positions.size());
625 
626  const auto comm = triangulation->get_communicator();
627 
629 
630  // Compute the global number of properties
631  const auto n_global_properties =
632  Utilities::MPI::sum(properties.size(), comm);
633 
634  // Gather the number of points per processor
635  const auto n_particles_per_proc =
636  Utilities::MPI::all_gather(comm, positions.size());
637 
638  // Calculate all starting points locally
639  std::vector<unsigned int> particle_start_indices(n_mpi_processes);
640 
641  unsigned int particle_start_index = get_next_free_particle_index();
642  for (unsigned int process = 0; process < particle_start_indices.size();
643  ++process)
644  {
645  particle_start_indices[process] = particle_start_index;
646  particle_start_index += n_particles_per_proc[process];
647  }
648 
649  // Get all local information
650  const auto cells_positions_and_index_maps =
652  positions,
653  global_bounding_boxes);
654 
655  // Unpack the information into several vectors:
656  // All cells that contain at least one particle
657  const auto &local_cells_containing_particles =
658  std::get<0>(cells_positions_and_index_maps);
659 
660  // The reference position of every particle in the local part of the
661  // triangulation.
662  const auto &local_reference_positions =
663  std::get<1>(cells_positions_and_index_maps);
664  // The original index in the positions vector for each particle in the
665  // local part of the triangulation
666  const auto &original_indices_of_local_particles =
667  std::get<2>(cells_positions_and_index_maps);
668  // The real spatial position of every particle in the local part of the
669  // triangulation.
670  const auto &local_positions = std::get<3>(cells_positions_and_index_maps);
671  // The MPI process that inserted each particle
672  const auto &calling_process_indices =
673  std::get<4>(cells_positions_and_index_maps);
674 
675  // Create the map of cpu to indices, indicating who sent us what particle
676  std::map<unsigned int, std::vector<unsigned int>>
677  original_process_to_local_particle_indices_tmp;
678  for (unsigned int i_cell = 0;
679  i_cell < local_cells_containing_particles.size();
680  ++i_cell)
681  {
682  for (unsigned int i_particle = 0;
683  i_particle < local_positions[i_cell].size();
684  ++i_particle)
685  {
686  const unsigned int local_id_on_calling_process =
687  original_indices_of_local_particles[i_cell][i_particle];
688  const unsigned int calling_process =
689  calling_process_indices[i_cell][i_particle];
690 
691  original_process_to_local_particle_indices_tmp[calling_process]
692  .push_back(local_id_on_calling_process);
693  }
694  }
695  std::map<unsigned int, IndexSet> original_process_to_local_particle_indices;
696  for (auto &process_and_particle_indices :
697  original_process_to_local_particle_indices_tmp)
698  {
699  const unsigned int calling_process = process_and_particle_indices.first;
700  original_process_to_local_particle_indices.insert(
701  {calling_process, IndexSet(n_particles_per_proc[calling_process])});
702  std::sort(process_and_particle_indices.second.begin(),
703  process_and_particle_indices.second.end());
704  original_process_to_local_particle_indices[calling_process].add_indices(
705  process_and_particle_indices.second.begin(),
706  process_and_particle_indices.second.end());
707  original_process_to_local_particle_indices[calling_process].compress();
708  }
709 
710  // A map from mpi process to properties, ordered as in the IndexSet.
711  // Notice that this ordering may be different from the ordering in the
712  // vectors above, since no local ordering is guaranteed by the
713  // distribute_compute_point_locations() call.
714  // This is only filled if n_global_properties is > 0
715  std::map<unsigned int, std::vector<std::vector<double>>>
716  locally_owned_properties_from_other_processes;
717 
718  // A map from mpi process to ids, ordered as in the IndexSet.
719  // Notice that this ordering may be different from the ordering in the
720  // vectors above, since no local ordering is guaranteed by the
721  // distribute_compute_point_locations() call.
722  // This is only filled if ids.size() is > 0
723  std::map<unsigned int, std::vector<types::particle_index>>
724  locally_owned_ids_from_other_processes;
725 
726  if (n_global_properties > 0 || !ids.empty())
727  {
728  // Gather whom I sent my own particles to, to decide whom to send
729  // the particle properties or the ids
730  auto send_to_cpu = Utilities::MPI::some_to_some(
731  comm, original_process_to_local_particle_indices);
732 
733  // Prepare the vector of properties to send
734  if (n_global_properties > 0)
735  {
736  std::map<unsigned int, std::vector<std::vector<double>>>
737  non_locally_owned_properties;
738 
739  for (const auto &it : send_to_cpu)
740  {
741  std::vector<std::vector<double>> properties_to_send(
742  it.second.n_elements(),
743  std::vector<double>(n_properties_per_particle()));
744  unsigned int index = 0;
745  for (const auto el : it.second)
746  properties_to_send[index++] = properties[el];
747  non_locally_owned_properties.insert(
748  {it.first, properties_to_send});
749  }
750 
751  // Send the non locally owned properties to each mpi process
752  // that needs them
753  locally_owned_properties_from_other_processes =
754  Utilities::MPI::some_to_some(comm, non_locally_owned_properties);
755 
757  locally_owned_properties_from_other_processes.size(),
758  original_process_to_local_particle_indices.size());
759  }
760 
761  if (!ids.empty())
762  {
763  std::map<unsigned int, std::vector<types::particle_index>>
764  non_locally_owned_ids;
765  for (const auto &it : send_to_cpu)
766  {
767  std::vector<types::particle_index> ids_to_send(
768  it.second.n_elements());
769  unsigned int index = 0;
770  for (const auto el : it.second)
771  ids_to_send[index++] = ids[el];
772  non_locally_owned_ids.insert({it.first, ids_to_send});
773  }
774 
775  // Send the non locally owned ids to each mpi process
776  // that needs them
777  locally_owned_ids_from_other_processes =
778  Utilities::MPI::some_to_some(comm, non_locally_owned_ids);
779 
780  AssertDimension(locally_owned_ids_from_other_processes.size(),
781  original_process_to_local_particle_indices.size());
782  }
783  }
784 
785  // Now fill up the actual particles
786  for (unsigned int i_cell = 0;
787  i_cell < local_cells_containing_particles.size();
788  ++i_cell)
789  {
790  for (unsigned int i_particle = 0;
791  i_particle < local_positions[i_cell].size();
792  ++i_particle)
793  {
794  const unsigned int local_id_on_calling_process =
795  original_indices_of_local_particles[i_cell][i_particle];
796 
797  const unsigned int calling_process =
798  calling_process_indices[i_cell][i_particle];
799 
800  const unsigned int index_within_set =
801  original_process_to_local_particle_indices[calling_process]
802  .index_within_set(local_id_on_calling_process);
803 
804  const unsigned int particle_id =
805  ids.empty() ?
806  local_id_on_calling_process +
807  particle_start_indices[calling_process] :
808  locally_owned_ids_from_other_processes[calling_process]
809  [index_within_set];
810 
811  auto particle_it =
812  insert_particle(local_positions[i_cell][i_particle],
813  local_reference_positions[i_cell][i_particle],
814  particle_id,
815  local_cells_containing_particles[i_cell]);
816 
817  if (n_global_properties > 0)
818  {
819  particle_it->set_properties(
820  locally_owned_properties_from_other_processes
821  [calling_process][index_within_set]);
822  }
823  }
824  }
825 
827 
828  return original_process_to_local_particle_indices;
829  }
830 
831 
832 
833  template <int dim, int spacedim>
834  std::map<unsigned int, IndexSet>
836  const std::vector<Particle<dim, spacedim>> &particles,
837  const std::vector<std::vector<BoundingBox<spacedim>>>
838  &global_bounding_boxes)
839  {
840  // Store the positions in a vector of points, the ids in a vector of ids,
841  // and the properties, if any, in a vector of vector of properties.
842  std::vector<Point<spacedim>> positions;
843  std::vector<std::vector<double>> properties;
844  std::vector<types::particle_index> ids;
845  positions.resize(particles.size());
846  ids.resize(particles.size());
847  if (n_properties_per_particle() > 0)
848  properties.resize(particles.size(),
849  std::vector<double>(n_properties_per_particle()));
850 
851  unsigned int i = 0;
852  for (const auto &p : particles)
853  {
854  positions[i] = p.get_location();
855  ids[i] = p.get_id();
856  if (p.has_properties())
857  properties[i] = {p.get_properties().begin(),
858  p.get_properties().end()};
859  ++i;
860  }
861 
862  return insert_global_particles(positions,
863  global_bounding_boxes,
864  properties,
865  ids);
866  }
867 
868 
869 
870  template <int dim, int spacedim>
873  {
875  }
876 
877 
878 
879  template <int dim, int spacedim>
882  {
884  }
885 
886 
887 
888  template <int dim, int spacedim>
891  {
893  }
894 
895 
896 
897  template <int dim, int spacedim>
898  unsigned int
900  {
901  return property_pool->n_properties_per_slot();
902  }
903 
904 
905 
906  template <int dim, int spacedim>
909  {
911  }
912 
913 
914 
915  template <int dim, int spacedim>
916  IndexSet
918  {
920  std::vector<types::particle_index> indices;
921  indices.reserve(n_locally_owned_particles());
922  for (const auto &p : *this)
923  indices.push_back(p.get_id());
924  set.add_indices(indices.begin(), indices.end());
925  set.compress();
926  return set;
927  }
928 
929 
930 
931  template <int dim, int spacedim>
934  {
935  return property_pool->n_slots();
936  }
937 
938 
939 
940  template <int dim, int spacedim>
941  void
943  std::vector<Point<spacedim>> &positions,
944  const bool add_to_output_vector)
945  {
946  // There should be one point per particle to gather
947  AssertDimension(positions.size(), n_locally_owned_particles());
948 
949  unsigned int i = 0;
950  for (auto it = begin(); it != end(); ++it, ++i)
951  {
952  if (add_to_output_vector)
953  positions[i] = positions[i] + it->get_location();
954  else
955  positions[i] = it->get_location();
956  }
957  }
958 
959 
960 
961  template <int dim, int spacedim>
962  void
964  const std::vector<Point<spacedim>> &new_positions,
965  const bool displace_particles)
966  {
967  // There should be one point per particle to fix the new position
968  AssertDimension(new_positions.size(), n_locally_owned_particles());
969 
970  unsigned int i = 0;
971  for (auto it = begin(); it != end(); ++it, ++i)
972  {
973  if (displace_particles)
974  it->set_location(it->get_location() + new_positions[i]);
975  else
976  it->set_location(new_positions[i]);
977  }
979  }
980 
981 
982 
983  template <int dim, int spacedim>
984  void
986  const Function<spacedim> &function,
987  const bool displace_particles)
988  {
989  // The function should have sufficient components to displace the
990  // particles
991  AssertDimension(function.n_components, spacedim);
992 
993  Vector<double> new_position(spacedim);
994  for (auto &particle : *this)
995  {
996  Point<spacedim> particle_location = particle.get_location();
997  function.vector_value(particle_location, new_position);
998  if (displace_particles)
999  for (unsigned int d = 0; d < spacedim; ++d)
1000  particle_location[d] += new_position[d];
1001  else
1002  for (unsigned int d = 0; d < spacedim; ++d)
1003  particle_location[d] = new_position[d];
1004  particle.set_location(particle_location);
1005  }
1007  }
1008 
1009 
1010 
1011  template <int dim, int spacedim>
1014  {
1015  return *property_pool;
1016  }
1017 
1018 
1019 
1020  namespace
1021  {
1031  template <int dim>
1032  bool
1033  compare_particle_association(
1034  const unsigned int a,
1035  const unsigned int b,
1036  const Tensor<1, dim> & particle_direction,
1037  const std::vector<Tensor<1, dim>> &center_directions)
1038  {
1039  const double scalar_product_a = center_directions[a] * particle_direction;
1040  const double scalar_product_b = center_directions[b] * particle_direction;
1041 
1042  // The function is supposed to return if a is before b. We are looking
1043  // for the alignment of particle direction and center direction,
1044  // therefore return if the scalar product of a is larger.
1045  return (scalar_product_a > scalar_product_b);
1046  }
1047  } // namespace
1048 
1049 
1050 
1051  template <int dim, int spacedim>
1052  void
1054  {
1055  Assert(triangulation != nullptr, ExcInternalError());
1056  Assert(particles.size() == triangulation->n_active_cells(),
1057  ExcInternalError());
1058 
1059  // TODO: The current algorithm only works for particles that are in
1060  // the local domain or in ghost cells, because it only knows the
1061  // subdomain_id of ghost cells, but not of artificial cells. This
1062  // can be extended using the distributed version of compute point
1063  // locations.
1064  // TODO: Extend this function to allow keeping particles on other
1065  // processes around (with an invalid cell).
1066 
1067  std::vector<particle_iterator> particles_out_of_cell;
1068 
1069  // Reserve some space for particles that need sorting to avoid frequent
1070  // re-allocation. Guess 25% of particles need sorting. Balance memory
1071  // overhead and performance.
1072  particles_out_of_cell.reserve(n_locally_owned_particles() / 4);
1073 
1074  // Now update the reference locations of the moved particles
1075  std::vector<Point<spacedim>> real_locations;
1076  std::vector<Point<dim>> reference_locations;
1077  real_locations.reserve(global_max_particles_per_cell);
1078  reference_locations.reserve(global_max_particles_per_cell);
1079 
1080  for (const auto &cell : triangulation->active_cell_iterators())
1081  {
1082  const unsigned int active_cell_index = cell->active_cell_index();
1083  const unsigned int n_pic = particles[active_cell_index].size();
1084 
1085  // Particles can be inserted into arbitrary cells, e.g. if their cell is
1086  // not known. However, for artificial cells we can not evaluate
1087  // the reference position of particles. Do not sort particles that are
1088  // not locally owned, because they will be sorted by the process that
1089  // owns them.
1090  if (cell->is_locally_owned() == false)
1091  {
1092  continue;
1093  }
1094 
1095  auto pic = particles_in_cell(cell);
1096 
1097  real_locations.clear();
1098  for (const auto &particle : pic)
1099  real_locations.push_back(particle.get_location());
1100 
1101  reference_locations.resize(n_pic);
1102  mapping->transform_points_real_to_unit_cell(cell,
1103  real_locations,
1104  reference_locations);
1105 
1106  auto particle = pic.begin();
1107  for (const auto &p_unit : reference_locations)
1108  {
1109  if (p_unit[0] == std::numeric_limits<double>::infinity() ||
1111  particles_out_of_cell.push_back(particle);
1112  else
1113  particle->set_reference_location(p_unit);
1114 
1115  ++particle;
1116  }
1117  }
1118 
1119  // There are three reasons why a particle is not in its old cell:
1120  // It moved to another cell, to another subdomain or it left the mesh.
1121  // Particles that moved to another cell are updated and moved inside the
1122  // particles vector, particles that moved to another domain are
1123  // collected in the moved_particles_domain vector. Particles that left
1124  // the mesh completely are ignored and removed.
1125  std::map<types::subdomain_id, std::vector<particle_iterator>>
1126  moved_particles;
1127  std::map<
1129  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1130  moved_cells;
1131 
1132  // We do not know exactly how many particles are lost, exchanged between
1133  // domains, or remain on this process. Therefore we pre-allocate
1134  // approximate sizes for these vectors. If more space is needed an
1135  // automatic and relatively fast (compared to other parts of this
1136  // algorithm) re-allocation will happen.
1137  using vector_size = typename std::vector<particle_iterator>::size_type;
1138 
1139  std::set<types::subdomain_id> ghost_owners;
1140  if (const auto parallel_triangulation =
1141  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1142  &*triangulation))
1143  ghost_owners = parallel_triangulation->ghost_owners();
1144 
1145  // Reserve some space for particles that need communication to avoid
1146  // frequent re-allocation. Guess 25% of particles out of their old cell need
1147  // communication. Balance memory overhead and performance.
1148  for (const auto &ghost_owner : ghost_owners)
1149  moved_particles[ghost_owner].reserve(particles_out_of_cell.size() / 4);
1150  for (const auto &ghost_owner : ghost_owners)
1151  moved_cells[ghost_owner].reserve(particles_out_of_cell.size() / 4);
1152 
1153  {
1154  // Create a map from vertices to adjacent cells using grid cache
1155  const std::vector<
1156  std::set<typename Triangulation<dim, spacedim>::active_cell_iterator>>
1157  &vertex_to_cells = triangulation_cache->get_vertex_to_cell_map();
1158 
1159  // Create a corresponding map of vectors from vertex to cell center
1160  // using grid cache
1161  const std::vector<std::vector<Tensor<1, spacedim>>>
1162  &vertex_to_cell_centers =
1163  triangulation_cache->get_vertex_to_cell_centers_directions();
1164 
1165  std::vector<unsigned int> neighbor_permutation;
1166 
1167  // Find the cells that the particles moved to.
1168  for (auto &out_particle : particles_out_of_cell)
1169  {
1170  // make a copy of the current cell, since we will modify the
1171  // variable current_cell in the following but we need the original in
1172  // the case the particle is not found
1173  auto current_cell = out_particle->get_surrounding_cell();
1174 
1175  // The cell the particle is in
1176  Point<dim> current_reference_position;
1177  bool found_cell = false;
1178 
1179  // Check if the particle is in one of the old cell's neighbors
1180  // that are adjacent to the closest vertex
1181  const unsigned int closest_vertex =
1182  GridTools::find_closest_vertex_of_cell<dim, spacedim>(
1183  current_cell, out_particle->get_location(), *mapping);
1184  Tensor<1, spacedim> vertex_to_particle =
1185  out_particle->get_location() - current_cell->vertex(closest_vertex);
1186  vertex_to_particle /= vertex_to_particle.norm();
1187 
1188  const unsigned int closest_vertex_index =
1189  current_cell->vertex_index(closest_vertex);
1190  const unsigned int n_neighbor_cells =
1191  vertex_to_cells[closest_vertex_index].size();
1192 
1193  neighbor_permutation.resize(n_neighbor_cells);
1194  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1195  neighbor_permutation[i] = i;
1196 
1197  const auto &cell_centers =
1198  vertex_to_cell_centers[closest_vertex_index];
1199  std::sort(neighbor_permutation.begin(),
1200  neighbor_permutation.end(),
1201  [&vertex_to_particle, &cell_centers](const unsigned int a,
1202  const unsigned int b) {
1203  return compare_particle_association(a,
1204  b,
1205  vertex_to_particle,
1206  cell_centers);
1207  });
1208 
1209  // Search all of the cells adjacent to the closest vertex of the
1210  // previous cell Most likely we will find the particle in them.
1211  for (unsigned int i = 0; i < n_neighbor_cells; ++i)
1212  {
1213  try
1214  {
1215  typename std::set<typename Triangulation<dim, spacedim>::
1216  active_cell_iterator>::const_iterator
1217  cell = vertex_to_cells[closest_vertex_index].begin();
1218 
1219  std::advance(cell, neighbor_permutation[i]);
1220  const Point<dim> p_unit =
1221  mapping->transform_real_to_unit_cell(
1222  *cell, out_particle->get_location());
1224  {
1225  current_cell = *cell;
1226  current_reference_position = p_unit;
1227  found_cell = true;
1228  break;
1229  }
1230  }
1231  catch (typename Mapping<dim>::ExcTransformationFailed &)
1232  {}
1233  }
1234 
1235  if (!found_cell)
1236  {
1237  // The particle is not in a neighbor of the old cell.
1238  // Look for the new cell in the whole local domain.
1239  // This case is rare.
1240  const std::pair<const typename Triangulation<dim, spacedim>::
1241  active_cell_iterator,
1242  Point<dim>>
1243  current_cell_and_position =
1244  GridTools::find_active_cell_around_point<>(
1245  *triangulation_cache, out_particle->get_location());
1246  current_cell = current_cell_and_position.first;
1247  current_reference_position = current_cell_and_position.second;
1248 
1249  if (current_cell.state() != IteratorState::valid)
1250  {
1251  // We can find no cell for this particle. It has left the
1252  // domain due to an integration error or an open boundary.
1253  // Signal the loss and move on.
1254  signals.particle_lost(out_particle,
1255  out_particle->get_surrounding_cell());
1256  continue;
1257  }
1258  }
1259 
1260  // If we are here, we found a cell and reference position for this
1261  // particle
1262  out_particle->set_reference_location(current_reference_position);
1263 
1264  // Reinsert the particle into our domain if we own its cell.
1265  // Mark it for MPI transfer otherwise
1266  if (current_cell->is_locally_owned())
1267  {
1268  const unsigned int active_cell_index =
1269  current_cell->active_cell_index();
1270 
1271  particles[active_cell_index].push_back(
1272  (*out_particle->particles_on_cell)
1273  [out_particle->particle_index_within_cell]);
1274 
1275  // Avoid deallocating the memory of this particle
1276  (*out_particle->particles_on_cell)
1277  [out_particle->particle_index_within_cell] =
1279  }
1280  else
1281  {
1282  moved_particles[current_cell->subdomain_id()].push_back(
1283  out_particle);
1284  moved_cells[current_cell->subdomain_id()].push_back(current_cell);
1285  }
1286  }
1287  }
1288 
1289  // Exchange particles between processors if we have more than one process
1290 #ifdef DEAL_II_WITH_MPI
1291  if (const auto parallel_triangulation =
1292  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1293  &*triangulation))
1294  {
1296  parallel_triangulation->get_communicator()) > 1)
1297  send_recv_particles(moved_particles, particles, moved_cells);
1298  }
1299 #endif
1300 
1301  // remove_particles also calls update_cached_numbers()
1302  remove_particles(particles_out_of_cell);
1303 
1304  // now make sure particle data is sorted in order of iteration
1305  std::vector<typename PropertyPool<dim, spacedim>::Handle> unsorted_handles;
1306  unsorted_handles.reserve(property_pool->n_registered_slots());
1307 
1308  typename PropertyPool<dim, spacedim>::Handle sorted_handle = 0;
1309  for (auto &particles_in_cell : particles)
1310  for (auto &particle : particles_in_cell)
1311  {
1312  unsorted_handles.push_back(particle);
1313  particle = sorted_handle++;
1314  }
1315 
1316  property_pool->sort_memory_slots(unsorted_handles);
1317 
1318  } // namespace Particles
1319 
1320 
1321 
1322  template <int dim, int spacedim>
1323  void
1325  const bool enable_cache)
1326  {
1327  // Nothing to do in serial computations
1328  const auto parallel_triangulation =
1329  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1330  &*triangulation);
1331  if (parallel_triangulation != nullptr)
1332  {
1334  parallel_triangulation->get_communicator()) == 1)
1335  return;
1336  }
1337  else
1338  return;
1339 
1340 #ifndef DEAL_II_WITH_MPI
1341  (void)enable_cache;
1342 #else
1343  // Clear ghost particles and their properties
1344  for (const auto &cell : triangulation->active_cell_iterators())
1345  if (cell->is_ghost())
1346  {
1347  // Clear particle properties
1348  for (auto &ghost_particle : particles[cell->active_cell_index()])
1349  property_pool->deregister_particle(ghost_particle);
1350 
1351  // Clear particles themselves
1352  particles[cell->active_cell_index()].clear();
1353  }
1354 
1355  // Clear ghost particles cache and invalidate it
1356  ghost_particles_cache.ghost_particles_by_domain.clear();
1357  ghost_particles_cache.valid = false;
1358 
1359 
1360  const std::set<types::subdomain_id> ghost_owners =
1361  parallel_triangulation->ghost_owners();
1362  for (const auto ghost_owner : ghost_owners)
1363  ghost_particles_cache.ghost_particles_by_domain[ghost_owner].reserve(
1365 
1366  const std::vector<std::set<unsigned int>> vertex_to_neighbor_subdomain =
1367  triangulation_cache->get_vertex_to_neighbor_subdomain();
1368 
1369  for (const auto &cell : triangulation->active_cell_iterators())
1370  {
1371  if (cell->is_locally_owned())
1372  {
1373  std::set<unsigned int> cell_to_neighbor_subdomain;
1374  for (const unsigned int v : cell->vertex_indices())
1375  {
1376  cell_to_neighbor_subdomain.insert(
1377  vertex_to_neighbor_subdomain[cell->vertex_index(v)].begin(),
1378  vertex_to_neighbor_subdomain[cell->vertex_index(v)].end());
1379  }
1380 
1381  if (cell_to_neighbor_subdomain.size() > 0)
1382  {
1383  const particle_iterator_range particle_range =
1384  particles_in_cell(cell);
1385 
1386  for (const auto domain : cell_to_neighbor_subdomain)
1387  {
1388  for (typename particle_iterator_range::iterator particle =
1389  particle_range.begin();
1390  particle != particle_range.end();
1391  ++particle)
1392  ghost_particles_cache.ghost_particles_by_domain[domain]
1393  .push_back(particle);
1394  }
1395  }
1396  }
1397  }
1398 
1400  ghost_particles_cache.ghost_particles_by_domain,
1401  particles,
1402  std::map<
1404  std::vector<
1406  enable_cache);
1407 #endif
1408  }
1409 
1410  template <int dim, int spacedim>
1411  void
1413  {
1414  // Nothing to do in serial computations
1415  const auto parallel_triangulation =
1416  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1417  &*triangulation);
1418  if (parallel_triangulation == nullptr ||
1420  parallel_triangulation->get_communicator()) == 1)
1421  {
1422  return;
1423  }
1424 
1425 
1426 #ifdef DEAL_II_WITH_MPI
1427  // First clear the current ghost_particle information
1428  // ghost_particles.clear();
1430  ExcMessage(
1431  "Ghost particles cannot be updated if they first have not been "
1432  "exchanged at least once with the cache enabled"));
1433 
1434 
1436  ghost_particles_cache.ghost_particles_by_domain, particles);
1437 #endif
1438  }
1439 
1440 
1441 
1442 #ifdef DEAL_II_WITH_MPI
1443  template <int dim, int spacedim>
1444  void
1446  const std::map<types::subdomain_id, std::vector<particle_iterator>>
1447  & particles_to_send,
1448  particle_container &received_particles,
1449  const std::map<
1452  & send_cells,
1453  const bool build_cache)
1454  {
1455  Assert(triangulation != nullptr, ExcInternalError());
1456  Assert(received_particles.size() == triangulation->n_active_cells(),
1457  ExcInternalError());
1458 
1459  ghost_particles_cache.valid = build_cache;
1460 
1461  const auto parallel_triangulation =
1462  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1463  &*triangulation);
1464  Assert(
1465  parallel_triangulation,
1466  ExcMessage(
1467  "This function is only implemented for parallel::TriangulationBase "
1468  "objects."));
1469 
1470  // Determine the communication pattern
1471  const std::set<types::subdomain_id> ghost_owners =
1472  parallel_triangulation->ghost_owners();
1473  const std::vector<types::subdomain_id> neighbors(ghost_owners.begin(),
1474  ghost_owners.end());
1475  const unsigned int n_neighbors = neighbors.size();
1476 
1477  if (send_cells.size() != 0)
1478  Assert(particles_to_send.size() == send_cells.size(), ExcInternalError());
1479 
1480  // If we do not know the subdomain this particle needs to be send to,
1481  // throw an error
1482  Assert(particles_to_send.find(numbers::artificial_subdomain_id) ==
1483  particles_to_send.end(),
1484  ExcInternalError());
1485 
1486  // TODO: Implement the shipping of particles to processes that are not
1487  // ghost owners of the local domain
1488  for (auto send_particles = particles_to_send.begin();
1489  send_particles != particles_to_send.end();
1490  ++send_particles)
1491  Assert(ghost_owners.find(send_particles->first) != ghost_owners.end(),
1492  ExcNotImplemented());
1493 
1494  std::size_t n_send_particles = 0;
1495  for (auto send_particles = particles_to_send.begin();
1496  send_particles != particles_to_send.end();
1497  ++send_particles)
1498  n_send_particles += send_particles->second.size();
1499 
1500  const unsigned int cellid_size = sizeof(CellId::binary_type);
1501 
1502  // Containers for the amount and offsets of data we will send
1503  // to other processors and the data itself.
1504  std::vector<unsigned int> n_send_data(n_neighbors, 0);
1505  std::vector<unsigned int> send_offsets(n_neighbors, 0);
1506  std::vector<char> send_data;
1507 
1508  Particle<dim, spacedim> test_particle;
1509  test_particle.set_property_pool(*property_pool);
1510 
1511  const unsigned int individual_particle_data_size =
1512  test_particle.serialized_size_in_bytes() +
1513  (size_callback ? size_callback() : 0);
1514 
1515  const unsigned int individual_total_particle_data_size =
1516  individual_particle_data_size + cellid_size;
1517 
1518  // Only serialize things if there are particles to be send.
1519  // We can not return early even if no particles
1520  // are send, because we might receive particles from other processes
1521  if (n_send_particles > 0)
1522  {
1523  // Allocate space for sending particle data
1524  send_data.resize(n_send_particles *
1525  individual_total_particle_data_size);
1526 
1527  void *data = static_cast<void *>(&send_data.front());
1528 
1529  // Serialize the data sorted by receiving process
1530  for (unsigned int i = 0; i < n_neighbors; ++i)
1531  {
1532  send_offsets[i] = reinterpret_cast<std::size_t>(data) -
1533  reinterpret_cast<std::size_t>(&send_data.front());
1534 
1535  const unsigned int n_particles_to_send =
1536  particles_to_send.at(neighbors[i]).size();
1537 
1538  Assert(static_cast<std::size_t>(n_particles_to_send) *
1539  individual_total_particle_data_size ==
1540  static_cast<std::size_t>(
1541  n_particles_to_send *
1542  individual_total_particle_data_size),
1543  ExcMessage("Overflow when trying to send particle "
1544  "data"));
1545 
1546  for (unsigned int j = 0; j < n_particles_to_send; ++j)
1547  {
1548  // If no target cells are given, use the iterator
1549  // information
1551  cell;
1552  if (send_cells.size() == 0)
1553  cell = particles_to_send.at(neighbors[i])[j]
1554  ->get_surrounding_cell();
1555  else
1556  cell = send_cells.at(neighbors[i])[j];
1557 
1558  const CellId::binary_type cellid =
1559  cell->id().template to_binary<dim>();
1560  memcpy(data, &cellid, cellid_size);
1561  data = static_cast<char *>(data) + cellid_size;
1562 
1563  data = particles_to_send.at(neighbors[i])[j]
1564  ->write_particle_data_to_memory(data);
1565  if (store_callback)
1566  data =
1567  store_callback(particles_to_send.at(neighbors[i])[j], data);
1568  }
1569  n_send_data[i] = n_particles_to_send;
1570  }
1571  }
1572 
1573  // Containers for the data we will receive from other processors
1574  std::vector<unsigned int> n_recv_data(n_neighbors);
1575  std::vector<unsigned int> recv_offsets(n_neighbors);
1576 
1577  {
1578  const int mpi_tag = Utilities::MPI::internal::Tags::
1580 
1581  std::vector<MPI_Request> n_requests(2 * n_neighbors);
1582  for (unsigned int i = 0; i < n_neighbors; ++i)
1583  {
1584  const int ierr = MPI_Irecv(&(n_recv_data[i]),
1585  1,
1586  MPI_UNSIGNED,
1587  neighbors[i],
1588  mpi_tag,
1589  parallel_triangulation->get_communicator(),
1590  &(n_requests[2 * i]));
1591  AssertThrowMPI(ierr);
1592  }
1593  for (unsigned int i = 0; i < n_neighbors; ++i)
1594  {
1595  const int ierr = MPI_Isend(&(n_send_data[i]),
1596  1,
1597  MPI_UNSIGNED,
1598  neighbors[i],
1599  mpi_tag,
1600  parallel_triangulation->get_communicator(),
1601  &(n_requests[2 * i + 1]));
1602  AssertThrowMPI(ierr);
1603  }
1604  const int ierr =
1605  MPI_Waitall(2 * n_neighbors, n_requests.data(), MPI_STATUSES_IGNORE);
1606  AssertThrowMPI(ierr);
1607  }
1608 
1609  // Determine how many particles and data we will receive
1610  unsigned int total_recv_data = 0;
1611  for (unsigned int neighbor_id = 0; neighbor_id < n_neighbors; ++neighbor_id)
1612  {
1613  recv_offsets[neighbor_id] = total_recv_data;
1614  total_recv_data +=
1615  n_recv_data[neighbor_id] * individual_total_particle_data_size;
1616  }
1617 
1618  // Set up the space for the received particle data
1619  std::vector<char> recv_data(total_recv_data);
1620 
1621  // Exchange the particle data between domains
1622  {
1623  std::vector<MPI_Request> requests(2 * n_neighbors);
1624  unsigned int send_ops = 0;
1625  unsigned int recv_ops = 0;
1626 
1627  const int mpi_tag = Utilities::MPI::internal::Tags::
1629 
1630  for (unsigned int i = 0; i < n_neighbors; ++i)
1631  if (n_recv_data[i] > 0)
1632  {
1633  const int ierr =
1634  MPI_Irecv(&(recv_data[recv_offsets[i]]),
1635  n_recv_data[i] * individual_total_particle_data_size,
1636  MPI_CHAR,
1637  neighbors[i],
1638  mpi_tag,
1639  parallel_triangulation->get_communicator(),
1640  &(requests[send_ops]));
1641  AssertThrowMPI(ierr);
1642  send_ops++;
1643  }
1644 
1645  for (unsigned int i = 0; i < n_neighbors; ++i)
1646  if (n_send_data[i] > 0)
1647  {
1648  const int ierr =
1649  MPI_Isend(&(send_data[send_offsets[i]]),
1650  n_send_data[i] * individual_total_particle_data_size,
1651  MPI_CHAR,
1652  neighbors[i],
1653  mpi_tag,
1654  parallel_triangulation->get_communicator(),
1655  &(requests[send_ops + recv_ops]));
1656  AssertThrowMPI(ierr);
1657  recv_ops++;
1658  }
1659  const int ierr =
1660  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
1661  AssertThrowMPI(ierr);
1662  }
1663 
1664  // Put the received particles into the domain if they are in the
1665  // triangulation
1666  const void *recv_data_it = static_cast<const void *>(recv_data.data());
1667 
1668  // Store the particle iterators in the cache
1669  auto &ghost_particles_iterators =
1670  ghost_particles_cache.ghost_particles_iterators;
1671 
1672  if (build_cache)
1673  {
1674  ghost_particles_iterators.clear();
1675 
1676  auto &send_pointers_particles = ghost_particles_cache.send_pointers;
1677  send_pointers_particles.assign(n_neighbors + 1, 0);
1678 
1679  for (unsigned int i = 0; i < n_neighbors; ++i)
1680  send_pointers_particles[i + 1] =
1681  send_pointers_particles[i] +
1682  n_send_data[i] * individual_particle_data_size;
1683 
1684  auto &recv_pointers_particles = ghost_particles_cache.recv_pointers;
1685  recv_pointers_particles.assign(n_neighbors + 1, 0);
1686 
1687  for (unsigned int i = 0; i < n_neighbors; ++i)
1688  recv_pointers_particles[i + 1] =
1689  recv_pointers_particles[i] +
1690  n_recv_data[i] * individual_particle_data_size;
1691 
1692  ghost_particles_cache.neighbors = neighbors;
1693 
1694  ghost_particles_cache.send_data.resize(
1695  ghost_particles_cache.send_pointers.back());
1696  ghost_particles_cache.recv_data.resize(
1697  ghost_particles_cache.recv_pointers.back());
1698  }
1699 
1700  while (reinterpret_cast<std::size_t>(recv_data_it) -
1701  reinterpret_cast<std::size_t>(recv_data.data()) <
1702  total_recv_data)
1703  {
1704  CellId::binary_type binary_cellid;
1705  memcpy(&binary_cellid, recv_data_it, cellid_size);
1706  const CellId id(binary_cellid);
1707  recv_data_it = static_cast<const char *>(recv_data_it) + cellid_size;
1708 
1710  triangulation->create_cell_iterator(id);
1711 
1712  const unsigned int active_cell_index = cell->active_cell_index();
1713  received_particles[active_cell_index].emplace_back(
1714  property_pool->register_particle());
1715  particle_iterator particle_it(
1716  received_particles,
1717  *property_pool,
1718  cell,
1719  received_particles[active_cell_index].size() - 1);
1720 
1721  recv_data_it =
1722  particle_it->read_particle_data_from_memory(recv_data_it);
1723 
1724  if (load_callback)
1725  recv_data_it = load_callback(particle_it, recv_data_it);
1726 
1727  if (build_cache) // TODO: is this safe?
1728  ghost_particles_iterators.push_back(particle_it);
1729  }
1730 
1731  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
1732  ExcMessage(
1733  "The amount of data that was read into new particles "
1734  "does not match the amount of data sent around."));
1735  }
1736 #endif
1737 
1738 
1739 
1740 #ifdef DEAL_II_WITH_MPI
1741  template <int dim, int spacedim>
1742  void
1744  const std::map<types::subdomain_id, std::vector<particle_iterator>>
1745  & particles_to_send,
1746  particle_container &updated_particles)
1747  {
1748  const auto parallel_triangulation =
1749  dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
1750  &*triangulation);
1751  Assert(
1752  parallel_triangulation,
1753  ExcMessage(
1754  "This function is only implemented for parallel::TriangulationBase "
1755  "objects."));
1756 
1757  Assert(updated_particles.size() == parallel_triangulation->n_active_cells(),
1758  ExcInternalError());
1759 
1760  const auto &neighbors = ghost_particles_cache.neighbors;
1761  const auto &send_pointers = ghost_particles_cache.send_pointers;
1762  const auto &recv_pointers = ghost_particles_cache.recv_pointers;
1763 
1764  std::vector<char> &send_data = ghost_particles_cache.send_data;
1765 
1766  // Fill data to send
1767  if (send_pointers.back() > 0)
1768  {
1769  void *data = static_cast<void *>(&send_data.front());
1770 
1771  // Serialize the data sorted by receiving process
1772  for (const auto i : neighbors)
1773  for (const auto &p : particles_to_send.at(i))
1774  {
1775  data = p->write_particle_data_to_memory(data);
1776  if (store_callback)
1777  data = store_callback(p, data);
1778  }
1779  }
1780 
1781  std::vector<char> &recv_data = ghost_particles_cache.recv_data;
1782 
1783  // Exchange the particle data between domains
1784  {
1785  std::vector<MPI_Request> requests(2 * neighbors.size());
1786  unsigned int send_ops = 0;
1787  unsigned int recv_ops = 0;
1788 
1789  const int mpi_tag = Utilities::MPI::internal::Tags::
1791 
1792  for (unsigned int i = 0; i < neighbors.size(); ++i)
1793  if ((recv_pointers[i + 1] - recv_pointers[i]) > 0)
1794  {
1795  const int ierr =
1796  MPI_Irecv(recv_data.data() + recv_pointers[i],
1797  recv_pointers[i + 1] - recv_pointers[i],
1798  MPI_CHAR,
1799  neighbors[i],
1800  mpi_tag,
1801  parallel_triangulation->get_communicator(),
1802  &(requests[send_ops]));
1803  AssertThrowMPI(ierr);
1804  send_ops++;
1805  }
1806 
1807  for (unsigned int i = 0; i < neighbors.size(); ++i)
1808  if ((send_pointers[i + 1] - send_pointers[i]) > 0)
1809  {
1810  const int ierr =
1811  MPI_Isend(send_data.data() + send_pointers[i],
1812  send_pointers[i + 1] - send_pointers[i],
1813  MPI_CHAR,
1814  neighbors[i],
1815  mpi_tag,
1816  parallel_triangulation->get_communicator(),
1817  &(requests[send_ops + recv_ops]));
1818  AssertThrowMPI(ierr);
1819  recv_ops++;
1820  }
1821  const int ierr =
1822  MPI_Waitall(send_ops + recv_ops, requests.data(), MPI_STATUSES_IGNORE);
1823  AssertThrowMPI(ierr);
1824  }
1825 
1826  // Put the received particles into the domain if they are in the
1827  // triangulation
1828  const void *recv_data_it = static_cast<const void *>(recv_data.data());
1829 
1830  // Gather ghost particle iterators from the cache
1831  auto &ghost_particles_iterators =
1832  ghost_particles_cache.ghost_particles_iterators;
1833 
1834  for (auto &recv_particle : ghost_particles_iterators)
1835  {
1836  // Update particle data using previously allocated memory space
1837  // for efficiency reasons
1838  recv_data_it =
1839  recv_particle->read_particle_data_from_memory(recv_data_it);
1840 
1841  if (load_callback)
1842  recv_data_it = load_callback(
1843  particle_iterator(updated_particles,
1844  *property_pool,
1845  recv_particle->cell,
1846  recv_particle->particle_index_within_cell),
1847  recv_data_it);
1848  }
1849 
1850  AssertThrow(recv_data_it == recv_data.data() + recv_data.size(),
1851  ExcMessage(
1852  "The amount of data that was read into new particles "
1853  "does not match the amount of data sent around."));
1854  }
1855 #endif
1856 
1857  template <int dim, int spacedim>
1858  void
1860  const std::function<std::size_t()> & size_callb,
1861  const std::function<void *(const particle_iterator &, void *)> &store_callb,
1862  const std::function<const void *(const particle_iterator &, const void *)>
1863  &load_callb)
1864  {
1865  size_callback = size_callb;
1866  store_callback = store_callb;
1867  load_callback = load_callb;
1868  }
1869 
1870 
1871  template <int dim, int spacedim>
1872  void
1874  {
1875  // First disconnect existing connections
1876  for (const auto &connection : tria_listeners)
1877  connection.disconnect();
1878 
1879  tria_listeners.clear();
1880 
1881  tria_listeners.push_back(triangulation->signals.create.connect([&]() {
1882  this->initialize(*(this->triangulation),
1883  *(this->mapping),
1884  this->property_pool->n_properties_per_slot());
1885  }));
1886 
1887  this->tria_listeners.push_back(
1888  this->triangulation->signals.clear.connect([&]() { this->clear(); }));
1889 
1890  // for distributed triangulations, connect to distributed signals
1892  *>(&(*triangulation)) != nullptr)
1893  {
1894  tria_listeners.push_back(
1895  triangulation->signals.post_distributed_refinement.connect(
1896  [&]() { this->post_mesh_change_action(); }));
1897  tria_listeners.push_back(
1898  triangulation->signals.post_distributed_repartition.connect(
1899  [&]() { this->post_mesh_change_action(); }));
1900  tria_listeners.push_back(
1901  triangulation->signals.post_distributed_load.connect(
1902  [&]() { this->post_mesh_change_action(); }));
1903  }
1904  else
1905  {
1906  tria_listeners.push_back(triangulation->signals.post_refinement.connect(
1907  [&]() { this->post_mesh_change_action(); }));
1908  }
1909  }
1910 
1911 
1912 
1913  template <int dim, int spacedim>
1914  void
1916  {
1917  Assert(triangulation != nullptr, ExcInternalError());
1918 
1919  const bool distributed_triangulation =
1920  dynamic_cast<
1922  &(*triangulation)) != nullptr;
1923  (void)distributed_triangulation;
1924 
1925  Assert(
1926  distributed_triangulation || number_of_locally_owned_particles == 0,
1927  ExcMessage(
1928  "Mesh refinement in a non-distributed triangulation is not supported "
1929  "by the ParticleHandler class. Either insert particles after mesh "
1930  "creation, or use a distributed triangulation."));
1931 
1932  // Resize the container if it is possible without
1933  // transferring particles
1935  particles.resize(triangulation->n_active_cells());
1936  }
1937 
1938 
1939 
1940  template <int dim, int spacedim>
1941  void
1943  {
1945  }
1946 
1947 
1948 
1949  template <int dim, int spacedim>
1950  void
1952  {
1954  }
1955 
1956 
1957 
1958  template <int dim, int spacedim>
1959  void
1961  {
1963  }
1964 
1965 
1966 
1967  template <int dim, int spacedim>
1968  void
1970  {
1972  *distributed_triangulation =
1975  *>(&(*triangulation)));
1976  (void)distributed_triangulation;
1977 
1978  Assert(
1979  distributed_triangulation != nullptr,
1980  ExcMessage(
1981  "Mesh refinement in a non-distributed triangulation is not supported "
1982  "by the ParticleHandler class. Either insert particles after mesh "
1983  "creation and do not refine afterwards, or use a distributed triangulation."));
1984 
1985 #ifdef DEAL_II_WITH_P4EST
1986  const auto callback_function =
1987  [this](
1989  & cell_iterator,
1990  const typename Triangulation<dim, spacedim>::CellStatus cell_status) {
1991  return this->pack_callback(cell_iterator, cell_status);
1992  };
1993 
1994  handle = distributed_triangulation->register_data_attach(
1995  callback_function, /*returns_variable_size_data=*/true);
1996 #endif
1997  }
1998 
1999 
2000 
2001  template <int dim, int spacedim>
2002  void
2004  {
2005  const bool serialization = false;
2006  notify_ready_to_unpack(serialization);
2007  }
2008 
2009 
2010 
2011  template <int dim, int spacedim>
2012  void
2014  {
2015  const bool serialization = true;
2016  notify_ready_to_unpack(serialization);
2017  }
2018 
2019 
2020  template <int dim, int spacedim>
2021  void
2023  const bool serialization)
2024  {
2025  notify_ready_to_unpack(serialization);
2026  }
2027 
2028 
2029 
2030  template <int dim, int spacedim>
2031  void
2033  const bool serialization)
2034  {
2036  *distributed_triangulation =
2039  *>(&(*triangulation)));
2040  (void)distributed_triangulation;
2041 
2042  Assert(
2043  distributed_triangulation != nullptr,
2044  ExcMessage(
2045  "Mesh refinement in a non-distributed triangulation is not supported "
2046  "by the ParticleHandler class. Either insert particles after mesh "
2047  "creation and do not refine afterwards, or use a distributed triangulation."));
2048 
2049  // First prepare container for insertion
2050  clear();
2051 
2052 #ifdef DEAL_II_WITH_P4EST
2053  // If we are resuming from a checkpoint, we first have to register the
2054  // store function again, to set the triangulation to the same state as
2055  // before the serialization. Only afterwards we know how to deserialize the
2056  // data correctly.
2057  if (serialization)
2059 
2060  // Check if something was stored and load it
2062  {
2063  const auto callback_function =
2064  [this](
2066  &cell_iterator,
2067  const typename Triangulation<dim, spacedim>::CellStatus cell_status,
2068  const boost::iterator_range<std::vector<char>::const_iterator>
2069  &range_iterator) {
2070  this->unpack_callback(cell_iterator, cell_status, range_iterator);
2071  };
2072 
2073  distributed_triangulation->notify_ready_to_unpack(handle,
2074  callback_function);
2075 
2076  // Reset handle and update global numbers.
2079  }
2080 #else
2081  (void)serialization;
2082 #endif
2083  }
2084 
2085 
2086 
2087  template <int dim, int spacedim>
2088  std::vector<char>
2090  const typename Triangulation<dim, spacedim>::cell_iterator &cell,
2091  const typename Triangulation<dim, spacedim>::CellStatus status) const
2092  {
2093  std::vector<particle_iterator> stored_particles_on_cell;
2094 
2095  switch (status)
2096  {
2099  // If the cell persist or is refined store all particles of the
2100  // current cell.
2101  {
2102  const unsigned int n_particles =
2103  particles[cell->active_cell_index()].size();
2104  stored_particles_on_cell.reserve(n_particles);
2105 
2106  for (unsigned int i = 0; i < n_particles; ++i)
2107  stored_particles_on_cell.push_back(
2109  }
2110  break;
2111 
2113  // If this cell is the parent of children that will be coarsened,
2114  // collect the particles of all children.
2115  {
2116  for (const auto &child : cell->child_iterators())
2117  {
2118  const unsigned int n_particles = n_particles_in_cell(child);
2119 
2120  stored_particles_on_cell.reserve(
2121  stored_particles_on_cell.size() + n_particles);
2122 
2123  for (unsigned int i = 0; i < n_particles; ++i)
2124  stored_particles_on_cell.push_back(
2126  }
2127  }
2128  break;
2129 
2130  default:
2131  Assert(false, ExcInternalError());
2132  break;
2133  }
2134 
2135  return pack_particles(stored_particles_on_cell);
2136  }
2137 
2138 
2139 
2140  template <int dim, int spacedim>
2141  void
2143  const typename Triangulation<dim, spacedim>::cell_iterator & cell,
2144  const typename Triangulation<dim, spacedim>::CellStatus status,
2145  const boost::iterator_range<std::vector<char>::const_iterator> &data_range)
2146  {
2147  const auto cell_to_store_particles =
2149  cell :
2150  cell->child(0);
2151 
2152  // deserialize particles and insert into local storage
2153  {
2154  const void *data = static_cast<const void *>(&(*data_range.begin()));
2155 
2156  while (data < &(*data_range.end()))
2157  insert_particle(data, cell_to_store_particles);
2158 
2159  Assert(
2160  data == &(*data_range.end()),
2161  ExcMessage(
2162  "The particle data could not be deserialized successfully. "
2163  "Check that when deserializing the particles you expect the same "
2164  "number of properties that were serialized."));
2165  }
2166 
2167  auto loaded_particles_on_cell = particles_in_cell(cell_to_store_particles);
2168 
2169  // now update particle storage location and properties if necessary
2170  switch (status)
2171  {
2173  {
2174  // all particles are correctly inserted
2175  }
2176  break;
2177 
2179  {
2180  // all particles are in correct cell, but their reference location
2181  // has changed
2182  for (auto &particle : loaded_particles_on_cell)
2183  {
2184  const Point<dim> p_unit =
2185  mapping->transform_real_to_unit_cell(cell_to_store_particles,
2186  particle.get_location());
2187  particle.set_reference_location(p_unit);
2188  }
2189  }
2190  break;
2191 
2193  {
2194  // we need to find the correct child to store the particles and
2195  // their reference location has changed
2196  const unsigned int stored_index =
2197  cell_to_store_particles->active_cell_index();
2198 
2199  // Cannot use range-based loop, because number of particles in cell
2200  // is going to change
2201  auto particle = loaded_particles_on_cell.begin();
2202  for (unsigned int i = 0; i < particles[stored_index].size();)
2203  {
2204  for (unsigned int child_index = 0;
2205  child_index < GeometryInfo<dim>::max_children_per_cell;
2206  ++child_index)
2207  {
2209  child = cell->child(child_index);
2210 
2211  try
2212  {
2213  const Point<dim> p_unit =
2214  mapping->transform_real_to_unit_cell(
2215  child, particle->get_location());
2217  {
2218  particle->set_reference_location(p_unit);
2219 
2220  // if the particle is not in child 0, we stored the
2221  // handle in the wrong place; move the handle and
2222  // redo the loop; otherwise move on to next particle
2223  if (child_index != 0)
2224  {
2225  particles[child->active_cell_index()].push_back(
2226  particles[stored_index][i]);
2227 
2228  particles[stored_index][i] =
2229  particles[stored_index].back();
2230  particles[stored_index].resize(
2231  particles[stored_index].size() - 1);
2232  }
2233  else
2234  {
2235  ++i;
2236  ++particle;
2237  }
2238  break;
2239  }
2240  }
2241  catch (typename Mapping<dim>::ExcTransformationFailed &)
2242  {}
2243  }
2244  }
2245  }
2246  break;
2247 
2248  default:
2249  Assert(false, ExcInternalError());
2250  break;
2251  }
2252  }
2253 } // namespace Particles
2254 
2255 #include "particle_handler.inst"
2256 
types::particle_index n_global_particles() const
IteratorRange< BaseIterator > make_iterator_range(const BaseIterator &begin, const typename identity< BaseIterator >::type &end)
static const unsigned int invalid_unsigned_int
Definition: types.h:196
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
void send_recv_particles_properties_and_location(const std::map< types::subdomain_id, std::vector< particle_iterator >> &particles_to_send, particle_container &received_particles)
types::particle_index n_locally_owned_particles() const
void get_particle_positions(VectorType &output_vector, const bool add_to_output_vector=false)
types::global_dof_index size_type
Definition: cuda_kernels.h:45
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:789
void remove_particle(const particle_iterator &particle)
PropertyPool< dim, spacedim > & get_property_pool() const
void copy_from(const ParticleHandler< dim, spacedim > &particle_handler)
void exchange_ghost_particles(const bool enable_ghost_cache=false)
types::particle_index number_of_locally_owned_particles
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)
void notify_ready_to_unpack(const bool serialization)
types::particle_index n_particles_in_cell(const typename Triangulation< dim, spacedim >::active_cell_iterator &cell) const
types::particle_index get_id() const
Definition: particle.h:546
IndexSet locally_owned_particle_ids() const
std::vector< boost::signals2::connection > tria_listeners
types::particle_index get_max_local_particle_index() const
void add_indices(const ForwardIterator &begin, const ForwardIterator &end)
Definition: index_set.h:1708
std::unique_ptr< PropertyPool< dim, spacedim > > property_pool
const Point< spacedim > & get_location() const
Definition: particle.h:519
types::particle_index next_free_particle_index
boost::signals2::signal< void(const typename Particles::ParticleIterator< dim, spacedim > &particle, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)> particle_lost
SmartPointer< const Triangulation< dim, spacedim >, ParticleHandler< dim, spacedim > > triangulation
STL namespace.
const ArrayView< double > get_properties()
Definition: particle.cc:330
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
void remove_particles(const std::vector< particle_iterator > &particles)
unsigned int global_max_particles_per_cell
std::array< unsigned int, 4 > binary_type
Definition: cell_id.h:80
particle_container particles
cell_iterator begin(const unsigned int level=0) const
Definition: tria.cc:13047
void set_property_pool(PropertyPool< dim, spacedim > &property_pool)
Definition: particle.h:564
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:110
std::size_t size() const
Definition: array_view.h:575
void initialize(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping, const unsigned int n_properties=0)
#define DEAL_II_PARTICLE_INDEX_MPI_TYPE
Definition: property_pool.h:72
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)
std::size_t serialized_size_in_bytes() const
Definition: particle.cc:287
unsigned int particle_index
Definition: property_pool.h:65
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)
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:1461
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:303
const Point< dim > & get_reference_location() const
Definition: particle.h:537
static ::ExceptionBase & ExcPointNotAvailableHere()
std::enable_if< std::is_convertible< VectorType *, Function< spacedim > * >::value==false >::type set_particle_positions(const VectorType &input_vector, const bool displace_particles=true)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
particle_iterator begin() const
ParticleHandler<dim, spacedim>::send_recv_particles.
Definition: mpi_tags.h:114
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:760
void send_recv_particles(const std::map< types::subdomain_id, std::vector< particle_iterator >> &particles_to_send, particle_container &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 >>(), const bool enable_cache=false)
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:13704
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
Definition: cell_id.h:70
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::unique_ptr< GridTools::Cache< dim, spacedim > > triangulation_cache
const types::subdomain_id artificial_subdomain_id
Definition: types.h:293
types::particle_index get_next_free_particle_index() const
boost::iterator_range< particle_iterator > particle_iterator_range
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:5516
ParticleIterator< dim, spacedim > particle_iterator
void advance(std::tuple< I1, I2 > &t, const unsigned int n)
#define AssertThrowMPI(error_code)
Definition: exceptions.h:1778
std::function< void *(const particle_iterator &, void *)> store_callback
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
const MPI_Comm & comm
std::vector< std::vector< typename PropertyPool< dim, spacedim >::Handle > > particle_container
std::vector< char > pack_callback(const typename Triangulation< dim, spacedim >::cell_iterator &cell, const typename Triangulation< dim, spacedim >::CellStatus status) const
static ::ExceptionBase & ExcNotImplemented()
Iterator points to a valid object.
SmartPointer< const Mapping< dim, spacedim >, ParticleHandler< dim, spacedim > > mapping
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
std::vector< T > all_gather(const MPI_Comm &comm, const T &object_to_send)
numbers::NumberTraits< Number >::real_type norm() const
IteratorState::IteratorStates state() const
std::map< unsigned int, T > some_to_some(const MPI_Comm &comm, const std::map< unsigned int, T > &objects_to_send)
internal::GhostParticlePartitioner< dim, spacedim > ghost_particles_cache
T max(const T &t, const MPI_Comm &mpi_communicator)
particle_iterator end() const
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:5684
void register_additional_store_load_functions(const std::function< std::size_t()> &size_callback, const std::function< void *(const particle_iterator &, void *)> &store_callback, const std::function< const void *(const particle_iterator &, const void *)> &load_callback)
static ::ExceptionBase & ExcInternalError()
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={})