particle_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2017 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_particles_particle_handler_h
#define dealii_particles_particle_handler_h

#include <deal.II/particles/particle.h>
#include <deal.II/particles/particle_iterator.h>
#include <deal.II/particles/property_pool.h>

#include <deal.II/distributed/tria.h>
#include <deal.II/base/subscriptor.h>
#include <deal.II/base/mpi.h>
#include <deal.II/base/array_view.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/fe/mapping.h>

#include <boost/serialization/map.hpp>
#include <boost/range/iterator_range.hpp>

DEAL_II_NAMESPACE_OPEN

#ifdef DEAL_II_WITH_P4EST

namespace Particles
{
  template <int dim, int spacedim=dim>
  class ParticleHandler: public Subscriptor
  {
  public:
    typedef ParticleIterator<dim,spacedim> particle_iterator;

    typedef boost::iterator_range<particle_iterator> particle_iterator_range;

    ParticleHandler();

    ParticleHandler(const parallel::distributed::Triangulation<dim,spacedim> &tria,
                    const Mapping<dim,spacedim> &mapping,
                    const unsigned int n_properties = 0);

    ~ParticleHandler();

    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(const bool serialization);

    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 data_offset;

#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

    void
    store_particles(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                    const typename Triangulation<dim,spacedim>::CellStatus status,
                    void *data) const;

    void
    load_particles(const typename Triangulation<dim,spacedim>::cell_iterator &cell,
                   const typename Triangulation<dim,spacedim>::CellStatus status,
                   const void *data);
  };

  /* -------------------------- 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;
  }
}

#endif // DEAL_II_WITH_P4EST

DEAL_II_NAMESPACE_CLOSE

#endif