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