particle_accessor.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2017 - 2025 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_particles_particle_accessor_h
#define dealii_particles_particle_accessor_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
 
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_iterator_base.h>
 
#include <deal.II/particles/particle.h>
#include <deal.II/particles/property_pool.h>
 
#include <boost/geometry/index/indexable.hpp>
#include <boost/serialization/array_wrapper.hpp>
 
#include <list>
 
 
DEAL_II_NAMESPACE_OPEN
 
namespace Particles
{
  // Forward declarations
#ifndef DOXYGEN
  template <int, int>
  class ParticleIterator;
  template <int, int>
  class ParticleHandler;
#endif

  template <int dim, int spacedim = dim>
  class ParticleAccessor
  {
  public:
    struct ParticlesInCell
    {
      ParticlesInCell() = default;

      ParticlesInCell(
        const std::vector<typename PropertyPool<dim, spacedim>::Handle>
          &particles,
        const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
        : particles(particles)
        , cell(cell)
      {}

      std::vector<typename PropertyPool<dim, spacedim>::Handle> particles;

      typename Triangulation<dim, spacedim>::active_cell_iterator cell;
    };

    using particle_container = std::list<ParticlesInCell>;

    void *
    write_particle_data_to_memory(void *data) const;
 

    const void *
    read_particle_data_from_memory(const void *data);

    void
    set_location(const Point<spacedim> &new_location);

    const Point<spacedim> &
    get_location() const;

    Point<spacedim> &
    get_location();

    void
    set_reference_location(const Point<dim> &new_reference_location);

    const Point<dim> &
    get_reference_location() const;

    void
    set_id(const types::particle_index &new_id);

    types::particle_index
    get_id() const;

    types::particle_index
    get_local_index() const;

    bool
    has_properties() const;

    void
    set_properties(const std::vector<double> &new_properties);

    void
    set_properties(const ArrayView<const double> &new_properties);

    void
    set_properties(const Tensor<1, dim> &new_properties);

    ArrayView<double>
    get_properties()
    {
      // The implementation is up here inside the class declaration because
      // NVCC (at least in 12.5 and 12.6) otherwise produce a compile error:
      //
      // error: no declaration matches ‘::ArrayView<__remove_cv(const
      // double)> ::Particles::ParticleAccessor<dim,
      // spacedim>::get_properties()’
      //
      // See https://github.com/dealii/dealii/issues/17148
      Assert(state() == IteratorState::valid, ExcInternalError());
 
      return property_pool->get_properties(get_handle());
    }
 

    ArrayView<const double>
    get_properties() const;

    std::size_t
    serialized_size_in_bytes() const;

    const typename Triangulation<dim, spacedim>::cell_iterator &
    get_surrounding_cell() const;

    template <class Archive>
    void
    save(Archive &ar, const unsigned int version) const;

    template <class Archive>
    void
    load(Archive &ar, const unsigned int version);
 
#ifdef DOXYGEN
    template <class Archive>
    void
    serialize(Archive &archive, const unsigned int version);
#else
    // This macro defines the serialize() method that is compatible with
    // the templated save() and load() method that have been implemented.
    BOOST_SERIALIZATION_SPLIT_MEMBER()
#endif

    void
    next();

    void
    prev();

    bool
    operator!=(const ParticleAccessor<dim, spacedim> &other) const;

    bool
    operator==(const ParticleAccessor<dim, spacedim> &other) const;

    IteratorState::IteratorStates
    state() const;
 
  private:
    ParticleAccessor();

    ParticleAccessor(
      const typename particle_container::iterator particles_in_cell,
      const PropertyPool<dim, spacedim>          &property_pool,
      const unsigned int                          particle_index_within_cell);

    const typename PropertyPool<dim, spacedim>::Handle &
    get_handle() const;

    typename PropertyPool<dim, spacedim>::Handle &
    get_handle();

    typename particle_container::iterator particles_in_cell;

    PropertyPool<dim, spacedim> *property_pool;

    unsigned int particle_index_within_cell;
 
    // Make ParticleIterator a friend to allow it constructing
    // ParticleAccessors.
    template <int, int>
    friend class ParticleIterator;
    template <int, int>
    friend class ParticleHandler;
  };
 
 
 
  template <int dim, int spacedim>
  template <class Archive>
  inline void
  ParticleAccessor<dim, spacedim>::load(Archive &ar, const unsigned int)
  {
    unsigned int n_properties = 0;
 
    Point<spacedim>                       location;
    Point<dim>                            reference_location;
    types::particle_index                 id;
    ar &location &reference_location &id &n_properties;
 
    set_location(location);
    set_reference_location(reference_location);
    set_id(id);
 
    if (n_properties > 0)
      {
        ArrayView<double> properties(get_properties());
        Assert(
          properties.size() == n_properties,
          ExcMessage(
            "This particle was serialized with " +
            std::to_string(n_properties) +
            " properties, but the new property handler provides space for " +
            std::to_string(properties.size()) +
            " properties. Deserializing a particle only works for matching property sizes."));
 
        ar &boost::serialization::make_array(properties.data(), n_properties);
      }
  }
 
 
 
  template <int dim, int spacedim>
  template <class Archive>
  inline void
  ParticleAccessor<dim, spacedim>::save(Archive &ar, const unsigned int) const
  {
    unsigned int n_properties = 0;
    if ((property_pool != nullptr) &&
        (get_handle() != PropertyPool<dim, spacedim>::invalid_handle))
      n_properties = get_properties().size();
 
    Point<spacedim>       location           = get_location();
    Point<dim>            reference_location = get_reference_location();
    types::particle_index id                 = get_id();
 
    ar &location &reference_location &id &n_properties;
 
    if (n_properties > 0)
      ar &boost::serialization::make_array(get_properties().data(),
                                           n_properties);
  }
 
 
  // ------------------------- inline functions ------------------------------
 
  template <int dim, int spacedim>
  inline ParticleAccessor<dim, spacedim>::ParticleAccessor()
    : particles_in_cell(typename particle_container::iterator())
    , property_pool(nullptr)
    , particle_index_within_cell(numbers::invalid_unsigned_int)
  {}
 
 
 
  template <int dim, int spacedim>
  inline ParticleAccessor<dim, spacedim>::ParticleAccessor(
    const typename particle_container::iterator particles_in_cell,
    const PropertyPool<dim, spacedim>          &property_pool,
    const unsigned int                          particle_index_within_cell)
    : particles_in_cell(particles_in_cell)
    , property_pool(const_cast<PropertyPool<dim, spacedim> *>(&property_pool))
    , particle_index_within_cell(particle_index_within_cell)
  {}
 
 
 
  template <int dim, int spacedim>
  inline const void *
  ParticleAccessor<dim, spacedim>::read_particle_data_from_memory(
    const void *data)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    const types::particle_index *id_data =
      static_cast<const types::particle_index *>(data);
    set_id(*id_data++);
    const double *pdata = reinterpret_cast<const double *>(id_data);
 
    Point<spacedim> location;
    for (unsigned int i = 0; i < spacedim; ++i)
      location[i] = *pdata++;
    set_location(location);
 
    Point<dim> reference_location;
    for (unsigned int i = 0; i < dim; ++i)
      reference_location[i] = *pdata++;
    set_reference_location(reference_location);
 
    // See if there are properties to load
    if (has_properties())
      {
        const ArrayView<double> particle_properties =
          property_pool->get_properties(get_handle());
        const unsigned int size = particle_properties.size();
        for (unsigned int i = 0; i < size; ++i)
          particle_properties[i] = *pdata++;
      }
 
    return static_cast<const void *>(pdata);
  }
 
 
 
  template <int dim, int spacedim>
  inline void *
  ParticleAccessor<dim, spacedim>::write_particle_data_to_memory(
    void *data) const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    types::particle_index *id_data = static_cast<types::particle_index *>(data);
    *id_data                       = get_id();
    ++id_data;
    double *pdata = reinterpret_cast<double *>(id_data);
 
    // Write location
    for (unsigned int i = 0; i < spacedim; ++i, ++pdata)
      *pdata = get_location()[i];
 
    // Write reference location
    for (unsigned int i = 0; i < dim; ++i, ++pdata)
      *pdata = get_reference_location()[i];
 
    // Write properties
    if (has_properties())
      {
        const ArrayView<double> particle_properties =
          property_pool->get_properties(get_handle());
        for (unsigned int i = 0; i < particle_properties.size(); ++i, ++pdata)
          *pdata = particle_properties[i];
      }
 
    return static_cast<void *>(pdata);
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_location(const Point<spacedim> &new_loc)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    property_pool->set_location(get_handle(), new_loc);
  }
 
 
 
  template <int dim, int spacedim>
  inline const Point<spacedim> &
  ParticleAccessor<dim, spacedim>::get_location() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return property_pool->get_location(get_handle());
  }
 
 
 
  template <int dim, int spacedim>
  inline Point<spacedim> &
  ParticleAccessor<dim, spacedim>::get_location()
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return property_pool->get_location(get_handle());
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_reference_location(
    const Point<dim> &new_loc)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    property_pool->set_reference_location(get_handle(), new_loc);
  }
 
 
 
  template <int dim, int spacedim>
  inline const Point<dim> &
  ParticleAccessor<dim, spacedim>::get_reference_location() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return property_pool->get_reference_location(get_handle());
  }
 
 
 
  template <int dim, int spacedim>
  inline types::particle_index
  ParticleAccessor<dim, spacedim>::get_id() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return property_pool->get_id(get_handle());
  }
 
 
 
  template <int dim, int spacedim>
  inline types::particle_index
  ParticleAccessor<dim, spacedim>::get_local_index() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return get_handle();
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_id(const types::particle_index &new_id)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    property_pool->set_id(get_handle(), new_id);
  }
 
 
 
  template <int dim, int spacedim>
  inline bool
  ParticleAccessor<dim, spacedim>::has_properties() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    // Particles always have a property pool associated with them,
    // but we can access properties only if there is a valid handle.
    // The only way a particle can have no valid handle if it has
    // been moved-from -- but that leaves an object in an invalid
    // state, and so we can just assert that that can't be the case.
    Assert((get_handle() != PropertyPool<dim, spacedim>::invalid_handle),
           ExcInternalError());
    return (property_pool->n_properties_per_slot() > 0);
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_properties(
    const std::vector<double> &new_properties)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    set_properties(
      ArrayView<const double>(new_properties.data(), new_properties.size()));
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_properties(
    const ArrayView<const double> &new_properties)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    const ArrayView<double> property_values =
      property_pool->get_properties(get_handle());
 
    Assert(new_properties.size() == property_values.size(),
           ExcMessage(
             "You are trying to assign properties with an incompatible length. "
             "The particle has space to store " +
             std::to_string(property_values.size()) +
             " properties, but you are trying to assign " +
             std::to_string(new_properties.size()) +
             " properties. This is not allowed."));
 
    if (property_values.size() > 0)
      std::copy(new_properties.begin(),
                new_properties.end(),
                property_values.begin());
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::set_properties(
    const Tensor<1, dim> &new_properties)
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    // A Tensor object is not an array, so we cannot just create an
    // ArrayView object for it. Rather, copy the data into a true
    // array and make the ArrayView from that.
    double array[dim];
    for (unsigned int d = 0; d < dim; ++d)
      array[d] = new_properties[d];
 
    set_properties(make_array_view(array));
  }
 
 
 
  template <int dim, int spacedim>
  inline ArrayView<const double>
  ParticleAccessor<dim, spacedim>::get_properties() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    return property_pool->get_properties(get_handle());
  }
 
 
 
  template <int dim, int spacedim>
  inline const typename Triangulation<dim, spacedim>::cell_iterator &
  ParticleAccessor<dim, spacedim>::get_surrounding_cell() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
    Assert(particles_in_cell->cell.state() == IteratorState::valid,
           ExcInternalError());
 
    return particles_in_cell->cell;
  }
 
 
 
  // template <int dim, int spacedim>
  // inline ArrayView<double>
  //   ParticleAccessor<dim, spacedim>::get_properties()
 
 
 
  template <int dim, int spacedim>
  inline std::size_t
  ParticleAccessor<dim, spacedim>::serialized_size_in_bytes() const
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    std::size_t size = sizeof(get_id()) +
                       sizeof(double) * spacedim + // get_location()
                       sizeof(double) * dim;       // get_reference_location()
 
    if (has_properties())
      {
        size += sizeof(double) * get_properties().size();
      }
    return size;
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::next()
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    ++particle_index_within_cell;
 
    if (particle_index_within_cell >= particles_in_cell->particles.size())
      {
        particle_index_within_cell = 0;
        ++particles_in_cell;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline void
  ParticleAccessor<dim, spacedim>::prev()
  {
    Assert(state() == IteratorState::valid, ExcInternalError());
 
    if (particle_index_within_cell > 0)
      --particle_index_within_cell;
    else
      {
        --particles_in_cell;
        particle_index_within_cell = particles_in_cell->particles.empty() ?
                                       0 :
                                       particles_in_cell->particles.size() - 1;
      }
  }
 
 
 
  template <int dim, int spacedim>
  inline bool
  ParticleAccessor<dim, spacedim>::operator!=(
    const ParticleAccessor<dim, spacedim> &other) const
  {
    return !(*this == other);
  }
 
 
 
  template <int dim, int spacedim>
  inline bool
  ParticleAccessor<dim, spacedim>::operator==(
    const ParticleAccessor<dim, spacedim> &other) const
  {
    return (property_pool == other.property_pool) &&
           (particles_in_cell == other.particles_in_cell) &&
           (particle_index_within_cell == other.particle_index_within_cell);
  }
 
 
 
  template <int dim, int spacedim>
  inline IteratorState::IteratorStates
  ParticleAccessor<dim, spacedim>::state() const
  {
    if (property_pool != nullptr &&
        particles_in_cell->cell.state() == IteratorState::valid &&
        particle_index_within_cell < particles_in_cell->particles.size())
      return IteratorState::valid;
    else if (property_pool != nullptr &&
             particles_in_cell->cell.state() == IteratorState::past_the_end &&
             particle_index_within_cell == 0)
      return IteratorState::past_the_end;
    else
      return IteratorState::invalid;
 
    return IteratorState::invalid;
  }
 
 
 
  template <int dim, int spacedim>
  inline typename PropertyPool<dim, spacedim>::Handle &
  ParticleAccessor<dim, spacedim>::get_handle()
  {
    return particles_in_cell->particles[particle_index_within_cell];
  }
 
 
 
  template <int dim, int spacedim>
  inline const typename PropertyPool<dim, spacedim>::Handle &
  ParticleAccessor<dim, spacedim>::get_handle() const
  {
    return particles_in_cell->particles[particle_index_within_cell];
  }
 
} // namespace Particles
 
DEAL_II_NAMESPACE_CLOSE
 
namespace boost
{
  namespace geometry
  {
    namespace index
    {
      template <int dim, int spacedim>
      struct indexable<::Particles::ParticleAccessor<dim, spacedim>>
      {
        using result_type = const ::Point<spacedim> &;
 
        result_type
        operator()(const ::Particles::ParticleAccessor<dim, spacedim>
                     &accessor) const
        {
          return accessor.get_location();
        }
      };
    } // namespace index
  }   // namespace geometry
} // namespace boost
 
#endif