obstacle_data_structure.hpp

Go to the documentation of this file.

#pragma once
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/point.h>
 
#include <deal.II/grid/filtered_iterator.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/tria_iterator_selector.h>
 
#include <deal.II/matrix_free/matrix_free.h>
 
#include <deal.II/numerics/solution_transfer.h>
 
#include <deal.II/particles/particle_handler.h>
#include <deal.II/particles/property_pool.h>
 
#include <meltpooldg/particles/particle_iterator.hpp>
#include <meltpooldg/post_processing/postprocessor.hpp>
#include <meltpooldg/utilities/triangulation_utils.hpp>
 
#include <boost/container/small_vector.hpp>
 
#include <functional>
#include <ranges>
#include <vector>
 
namespace MeltPoolDG
{
  template <int dim, typename number, typename ObstacleType>
  struct CellListParticleHandler;
 
 
  template <int dim, typename number, typename ObstacleType>
  struct NeighborListUpdateTracker
  {
  public:
    explicit NeighborListUpdateTracker(
      const CellListParticleHandler<dim, number, ObstacleType> &obstacle_data_structure);
 
    bool
    update_required() const;
 
    void
    reinit_after_update(const number new_skin_thickness                 = 0,
                        const number new_max_displacement_before_update = 0);
 
  private:
    std::vector<dealii::Point<dim, number>> previous_obstacle_locations;
 
    const CellListParticleHandler<dim, number, ObstacleType> &obstacle_data_structure;
 
    number skin_thickness;
 
    number max_displacement_before_update = 0;
  };
 
 
  template <int dim, typename number, typename ObstacleType>
  struct CellListParticleHandler
  {
  public:
    enum class NotifyEvent
    {
      ObserverInitialization,
 
      UpdateGhostParticleProperties,
 
      SortParticlesIntoSubdomainsAndCells
    };
 
    explicit CellListParticleHandler(const dealii::Triangulation<dim> &triangulation,
                                     const dealii::Mapping<dim>       &mapping);
 
    ~CellListParticleHandler();
 
    void
    initialize();
 
    template <typename ObstacleContainer>
    void
    get_obstacles_in_cell(const dealii::TriaIterator<dealii::CellAccessor<dim>> &cell,
                          ObstacleContainer                                     &particles);
 
    void
    broadcast_global_particles() const;
 
    void
    prepare_for_coarsening_and_refinement();
 
    void
    unpack_after_coarsening_and_refinement();
 
    void
    insert_global_particles(const std::vector<dealii::Point<dim, number>> &obstacle_locations,
                            const std::vector<std::vector<number>>        &obstacle_properties);
 
    void
    register_particle_output(Postprocessor<dim, number> &postprocessor) const;
 
    std::ranges::subrange<ParticleIterator<dim, number>>
    locally_owned_particle_range() const;
 
    std::ranges::subrange<MeltPoolDG::ParticleIterator<dim, number>>
    ghost_particle_range() const;
 
    typename std::ranges::subrange<ParticleIterator<dim, number>>
    particles_in_cell(typename dealii::Triangulation<dim>::active_cell_iterator cell) const;
 
    boost::container::small_vector<DEMParticleAccessor<dim, number>, 3 * dim>
    find_particles_in_neighborhood(const DEMParticleAccessor<dim, number> &particle,
                                   const number                            relative_tolerance);
 
    void
    compress();
 
    void
    update_ghost_particle_properties();
 
    void
    sort_particles_into_subdomains_and_cells();
 
    void
    auto_update_particle_cache();
 
    void
    prepare_for_serialization();
 
    template <class Archive>
    void
    serialize(Archive &ar, const unsigned int version);
 
    void
    deserialize();
 
    unsigned int
    n_global_particles() const;
 
    unsigned int
    n_locally_owned_particles() const;
 
    unsigned int
    n_ghost_particles() const;
 
    MPI_Comm
    get_mpi_communicator() const;
 
    void
    subscribe(std::function<void(CellListParticleHandler &, const NotifyEvent)> callback);
 
  private:
    dealii::Particles::ParticleHandler<dim> obstacle_handler;
 
    mutable dealii::Particles::PropertyPool<dim> ghost_particles_property_pool;
 
    MPI_Comm mpi_communicator = MPI_COMM_WORLD;
 
    std::unordered_map<dealii::types::particle_index, dealii::Particles::ParticleIterator<dim>>
      particle_id_to_iterator_cache;
 
    boost::signals2::signal<void(CellListParticleHandler &, const NotifyEvent &)> notify_signal;
 
    struct
    {
      std::vector<std::vector<dealii::Particles::ParticleIterator<dim>>> locally_owned_particles;
 
      std::vector<std::vector<typename dealii::Particles::PropertyPool<dim>::Handle>>
        ghost_particles;
 
      int cell_level = 0;
 
      number global_max_particle_radius = 0;
    } cell_particle_cache;
 
    struct TriangulationLevelCache
    {
      TriangulationLevelCache(const dealii::Triangulation<dim> &triangulation)
        : communication_pattern(triangulation)
      {}
 
      LevelCommunicationPattern<dim> communication_pattern;
 
      LevelAdjacentCellsCache<dim> adjacent_cells;
 
      int level = 0;
    } triangulation_level_cache;
 
    struct GhostParticleUpdateCache
    {
      void
      reset(MPI_Comm mpi_communicator)
      {
        particles_to_send.clear();
        particles_to_send.resize(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator));
        n_particles_to_receive.clear();
        n_particles_to_receive.resize(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator), 0);
        rank_ghost_particle_start_handle.clear();
        rank_ghost_particle_start_handle.resize(
          dealii::Utilities::MPI::n_mpi_processes(mpi_communicator), 0);
      }
 
      // index is the rank to whom to send the corresponding particles. Must be sorted with
      // respect to global particle id.
      std::vector<std::vector<dealii::Particles::ParticleIterator<dim>>> particles_to_send;
 
      // index is the rank from which to receive the corresponding particles
      std::vector<unsigned int> n_particles_to_receive;
 
      // index is the rank from which to receive the corresponding particles. The value is the
      // first handle of the ghost particles received from that rank in the global property pool.
      std::vector<unsigned int> rank_ghost_particle_start_handle;
 
      std::vector<std::vector<typename dealii::Particles::PropertyPool<dim>::Handle>>
        rank_to_ghost_handle;
    } ghost_particle_update_cache;
 
    NeighborListUpdateTracker<dim, number, ObstacleType> neighbor_list_update_tracker;
 
    void
    deregister_property_pool() const;
 
    dealii::TriaIterator<dealii::CellAccessor<dim>>
    find_particle_cache_cell(const dealii::TriaIterator<dealii::CellAccessor<dim>> &cell) const;
 
    template <typename ObstacleContainer>
    void
    find_relevant_particles(const dealii::TriaIterator<dealii::CellAccessor<dim>> &cell,
                            ObstacleContainer                                     &particles);
 
    number
    compute_max_particle_radius() const;
 
    void
    notify(const NotifyEvent &event);
  };
 
  template <int dim, typename number, typename ObstacleType>
  template <class Archive>
  void
  CellListParticleHandler<dim, number, ObstacleType>::serialize(Archive           &ar,
                                                                const unsigned int version)
  {
    obstacle_handler.serialize(ar, version);
  }
 
  template <int dim, typename number, typename ObstacleType>
  template <typename ObstacleContainer>
  void
  CellListParticleHandler<dim, number, ObstacleType>::get_obstacles_in_cell(
    const dealii::TriaIterator<dealii::CellAccessor<dim>> &cell,
    ObstacleContainer                                     &particles)
  {
    std::size_t current_size = particles.size();
 
    // Find all relevant particles for the given cell and add them to the provided container.
    find_relevant_particles(find_particle_cache_cell(cell), particles);
 
    // Remove duplicates from the container while ensuring that the particles which were already
    // present before the function call are not removed.
    auto old_end   = particles.begin() + current_size;
    auto write_ptr = old_end;
    for (auto read_ptr = old_end; read_ptr != particles.end(); ++read_ptr)
      {
        auto found_it = std::find_if(particles.begin(), old_end, [&read_ptr](const auto &particle) {
          return particle.id() == read_ptr->id();
        });
        if (found_it == old_end)
          {
            if (write_ptr->id() != read_ptr->id())
              *write_ptr = std::move(*read_ptr);
            ++write_ptr;
          }
      }
    particles.erase(write_ptr, particles.end());
  }
 
  template <int dim, typename number, typename ObstacleType>
  template <typename ObstacleContainer>
  void
  CellListParticleHandler<dim, number, ObstacleType>::find_relevant_particles(
    const dealii::TriaIterator<dealii::CellAccessor<dim>> &cell,
    ObstacleContainer                                     &particles)
  {
    Assert(cell->level() == cell_particle_cache.cell_level,
           dealii::ExcMessage(
             "The function only supports searching for particles in cells on the same level as the "
             "level on which particles are cached."));
 
    const auto push_back_particles_from_cell = [&](const auto &current_cell) {
      for (dealii::Particles::ParticleIterator<dim> &particle :
           cell_particle_cache.locally_owned_particles[current_cell->index()])
        {
          Assert(particle.state() == dealii::IteratorState::valid, dealii::ExcInternalError());
          particles.emplace_back(*particle);
        }
 
      for (const typename dealii::Particles::PropertyPool<dim>::Handle particle_handle :
           cell_particle_cache.ghost_particles[current_cell->index()])
        {
          Assert(particle_handle < ghost_particles_property_pool.n_registered_slots(),
                 dealii::ExcInternalError());
          particles.emplace_back(ghost_particles_property_pool, particle_handle);
        }
    };
 
    // Adjacent cells
    for (const dealii::TriaIterator<dealii::CellAccessor<dim>> &current_cell :
         triangulation_level_cache.adjacent_cells.get_adjacent_cells(cell))
      {
        push_back_particles_from_cell(current_cell);
      }
 
    // Current cell
    push_back_particles_from_cell(cell);
  }
} // namespace MeltPoolDG