particle_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 - 2019 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/array_view.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/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

#ifdef DEAL_II_WITH_P4EST

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 parallel::distributed::Triangulation<dim, spacedim> &tria,
      const Mapping<dim, spacedim> &                             mapping,
      const unsigned int n_properties = 0);

    virtual ~ParticleHandler() override = default;

    void
    initialize(const parallel::distributed::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);

    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;

    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 parallel::distributed::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;

#  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;
  }
} // namespace Particles

#endif // DEAL_II_WITH_P4EST

DEAL_II_NAMESPACE_CLOSE

#endif