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