particle_handler.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2017 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_particles_particle_handler_h
#define dealii_particles_particle_handler_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/function.h>
#include <deal.II/base/observer_pointer.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/fe/mapping.h>
 
#include <deal.II/grid/grid_tools_cache.h>
 
#include <deal.II/particles/particle.h>
#include <deal.II/particles/particle_iterator.h>
#include <deal.II/particles/partitioner.h>
#include <deal.II/particles/property_pool.h>
 
#include <boost/range/iterator_range.hpp>
#include <boost/serialization/map.hpp>
 
DEAL_II_NAMESPACE_OPEN
 
namespace Particles
{
  template <int dim, int spacedim = dim>
  class ParticleHandler : public EnableObserverPointer
  {
  public:
    using particle_iterator = ParticleIterator<dim, spacedim>;
 
    using particle_iterator_range = boost::iterator_range<particle_iterator>;
 
    using particle_container =
      typename ParticleAccessor<dim, spacedim>::particle_container;
 
    ParticleHandler();
 
    ParticleHandler(const Triangulation<dim, spacedim> &tria,
                    const Mapping<dim, spacedim>       &mapping,
                    const unsigned int                  n_properties = 0);
 
    virtual ~ParticleHandler();
 
    void
    initialize(const Triangulation<dim, spacedim> &tria,
               const Mapping<dim, spacedim>       &mapping,
               const unsigned int                  n_properties = 0);
 
    void
    copy_from(const ParticleHandler<dim, spacedim> &particle_handler);
 
    void
    clear();
 
    void
    clear_particles();
 
    void
    reserve(const std::size_t n_particles);
 
    void
    update_cached_numbers();
 
    particle_iterator
    begin() const;
 
    particle_iterator
    begin();
 
    particle_iterator
    end() const;
 
    particle_iterator
    end();
 
    particle_iterator
    begin_ghost() const;
 
    particle_iterator
    begin_ghost();
 
    particle_iterator
    end_ghost() const;
 
    particle_iterator
    end_ghost();
 
    types::particle_index
    n_particles_in_cell(
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
      const;
 
    particle_iterator_range
    particles_in_cell(
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell);
 
    particle_iterator_range
    particles_in_cell(
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
      const;
 
    void
    remove_particle(const particle_iterator &particle);
 
    void
    remove_particles(const std::vector<particle_iterator> &particles);
 
    particle_iterator
    insert_particle(
      const Particle<dim, spacedim> &particle,
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell);
 
    particle_iterator
    insert_particle(
      const Point<spacedim>      &position,
      const Point<dim>           &reference_position,
      const types::particle_index particle_index,
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
      const ArrayView<const double> &properties = {});
 
    void
    insert_particles(
      const std::multimap<
        typename Triangulation<dim, spacedim>::active_cell_iterator,
        Particle<dim, spacedim>> &particles);
 
    void
    insert_particles(const std::vector<Point<spacedim>> &positions);
 
    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        = {});
 
    std::map<unsigned int, IndexSet>
    insert_global_particles(
      const std::vector<Particle<dim, spacedim>> &particles,
      const std::vector<std::vector<BoundingBox<spacedim>>>
        &global_bounding_boxes);
 
    template <typename VectorType>
    std::enable_if_t<
      std::is_convertible_v<VectorType *, Function<spacedim> *> == false>
    set_particle_positions(const VectorType &input_vector,
                           const bool        displace_particles = true);
 
    void
    set_particle_positions(const std::vector<Point<spacedim>> &new_positions,
                           const bool displace_particles = true);
 
 
    void
    set_particle_positions(const Function<spacedim> &function,
                           const bool                displace_particles = true);
 
    template <typename VectorType>
    void
    get_particle_positions(VectorType &output_vector,
                           const bool  add_to_output_vector = false);
 
    void
    get_particle_positions(std::vector<Point<spacedim>> &positions,
                           const bool add_to_output_vector = false);
 
    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);
 
    types::particle_index
    n_global_particles() const;
 
    types::particle_index
    n_global_max_particles_per_cell() const;
 
    types::particle_index
    n_locally_owned_particles() const;
 
    types::particle_index
    get_next_free_particle_index() const;
 
    IndexSet
    locally_owned_particle_ids() const;
 
    unsigned int
    n_properties_per_particle() const;
 
    types::particle_index
    get_max_local_particle_index() const;
 
    PropertyPool<dim, spacedim> &
    get_property_pool() const;
 
    void
    sort_particles_into_subdomains_and_cells();
 
    void
    exchange_ghost_particles(const bool enable_ghost_cache = false);
 
    void
    update_ghost_particles();
 
    void
    prepare_for_coarsening_and_refinement();
 
    void
    unpack_after_coarsening_and_refinement();
 
    void
    prepare_for_serialization();
 
    void
    deserialize();
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);
 
    struct Signals
    {
      boost::signals2::signal<void(
        const typename Particles::ParticleIterator<dim, spacedim> &particle,
        const typename Triangulation<dim, spacedim>::active_cell_iterator
          &cell)>
        particle_lost;
    };
 
    mutable Signals signals;
 
  private:
    particle_iterator
    insert_particle(
      const void                                                       *&data,
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell);
 
    particle_iterator
    insert_particle(
      const typename PropertyPool<dim, spacedim>::Handle
        tria_attached_data_index,
      const typename Triangulation<dim, spacedim>::cell_iterator &cell);
 
    void
    reset_particle_container(particle_container &particles);
 
    ObserverPointer<const Triangulation<dim, spacedim>,
                    ParticleHandler<dim, spacedim>>
      triangulation;
 
    ObserverPointer<const Mapping<dim, spacedim>,
                    ParticleHandler<dim, spacedim>>
      mapping;
 
    std::unique_ptr<PropertyPool<dim, spacedim>> property_pool;
 
    particle_container particles;
 
    const typename particle_container::iterator owned_particles_end;
 
    std::vector<typename particle_container::iterator> cells_to_particle_cache;
 
    types::particle_index global_number_of_particles;
 
    types::particle_index number_of_locally_owned_particles;
 
    unsigned int global_max_particles_per_cell;
 
    types::particle_index next_free_particle_index;
 
    std::function<std::size_t()> size_callback;
 
    std::function<void *(const particle_iterator &, void *)> store_callback;
 
    std::function<const void *(const particle_iterator &, const void *)>
      load_callback;
 
    unsigned int tria_attached_data_index;
 
    const double tolerance_inside_cell = 1e-12;
 
    std::unique_ptr<GridTools::Cache<dim, spacedim>> triangulation_cache;
 
#ifdef DEAL_II_WITH_MPI
    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);
 
    void
    send_recv_particles_properties_and_location(
      const std::map<types::subdomain_id, std::vector<particle_iterator>>
        &particles_to_send);
 
#endif
 
    internal::GhostParticlePartitioner<dim, spacedim> ghost_particles_cache;
 
    void
    connect_to_triangulation_signals();
 
    std::vector<boost::signals2::connection> tria_listeners;
 
    void
    post_mesh_change_action();
 
    void
    register_data_attach();
 
    void
    notify_ready_to_unpack(const bool serialization);
 
    std::vector<char>
    pack_callback(
      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
      const CellStatus                                            status) 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);
 
    typename particle_container::iterator
    particle_container_owned_begin() const;
 
    typename particle_container::iterator
    particle_container_owned_end() const;
 
    typename particle_container::iterator
    particle_container_ghost_begin() const;
 
    typename particle_container::iterator
    particle_container_ghost_end() const;
  };
 
 
 
  /* ---------------------- inline and template functions ------------------
   */
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::begin() const
  {
    return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::begin()
  {
    return particle_iterator(particle_container_owned_begin(),
                             *property_pool,
                             0);
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::end() const
  {
    return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::end()
  {
    return particle_iterator(particle_container_owned_end(), *property_pool, 0);
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::begin_ghost() const
  {
    return (const_cast<ParticleHandler<dim, spacedim> *>(this))->begin_ghost();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::begin_ghost()
  {
    return particle_iterator(particle_container_ghost_begin(),
                             *property_pool,
                             0);
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::end_ghost() const
  {
    return (const_cast<ParticleHandler<dim, spacedim> *>(this))->end_ghost();
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_iterator
  ParticleHandler<dim, spacedim>::end_ghost()
  {
    return particle_iterator(particle_container_ghost_end(), *property_pool, 0);
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_container::iterator
  ParticleHandler<dim, spacedim>::particle_container_owned_begin() const
  {
    // We should always have at least the three anchor entries in the list of
    // particles
    Assert(!particles.empty(), ExcInternalError());
    typename particle_container::iterator begin =
      const_cast<particle_container &>(particles).begin();
    return ++begin;
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_container::iterator
  ParticleHandler<dim, spacedim>::particle_container_owned_end() const
  {
    Assert(!particles.empty(), ExcInternalError());
    return owned_particles_end;
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_container::iterator
  ParticleHandler<dim, spacedim>::particle_container_ghost_begin() const
  {
    Assert(!particles.empty(), ExcInternalError());
    typename particle_container::iterator begin = owned_particles_end;
    return ++begin;
  }
 
 
 
  template <int dim, int spacedim>
  inline typename ParticleHandler<dim, spacedim>::particle_container::iterator
  ParticleHandler<dim, spacedim>::particle_container_ghost_end() const
  {
    Assert(!particles.empty(), ExcInternalError());
    typename particle_container::iterator end =
      const_cast<particle_container &>(particles).end();
    return --end;
  }
 
 
 
  template <int dim, int spacedim>
  template <class Archive>
  inline void
  ParticleHandler<dim, spacedim>::serialize(Archive &ar, const unsigned int)
  {
    // Note that we do not serialize the particle data itself (i.e., the
    // 'particles' member variable). Instead we
    // use the serialization functionality of the triangulation class, because
    // this guarantees that data is immediately shipped to new processes if
    // the domain is distributed differently after resuming from a checkpoint.
    //
    // See @ref step_83 "step-83" for how to serialize ParticleHandler objects.
    ar &global_number_of_particles &global_max_particles_per_cell
      &next_free_particle_index;
  }
 
 
 
  template <int dim, int spacedim>
  template <typename VectorType>
  inline std::enable_if_t<
    std::is_convertible_v<VectorType *, Function<spacedim> *> == false>
  ParticleHandler<dim, spacedim>::set_particle_positions(
    const VectorType &input_vector,
    const bool        displace_particles)
  {
    AssertDimension(input_vector.size(),
                    get_next_free_particle_index() * spacedim);
    for (auto &p : *this)
      {
        Point<spacedim> &position = p.get_location();
        const auto       id       = p.get_id();
 
        if (displace_particles)
          for (unsigned int i = 0; i < spacedim; ++i)
            position[i] += input_vector[id * spacedim + i];
        else
          for (unsigned int i = 0; i < spacedim; ++i)
            position[i] = input_vector[id * spacedim + i];
      }
    sort_particles_into_subdomains_and_cells();
  }
 
 
 
  template <int dim, int spacedim>
  template <typename VectorType>
  inline void
  ParticleHandler<dim, spacedim>::get_particle_positions(
    VectorType &output_vector,
    const bool  add_to_output_vector)
  {
    AssertDimension(output_vector.size(),
                    get_next_free_particle_index() * spacedim);
    for (const auto &p : *this)
      {
        auto       point = p.get_location();
        const auto id    = p.get_id();
        if (add_to_output_vector)
          for (unsigned int i = 0; i < spacedim; ++i)
            output_vector[id * spacedim + i] += point[i];
        else
          for (unsigned int i = 0; i < spacedim; ++i)
            output_vector[id * spacedim + i] = point[i];
      }
    if (add_to_output_vector)
      output_vector.compress(VectorOperation::add);
    else
      output_vector.compress(VectorOperation::insert);
  }
 
} // namespace Particles
 
DEAL_II_NAMESPACE_CLOSE
 
#endif