particle_handler.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2020 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE.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/function.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/subscriptor.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/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 Subscriptor
  {
  public:
    using particle_iterator = ParticleIterator<dim, spacedim>;
 
    using particle_iterator_range = boost::iterator_range<particle_iterator>;
 
    ParticleHandler();
 
    ParticleHandler(const Triangulation<dim, spacedim> &tria,
                    const Mapping<dim, spacedim> &      mapping,
                    const unsigned int                  n_properties = 0);
 
    virtual ~ParticleHandler() override = default;
 
    void
    initialize(const Triangulation<dim, spacedim> &tria,
               const Mapping<dim, spacedim> &      mapping,
               const unsigned int                  n_properties = 0);
 
    void
    clear();
 
    void
    clear_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();
 
    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);
 
    particle_iterator
    insert_particle(
      const Particle<dim, spacedim> &particle,
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell);
 
    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 = {});
 
    template <class VectorType>
    typename std::enable_if<
      std::is_convertible<VectorType *, Function<spacedim> *>::value ==
      false>::type
    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 <class 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_relevant_ids() const;
 
    unsigned int
    n_properties_per_particle() const;
 
    PropertyPool &
    get_property_pool() const;
 
    unsigned int
    n_particles_in_cell(
      const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
      const;
 
    void
    sort_particles_into_subdomains_and_cells();
 
    void
    exchange_ghost_particles();
 
    void
    register_store_callback_function();
 
    void
    register_load_callback_function(const bool serialization);
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);
 
  private:
    SmartPointer<const Triangulation<dim, spacedim>,
                 ParticleHandler<dim, spacedim>>
      triangulation;
 
    SmartPointer<const Mapping<dim, spacedim>, ParticleHandler<dim, spacedim>>
      mapping;
 
    std::multimap<internal::LevelInd, Particle<dim, spacedim>> particles;
 
    std::multimap<internal::LevelInd, Particle<dim, spacedim>> ghost_particles;
 
    types::particle_index global_number_of_particles;
 
    unsigned int global_max_particles_per_cell;
 
    types::particle_index next_free_particle_index;
 
    std::unique_ptr<PropertyPool> property_pool;
 
    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 handle;
 
    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,
      std::multimap<internal::LevelInd, Particle<dim, spacedim>>
        &received_particles,
      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>>());
#endif
 
    std::vector<char>
    store_particles(
      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
      const typename Triangulation<dim, spacedim>::CellStatus     status) const;
 
    void
    load_particles(
      const typename Triangulation<dim, spacedim>::cell_iterator &cell,
      const typename Triangulation<dim, spacedim>::CellStatus     status,
      const boost::iterator_range<std::vector<char>::const_iterator>
        &data_range);
  };
 
 
 
  /* ---------------------- inline and template functions ------------------
   */
 
 
 
  template <int dim, int spacedim>
  template <class Archive>
  void
  ParticleHandler<dim, spacedim>::serialize(Archive &ar, const unsigned int)
  {
    // Note that we do not serialize the particle data itself. 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.
    ar //&particles
      &global_number_of_particles &global_max_particles_per_cell
        &                          next_free_particle_index;
  }
 
 
 
  template <int dim, int spacedim>
  template <class VectorType>
  typename std::enable_if<
    std::is_convertible<VectorType *, Function<spacedim> *>::value ==
    false>::type
  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)
      {
        auto       new_point(displace_particles ? p.get_location() :
                                            Point<spacedim>());
        const auto id = p.get_id();
        for (unsigned int i = 0; i < spacedim; ++i)
          new_point[i] += input_vector[id * spacedim + i];
        p.set_location(new_point);
      }
    sort_particles_into_subdomains_and_cells();
  }
 
 
 
  template <int dim, int spacedim>
  template <class VectorType>
  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