Loading [MathJax]/extensions/TeX/newcommand.js
 deal.II version GIT relicensing-3075-gc235bd4825 2025-04-15 08:10: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_accessor.h
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
15#ifndef dealii_particles_particle_accessor_h
16#define dealii_particles_particle_accessor_h
17
18#include <deal.II/base/config.h>
19
21
22#include <deal.II/grid/tria.h>
24
27
28#include <boost/geometry/index/indexable.hpp>
29#include <boost/serialization/array_wrapper.hpp>
30
31#include <list>
32
33
35
36namespace Particles
37{
38 // Forward declarations
39#ifndef DOXYGEN
40 template <int, int>
41 class ParticleIterator;
42 template <int, int>
43 class ParticleHandler;
44#endif
45
49 template <int dim, int spacedim = dim>
51 {
52 public:
87 {
91 ParticlesInCell() = default;
92
103
107 std::vector<typename PropertyPool<dim, spacedim>::Handle> particles;
108
113 };
114
118 using particle_container = std::list<ParticlesInCell>;
119
123 void *
125
126
130 const void *
132
153 void
154 set_location(const Point<spacedim> &new_location);
155
161 const Point<spacedim> &
163
187
208 void
209 set_reference_location(const Point<dim> &new_reference_location);
210
214 const Point<dim> &
216
220 void
222
227 get_id() const;
228
249
254 bool
256
277 void
278 set_properties(const std::vector<double> &new_properties);
279
300 void
302
326 void
327 set_properties(const Tensor<1, dim> &new_properties);
328
336 {
337 // The implementation is up here inside the class declaration because
338 // NVCC (at least in 12.5 and 12.6) otherwise produce a compile error:
339 //
340 // error: no declaration matches ‘::ArrayView<__remove_cv(const
341 // double)> ::Particles::ParticleAccessor<dim,
342 // spacedim>::get_properties()’
343 //
344 // See https://github.com/dealii/dealii/issues/17148
346
348 }
349
350
358
364 std::size_t
366
372
378 template <class Archive>
379 void
380 save(Archive &ar, const unsigned int version) const;
381
387 template <class Archive>
388 void
389 load(Archive &ar, const unsigned int version);
390
391#ifdef DOXYGEN
397 template <class Archive>
398 void
399 serialize(Archive &archive, const unsigned int version);
400#else
401 // This macro defines the serialize() method that is compatible with
402 // the templated save() and load() method that have been implemented.
403 BOOST_SERIALIZATION_SPLIT_MEMBER()
404#endif
405
409 void
411
415 void
417
421 bool
423
427 bool
429
434 state() const;
435
436 private:
441
449 const typename particle_container::iterator particles_in_cell,
451 const unsigned int particle_index_within_cell);
452
459 get_handle() const;
460
466
471 typename particle_container::iterator particles_in_cell;
472
477
482
483 // Make ParticleIterator a friend to allow it constructing
484 // ParticleAccessors.
485 template <int, int>
486 friend class ParticleIterator;
487 template <int, int>
488 friend class ParticleHandler;
489 };
490
491
492
493 template <int dim, int spacedim>
494 template <class Archive>
495 inline void
496 ParticleAccessor<dim, spacedim>::load(Archive &ar, const unsigned int)
497 {
498 unsigned int n_properties = 0;
499
500 Point<spacedim> location;
501 Point<dim> reference_location;
503 ar &location &reference_location &id &n_properties;
504
505 set_location(location);
506 set_reference_location(reference_location);
507 set_id(id);
508
509 if (n_properties > 0)
510 {
511 ArrayView<double> properties(get_properties());
512 Assert(
513 properties.size() == n_properties,
515 "This particle was serialized with " +
516 std::to_string(n_properties) +
517 " properties, but the new property handler provides space for " +
518 std::to_string(properties.size()) +
519 " properties. Deserializing a particle only works for matching property sizes."));
520
521 ar &boost::serialization::make_array(properties.data(), n_properties);
522 }
523 }
524
525
526
527 template <int dim, int spacedim>
528 template <class Archive>
529 inline void
530 ParticleAccessor<dim, spacedim>::save(Archive &ar, const unsigned int) const
531 {
532 unsigned int n_properties = 0;
533 if ((property_pool != nullptr) &&
535 n_properties = get_properties().size();
536
537 Point<spacedim> location = get_location();
538 Point<dim> reference_location = get_reference_location();
539 types::particle_index id = get_id();
540
541 ar &location &reference_location &id &n_properties;
542
543 if (n_properties > 0)
544 ar &boost::serialization::make_array(get_properties().data(),
545 n_properties);
546 }
547
548
549 // ------------------------- inline functions ------------------------------
550
551 template <int dim, int spacedim>
553 : particles_in_cell(typename particle_container::iterator())
554 , property_pool(nullptr)
555 , particle_index_within_cell(numbers::invalid_unsigned_int)
556 {}
557
558
559
560 template <int dim, int spacedim>
562 const typename particle_container::iterator particles_in_cell,
563 const PropertyPool<dim, spacedim> &property_pool,
564 const unsigned int particle_index_within_cell)
565 : particles_in_cell(particles_in_cell)
566 , property_pool(const_cast<PropertyPool<dim, spacedim> *>(&property_pool))
567 , particle_index_within_cell(particle_index_within_cell)
568 {}
569
570
571
572 template <int dim, int spacedim>
573 inline const void *
575 const void *data)
576 {
578
579 const types::particle_index *id_data =
580 static_cast<const types::particle_index *>(data);
581 set_id(*id_data++);
582 const double *pdata = reinterpret_cast<const double *>(id_data);
583
584 Point<spacedim> location;
585 for (unsigned int i = 0; i < spacedim; ++i)
586 location[i] = *pdata++;
587 set_location(location);
588
589 Point<dim> reference_location;
590 for (unsigned int i = 0; i < dim; ++i)
591 reference_location[i] = *pdata++;
592 set_reference_location(reference_location);
593
594 // See if there are properties to load
595 if (has_properties())
596 {
597 const ArrayView<double> particle_properties =
598 property_pool->get_properties(get_handle());
599 const unsigned int size = particle_properties.size();
600 for (unsigned int i = 0; i < size; ++i)
601 particle_properties[i] = *pdata++;
602 }
603
604 return static_cast<const void *>(pdata);
605 }
606
607
608
609 template <int dim, int spacedim>
610 inline void *
612 void *data) const
613 {
615
616 types::particle_index *id_data = static_cast<types::particle_index *>(data);
617 *id_data = get_id();
618 ++id_data;
619 double *pdata = reinterpret_cast<double *>(id_data);
620
621 // Write location
622 for (unsigned int i = 0; i < spacedim; ++i, ++pdata)
623 *pdata = get_location()[i];
624
625 // Write reference location
626 for (unsigned int i = 0; i < dim; ++i, ++pdata)
627 *pdata = get_reference_location()[i];
628
629 // Write properties
630 if (has_properties())
631 {
632 const ArrayView<double> particle_properties =
633 property_pool->get_properties(get_handle());
634 for (unsigned int i = 0; i < particle_properties.size(); ++i, ++pdata)
635 *pdata = particle_properties[i];
636 }
637
638 return static_cast<void *>(pdata);
639 }
640
641
642
643 template <int dim, int spacedim>
644 inline void
646 {
648
649 property_pool->set_location(get_handle(), new_loc);
650 }
651
652
653
654 template <int dim, int spacedim>
655 inline const Point<spacedim> &
657 {
659
660 return property_pool->get_location(get_handle());
661 }
662
663
664
665 template <int dim, int spacedim>
666 inline Point<spacedim> &
668 {
670
671 return property_pool->get_location(get_handle());
672 }
673
674
675
676 template <int dim, int spacedim>
677 inline void
679 const Point<dim> &new_loc)
680 {
682
683 property_pool->set_reference_location(get_handle(), new_loc);
684 }
685
686
687
688 template <int dim, int spacedim>
689 inline const Point<dim> &
691 {
693
694 return property_pool->get_reference_location(get_handle());
695 }
696
697
698
699 template <int dim, int spacedim>
702 {
704
705 return property_pool->get_id(get_handle());
706 }
707
708
709
710 template <int dim, int spacedim>
713 {
715
716 return get_handle();
717 }
718
719
720
721 template <int dim, int spacedim>
722 inline void
724 {
726
727 property_pool->set_id(get_handle(), new_id);
728 }
729
730
731
732 template <int dim, int spacedim>
733 inline bool
735 {
737
738 // Particles always have a property pool associated with them,
739 // but we can access properties only if there is a valid handle.
740 // The only way a particle can have no valid handle if it has
741 // been moved-from -- but that leaves an object in an invalid
742 // state, and so we can just assert that that can't be the case.
745 return (property_pool->n_properties_per_slot() > 0);
746 }
747
748
749
750 template <int dim, int spacedim>
751 inline void
753 const std::vector<double> &new_properties)
754 {
756
757 set_properties(
758 ArrayView<const double>(new_properties.data(), new_properties.size()));
759 }
760
761
762
763 template <int dim, int spacedim>
764 inline void
766 const ArrayView<const double> &new_properties)
767 {
769
770 const ArrayView<double> property_values =
771 property_pool->get_properties(get_handle());
772
773 Assert(new_properties.size() == property_values.size(),
775 "You are trying to assign properties with an incompatible length. "
776 "The particle has space to store " +
777 std::to_string(property_values.size()) +
778 " properties, but you are trying to assign " +
779 std::to_string(new_properties.size()) +
780 " properties. This is not allowed."));
781
782 if (property_values.size() > 0)
783 std::copy(new_properties.begin(),
784 new_properties.end(),
785 property_values.begin());
786 }
787
788
789
790 template <int dim, int spacedim>
791 inline void
793 const Tensor<1, dim> &new_properties)
794 {
796
797 // A Tensor object is not an array, so we cannot just create an
798 // ArrayView object for it. Rather, copy the data into a true
799 // array and make the ArrayView from that.
800 double array[dim];
801 for (unsigned int d = 0; d < dim; ++d)
802 array[d] = new_properties[d];
803
804 set_properties(make_array_view(array));
805 }
806
807
808
809 template <int dim, int spacedim>
812 {
814
815 return property_pool->get_properties(get_handle());
816 }
817
818
819
820 template <int dim, int spacedim>
821 inline const typename Triangulation<dim, spacedim>::cell_iterator &
823 {
825 Assert(particles_in_cell->cell.state() == IteratorState::valid,
827
828 return particles_in_cell->cell;
829 }
830
831
832
833 // template <int dim, int spacedim>
834 // inline ArrayView<double>
835 // ParticleAccessor<dim, spacedim>::get_properties()
836
837
838
839 template <int dim, int spacedim>
840 inline std::size_t
842 {
844
845 std::size_t size = sizeof(get_id()) +
846 sizeof(double) * spacedim + // get_location()
847 sizeof(double) * dim; // get_reference_location()
848
849 if (has_properties())
850 {
851 size += sizeof(double) * get_properties().size();
852 }
853 return size;
854 }
855
856
857
858 template <int dim, int spacedim>
859 inline void
861 {
863
864 ++particle_index_within_cell;
865
866 if (particle_index_within_cell >= particles_in_cell->particles.size())
867 {
868 particle_index_within_cell = 0;
869 ++particles_in_cell;
870 }
871 }
872
873
874
875 template <int dim, int spacedim>
876 inline void
878 {
880
881 if (particle_index_within_cell > 0)
882 --particle_index_within_cell;
883 else
884 {
885 --particles_in_cell;
886 particle_index_within_cell = particles_in_cell->particles.empty() ?
887 0 :
888 particles_in_cell->particles.size() - 1;
889 }
890 }
891
892
893
894 template <int dim, int spacedim>
895 inline bool
897 const ParticleAccessor<dim, spacedim> &other) const
898 {
899 return !(*this == other);
900 }
901
902
903
904 template <int dim, int spacedim>
905 inline bool
907 const ParticleAccessor<dim, spacedim> &other) const
908 {
909 return (property_pool == other.property_pool) &&
910 (particles_in_cell == other.particles_in_cell) &&
911 (particle_index_within_cell == other.particle_index_within_cell);
912 }
913
914
915
916 template <int dim, int spacedim>
919 {
920 if (property_pool != nullptr &&
921 particles_in_cell->cell.state() == IteratorState::valid &&
922 particle_index_within_cell < particles_in_cell->particles.size())
924 else if (property_pool != nullptr &&
925 particles_in_cell->cell.state() == IteratorState::past_the_end &&
926 particle_index_within_cell == 0)
928 else
930
932 }
933
934
935
936 template <int dim, int spacedim>
939 {
940 return particles_in_cell->particles[particle_index_within_cell];
941 }
942
943
944
945 template <int dim, int spacedim>
946 inline const typename PropertyPool<dim, spacedim>::Handle &
948 {
949 return particles_in_cell->particles[particle_index_within_cell];
950 }
951
952} // namespace Particles
953
955
956namespace boost
957{
958 namespace geometry
959 {
960 namespace index
961 {
966 template <int dim, int spacedim>
967 struct indexable<::Particles::ParticleAccessor<dim, spacedim>>
968 {
973 using result_type = const ::Point<spacedim> &;
974
976 operator()(const ::Particles::ParticleAccessor<dim, spacedim>
977 &accessor) const
978 {
979 return accessor.get_location();
980 }
981 };
982 } // namespace index
983 } // namespace geometry
984} // namespace boost
985
986#endif
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:954
iterator begin() const
Definition array_view.h:707
iterator end() const
Definition array_view.h:716
value_type * data() const noexcept
Definition array_view.h:666
std::size_t size() const
Definition array_view.h:689
particle_container::iterator particles_in_cell
const Point< spacedim > & get_location() const
types::particle_index get_id() const
void save(Archive &ar, const unsigned int version) const
types::particle_index get_local_index() const
std::list< ParticlesInCell > particle_container
const void * read_particle_data_from_memory(const void *data)
void load(Archive &ar, const unsigned int version)
const Triangulation< dim, spacedim >::cell_iterator & get_surrounding_cell() const
IteratorState::IteratorStates state() const
bool operator!=(const ParticleAccessor< dim, spacedim > &other) const
void set_properties(const Tensor< 1, dim > &new_properties)
void set_location(const Point< spacedim > &new_location)
PropertyPool< dim, spacedim >::Handle & get_handle()
void serialize(Archive &archive, const unsigned int version)
ArrayView< double > get_properties()
void set_reference_location(const Point< dim > &new_reference_location)
const Point< dim > & get_reference_location() const
std::size_t serialized_size_in_bytes() const
const PropertyPool< dim, spacedim >::Handle & get_handle() const
void set_properties(const std::vector< double > &new_properties)
Point< spacedim > & get_location()
ParticleAccessor(const typename particle_container::iterator particles_in_cell, const PropertyPool< dim, spacedim > &property_pool, const unsigned int particle_index_within_cell)
PropertyPool< dim, spacedim > * property_pool
void * write_particle_data_to_memory(void *data) const
void set_properties(const ArrayView< const double > &new_properties)
ArrayView< const double > get_properties() const
bool operator==(const ParticleAccessor< dim, spacedim > &other) const
void set_id(const types::particle_index &new_id)
ArrayView< double, ::MemorySpace::Host > get_properties(const Handle handle)
Definition point.h:113
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:35
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:36
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
std::vector< index_type > data
Definition mpi.cc:746
std::size_t size
Definition mpi.cc:745
@ past_the_end
Iterator reached end of container.
@ valid
Iterator points to a valid object.
@ invalid
Iterator is invalid, probably due to an error.
Triangulation< dim, spacedim >::active_cell_iterator cell
ParticlesInCell(const std::vector< typename PropertyPool< dim, spacedim >::Handle > &particles, const typename Triangulation< dim, spacedim >::active_cell_iterator &cell)
std::vector< typename PropertyPool< dim, spacedim >::Handle > particles
result_type operator()(const ::Particles::ParticleAccessor< dim, spacedim > &accessor) const