agglomeration_handler.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) XXXX - YYYY by the polyDEAL authors
//
// This file is part of the polyDEAL library.
//
// Detailed license information governing the source code
// can be found in LICENSE.md at the top level directory.
//
// -----------------------------------------------------------------------------
#ifndef agglomeration_handler_h
#define agglomeration_handler_h
 
#include <deal.II/base/mpi.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/quadrature.h>
#include <deal.II/base/subscriptor.h>
 
#include <deal.II/distributed/shared_tria.h>
#include <deal.II/distributed/tria.h>
 
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/dofs/dof_tools.h>
 
#include <deal.II/fe/fe_dgp.h>
#include <deal.II/fe/fe_dgq.h>
#include <deal.II/fe/fe_nothing.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_fe_field.h>
#include <deal.II/fe/mapping_q.h>
 
#include <deal.II/grid/grid_tools_cache.h>
#include <deal.II/grid/tria.h>
 
#include <deal.II/hp/fe_collection.h>
 
#include <deal.II/lac/dynamic_sparsity_pattern.h>
#include <deal.II/lac/la_parallel_vector.h>
#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/sparsity_pattern.h>
#include <deal.II/lac/trilinos_sparse_matrix.h>
#include <deal.II/lac/vector.h>
 
#include <deal.II/meshworker/scratch_data.h>
 
#include <deal.II/non_matching/fe_immersed_values.h>
#include <deal.II/non_matching/immersed_surface_quadrature.h>
 
#include <agglomeration_iterator.h>
#include <agglomerator.h>
#include <mapping_box.h>
#include <utils.h>
 
#include <fstream>
 
using namespace dealii;
 
// Forward declarations
template <int dim, int spacedim>
class AgglomerationHandler;
 
namespace dealii
{
  namespace internal
  {
    template <int, int>
    class AgglomerationHandlerImplementation;
  } // namespace internal
} // namespace dealii
 
 

namespace dealii
{
  namespace internal
  {
    template <int dim, int spacedim>
    class PolytopeCache
    {
    public:
      PolytopeCache() = default;

      ~PolytopeCache() = default;
 
      void
      clear()
      {
        // clear all the members
        cell_face_at_boundary.clear();
        interface.clear();
        visited_cell_and_faces.clear();
      }

      mutable std::set<std::pair<types::global_cell_index, unsigned int>>
        visited_cell_and_faces;
 
 
      mutable std::set<std::pair<CellId, unsigned int>>
        visited_cell_and_faces_id;
 
 

      mutable std::map<
        std::pair<types::global_cell_index, unsigned int>,
        std::pair<bool,
                  typename Triangulation<dim, spacedim>::active_cell_iterator>>
        cell_face_at_boundary;

      mutable std::map<std::pair<CellId, unsigned int>, CellId>
        ghosted_master_id;

      mutable std::map<
        std::pair<CellId, CellId>,
        std::vector<
          std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator,
                    unsigned int>>>
        interface;
    };
  } // namespace internal
} // namespace dealii
 

template <int dim, int spacedim = dim>
class AgglomerationHandler : public Subscriptor
{
public:
  using agglomeration_iterator = AgglomerationIterator<dim, spacedim>;
 
  using AgglomerationContainer =
    typename AgglomerationIterator<dim, spacedim>::AgglomerationContainer;
 
 
  enum CellAgglomerationType
  {
    master = 0,
    slave  = 1
  };
 
 
 
  explicit AgglomerationHandler(
    const GridTools::Cache<dim, spacedim> &cached_tria);
 
  AgglomerationHandler() = default;
 
  ~AgglomerationHandler()
  {
    // disconnect the signal
    tria_listener.disconnect();
  }

  agglomeration_iterator
  begin() const;

  agglomeration_iterator
  begin();

  agglomeration_iterator
  end() const;

  agglomeration_iterator
  end();

  agglomeration_iterator
  last();

  IteratorRange<agglomeration_iterator>
  polytope_iterators() const;
 
  template <int, int>
  friend class AgglomerationIterator;
 
  template <int, int>
  friend class AgglomerationAccessor;

  void
  distribute_agglomerated_dofs(const FiniteElement<dim> &fe_space);

  void
  distribute_agglomerated_dofs(
    const hp::FECollection<dim, spacedim> &fe_collection_in);

  void
  initialize_fe_values(
    const Quadrature<dim>     &cell_quadrature = QGauss<dim>(1),
    const UpdateFlags         &flags           = UpdateFlags::update_default,
    const Quadrature<dim - 1> &face_quadrature = QGauss<dim - 1>(1),
    const UpdateFlags         &face_flags      = UpdateFlags::update_default);

  void
  initialize_fe_values(
    const hp::QCollection<dim> &cell_qcollection =
      hp::QCollection<dim>(QGauss<dim>(1)),
    const UpdateFlags              &flags = UpdateFlags::update_default,
    const hp::QCollection<dim - 1> &face_qcollection =
      hp::QCollection<dim - 1>(QGauss<dim - 1>(1)),
    const UpdateFlags &face_flags = UpdateFlags::update_default);

  template <typename SparsityPatternType, typename Number = double>
  void
  create_agglomeration_sparsity_pattern(
    SparsityPatternType             &sparsity_pattern,
    const AffineConstraints<Number> &constraints = AffineConstraints<Number>(),
    const bool                       keep_constrained_dofs = true,
    const types::subdomain_id subdomain_id = numbers::invalid_subdomain_id);

  agglomeration_iterator
  define_agglomerate(const AgglomerationContainer &cells);

  agglomeration_iterator
  define_agglomerate(const AgglomerationContainer &cells,
                     const unsigned int            fecollection_size);

  void
  define_agglomerate_with_check(const AgglomerationContainer &cells);
 
  inline const Triangulation<dim, spacedim> &
  get_triangulation() const;
 
  inline const FiniteElement<dim, spacedim> &
  get_fe() const;
 
  inline const Mapping<dim> &
  get_mapping() const;
 
  inline const MappingBox<dim> &
  get_agglomeration_mapping() const;
 
  inline const std::vector<BoundingBox<dim>> &
  get_local_bboxes() const;

  double
  get_mesh_size() const;
 
  inline types::global_cell_index
  cell_to_polytope_index(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
    const;
 
  inline decltype(auto)
  get_interface() const;

  template <typename CellIterator>
  inline bool
  is_master_cell(const CellIterator &cell) const;

  inline const std::vector<
    typename Triangulation<dim, spacedim>::active_cell_iterator> &
  get_slaves_of_idx(types::global_cell_index idx) const;
 
 
  inline const LinearAlgebra::distributed::Vector<float> &
  get_relationships() const;

  inline std::vector<
    typename Triangulation<dim, spacedim>::active_cell_iterator>
  get_agglomerate(
    const typename Triangulation<dim, spacedim>::active_cell_iterator
      &master_cell) const;

  template <class StreamType>
  void
  print_agglomeration(StreamType &out)
  {
    for (const auto &cell : tria->active_cell_iterators())
      out << "Cell with index: " << cell->active_cell_index()
          << " has associated value: "
          << master_slave_relationships[cell->global_active_cell_index()]
          << std::endl;
  }

  inline const DoFHandler<dim, spacedim> &
  get_dof_handler() const;

  unsigned int
  n_agglomerates() const;

  unsigned int
  n_agglomerated_faces_per_cell(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
    const;

  const FEValues<dim, spacedim> &
  reinit(const AgglomerationIterator<dim, spacedim> &polytope) const;

  const FEValuesBase<dim, spacedim> &
  reinit(const AgglomerationIterator<dim, spacedim> &polytope,
         const unsigned int                          face_index) const;

  std::pair<const FEValuesBase<dim, spacedim> &,
            const FEValuesBase<dim, spacedim> &>
  reinit_interface(const AgglomerationIterator<dim, spacedim> &polytope_in,
                   const AgglomerationIterator<dim, spacedim> &neigh_polytope,
                   const unsigned int                          local_in,
                   const unsigned int local_outside) const;

  Quadrature<dim>
  agglomerated_quadrature(
    const AgglomerationContainer &cells,
    const typename Triangulation<dim, spacedim>::active_cell_iterator
      &master_cell) const;
 

  inline bool
  at_boundary(
    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
    const unsigned int                                              f) const;
 
  inline unsigned int
  n_dofs_per_cell() const noexcept;
 
  inline types::global_dof_index
  n_dofs() const noexcept;
 
 

  inline const std::vector<typename Triangulation<dim>::active_face_iterator> &
  polytope_boundary(
    const typename Triangulation<dim>::active_cell_iterator &cell);
 

  DoFHandler<dim, spacedim> agglo_dh;

  DoFHandler<dim, spacedim> output_dh;
 
  std::unique_ptr<MappingBox<dim>> box_mapping;

  void
  setup_ghost_polytopes();
 
  void
  exchange_interface_values();
 
  // TODO: move it to private interface
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<Point<spacedim>>>>
    recv_qpoints;
 
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<double>>>
    recv_jxws;
 
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<Tensor<1, spacedim>>>>
    recv_normals;
 
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<std::vector<double>>>>
    recv_values;
 
  mutable std::map<types::subdomain_id,
                   std::map<std::pair<CellId, unsigned int>,
                            std::vector<std::vector<Tensor<1, spacedim>>>>>
    recv_gradients;

  inline const typename DoFHandler<dim, spacedim>::active_cell_iterator
  polytope_to_dh_iterator(const types::global_cell_index polytope_index) const;

  template <typename RtreeType>
  void
  connect_hierarchy(const CellsAgglomerator<dim, RtreeType> &agglomerator);

  inline const hp::FECollection<dim, spacedim> &
  get_fe_collection() const;

  inline bool
  used_fe_collection() const;
 
private:
  void
  initialize_agglomeration_data(
    const std::unique_ptr<GridTools::Cache<dim, spacedim>> &cache_tria);
 
  void
  update_agglomerate(
    AgglomerationContainer &polytope,
    const typename Triangulation<dim, spacedim>::active_cell_iterator
      &master_cell);

  void
  connect_to_tria_signals()
  {
    // First disconnect existing connections
    tria_listener.disconnect();
    tria_listener = tria->signals.any_change.connect(
      [&]() { this->initialize_agglomeration_data(this->cached_tria); });
  }

 
  inline typename Triangulation<dim, spacedim>::active_cell_iterator &
  is_slave_cell_of(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell);

  void
  create_bounding_box(const AgglomerationContainer &polytope);
 
 
  inline types::global_cell_index
  get_master_idx_of_cell(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
    const;

  inline bool
  are_cells_agglomerated(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
    const typename Triangulation<dim, spacedim>::active_cell_iterator
      &other_cell) const;

  void
  initialize_hp_structure();
 

  const FEValuesBase<dim, spacedim> &
  reinit_master(
    const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
    const unsigned int                                              face_number,
    std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
      &agglo_isv_ptr) const;
 

  template <typename CellIterator>
  inline bool
  is_slave_cell(const CellIterator &cell) const;
 

  void
  setup_connectivity_of_agglomeration();
 

  unsigned int n_agglomerations;
 

  LinearAlgebra::distributed::Vector<float> master_slave_relationships;

  std::map<types::global_cell_index,
           typename Triangulation<dim, spacedim>::active_cell_iterator>
    master_slave_relationships_iterators;
 
  using ScratchData = MeshWorker::ScratchData<dim, spacedim>;
 
  mutable std::vector<types::global_cell_index> number_of_agglomerated_faces;

  mutable std::map<
    const typename Triangulation<dim, spacedim>::active_cell_iterator,
    std::vector<typename Triangulation<dim>::active_face_iterator>>
    polygon_boundary;
 

  std::vector<BoundingBox<spacedim>> bboxes;

 
 
  // n_faces
  mutable std::map<types::subdomain_id, std::map<CellId, unsigned int>>
    local_n_faces;
 
  mutable std::map<types::subdomain_id, std::map<CellId, unsigned int>>
    recv_n_faces;
 
 
  // CellId (including slaves)
  mutable std::map<types::subdomain_id, std::map<CellId, CellId>>
    local_cell_ids_neigh_cell;
 
  mutable std::map<types::subdomain_id, std::map<CellId, CellId>>
    recv_cell_ids_neigh_cell;
 
 
  // send to neighborign rank the information that
  // - current polytope id
  // - face f
  // has the following neighboring id.
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::map<unsigned int, CellId>>>
    local_ghosted_master_id;
 
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::map<unsigned int, CellId>>>
    recv_ghosted_master_id;
 
  // CellIds from neighboring rank
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::map<unsigned int, bool>>>
    local_bdary_info;
 
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::map<unsigned int, bool>>>
    recv_bdary_info;
 
  // Exchange neighboring bounding boxes
  mutable std::map<types::subdomain_id, std::map<CellId, BoundingBox<dim>>>
    local_ghosted_bbox;
 
  mutable std::map<types::subdomain_id, std::map<CellId, BoundingBox<dim>>>
    recv_ghosted_bbox;
 
  // Exchange DoF indices with ghosted polytopes
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::vector<types::global_dof_index>>>
    local_ghost_dofs;
 
  mutable std::map<types::subdomain_id,
                   std::map<CellId, std::vector<types::global_dof_index>>>
    recv_ghost_dofs;
 
  // Exchange qpoints
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<Point<spacedim>>>>
    local_qpoints;
 
  // Exchange jxws
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<double>>>
    local_jxws;
 
  // Exchange normals
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<Tensor<1, spacedim>>>>
    local_normals;
 
  // Exchange values
  mutable std::map<
    types::subdomain_id,
    std::map<std::pair<CellId, unsigned int>, std::vector<std::vector<double>>>>
    local_values;
 
  mutable std::map<types::subdomain_id,
                   std::map<std::pair<CellId, unsigned int>,
                            std::vector<std::vector<Tensor<1, spacedim>>>>>
    local_gradients;
 
 

 
  ObserverPointer<const Triangulation<dim, spacedim>> tria;
 
  ObserverPointer<const Mapping<dim, spacedim>> mapping;
 
  std::unique_ptr<GridTools::Cache<dim, spacedim>> cached_tria;
 
  const MPI_Comm communicator;
 
  // The FiniteElement space we have on each cell. Currently supported types are
  // FE_DGQ and FE_DGP elements.
  std::unique_ptr<FiniteElement<dim>> fe;
 
  hp::FECollection<dim, spacedim> fe_collection;

  LinearAlgebra::distributed::Vector<double> euler_vector;
 

  mutable std::unique_ptr<ScratchData> standard_scratch;

  mutable std::unique_ptr<ScratchData> agglomerated_scratch;
 
 
  mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
    agglomerated_isv;
 
  mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
    agglomerated_isv_neigh;
 
  mutable std::unique_ptr<NonMatching::FEImmersedSurfaceValues<spacedim>>
    agglomerated_isv_bdary;
 
  boost::signals2::connection tria_listener;
 
  UpdateFlags agglomeration_flags;
 
  const UpdateFlags internal_agglomeration_flags =
    update_values | update_gradients | update_JxW_values |
    update_quadrature_points;
 
  UpdateFlags agglomeration_face_flags;
 
  const UpdateFlags internal_agglomeration_face_flags =
    update_quadrature_points | update_normal_vectors | update_values |
    update_gradients | update_JxW_values | update_inverse_jacobians;
 
  Quadrature<dim> agglomeration_quad;
 
  Quadrature<dim - 1> agglomeration_face_quad;
 
  // Associate the master cell to the slaves.
  std::unordered_map<
    types::global_cell_index,
    std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>>
    master2slaves;
 
  // Map the master cell index with the polytope index
  std::map<types::global_cell_index, types::global_cell_index> master2polygon;
 
 
  std::vector<typename Triangulation<dim>::active_cell_iterator>
    master_disconnected;
 
  // Dummy FiniteElement objects needed only to generate quadratures

  FE_Nothing<dim, spacedim> dummy_fe;

  std::unique_ptr<FEValues<dim, spacedim>> no_values;

  std::unique_ptr<FEFaceValues<dim, spacedim>> no_face_values;

  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
    master_cells_container;
 
  friend class internal::AgglomerationHandlerImplementation<dim, spacedim>;
 
  internal::PolytopeCache<dim, spacedim> polytope_cache;

  bool hybrid_mesh;
 
  std::map<std::pair<types::global_cell_index, types::global_cell_index>,
           std::vector<types::global_cell_index>>
    parent_child_info;
 
  unsigned int present_extraction_level;
 
  // Support for hp::FECollection
  bool is_hp_collection = false; // Indicates whether hp::FECollection is used
  std::unique_ptr<hp::FECollection<dim, spacedim>>
    hp_fe_collection; // External input FECollection
 
  // Stores quadrature rules; these QCollections should have the same size as
  // hp_fe_collection
  hp::QCollection<dim>     agglomeration_quad_collection;
  hp::QCollection<dim - 1> agglomeration_face_quad_collection;
 
  hp::MappingCollection<dim>
    mapping_collection; // Contains only one mapping object
  hp::FECollection<dim, spacedim>
    dummy_fe_collection; // Similar to dummy_fe, but as an FECollection
                         // containing only dummy_fe
  // Note: The above two variables provide an hp::FECollection interface but
  // actually contain only one element each.
 
  // Analogous to no_values and no_face_values, but used when different cells
  // employ different FEs or quadratures
  std::unique_ptr<hp::FEValues<dim, spacedim>>     hp_no_values;
  std::unique_ptr<hp::FEFaceValues<dim, spacedim>> hp_no_face_values;
};
 
 
 
// ------------------------------ inline functions -------------------------
template <int dim, int spacedim>
inline const FiniteElement<dim, spacedim> &
AgglomerationHandler<dim, spacedim>::get_fe() const
{
  return *fe;
}
 
 
 
template <int dim, int spacedim>
inline const Mapping<dim> &
AgglomerationHandler<dim, spacedim>::get_mapping() const
{
  return *mapping;
}
 
 
 
template <int dim, int spacedim>
inline const MappingBox<dim> &
AgglomerationHandler<dim, spacedim>::get_agglomeration_mapping() const
{
  return *box_mapping;
}
 
 
 
template <int dim, int spacedim>
inline const Triangulation<dim, spacedim> &
AgglomerationHandler<dim, spacedim>::get_triangulation() const
{
  return *tria;
}
 
 
template <int dim, int spacedim>
inline const std::vector<BoundingBox<dim>> &
AgglomerationHandler<dim, spacedim>::get_local_bboxes() const
{
  return bboxes;
}
 
 
 
template <int dim, int spacedim>
inline types::global_cell_index
AgglomerationHandler<dim, spacedim>::cell_to_polytope_index(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell) const
{
  return master2polygon.at(cell->active_cell_index());
}
 
 
 
template <int dim, int spacedim>
inline decltype(auto)
AgglomerationHandler<dim, spacedim>::get_interface() const
{
  return polytope_cache.interface;
}
 
 
 
template <int dim, int spacedim>
inline const LinearAlgebra::distributed::Vector<float> &
AgglomerationHandler<dim, spacedim>::get_relationships() const
{
  return master_slave_relationships;
}
 
 
 
template <int dim, int spacedim>
inline std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>
AgglomerationHandler<dim, spacedim>::get_agglomerate(
  const typename Triangulation<dim, spacedim>::active_cell_iterator
    &master_cell) const
{
  Assert(is_master_cell(master_cell), ExcInternalError());
  auto agglomeration = get_slaves_of_idx(master_cell->active_cell_index());
  agglomeration.push_back(master_cell);
  return agglomeration;
}
 
 
 
template <int dim, int spacedim>
inline const DoFHandler<dim, spacedim> &
AgglomerationHandler<dim, spacedim>::get_dof_handler() const
{
  return agglo_dh;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<
  typename Triangulation<dim, spacedim>::active_cell_iterator> &
AgglomerationHandler<dim, spacedim>::get_slaves_of_idx(
  types::global_cell_index idx) const
{
  return master2slaves.at(idx);
}
 
 
 
template <int dim, int spacedim>
template <typename CellIterator>
inline bool
AgglomerationHandler<dim, spacedim>::is_master_cell(
  const CellIterator &cell) const
{
  return master_slave_relationships[cell->global_active_cell_index()] == -1;
}
 
 

template <int dim, int spacedim>
template <typename CellIterator>
inline bool
AgglomerationHandler<dim, spacedim>::is_slave_cell(
  const CellIterator &cell) const
{
  return master_slave_relationships[cell->global_active_cell_index()] >= 0;
}
 
 
 
template <int dim, int spacedim>
inline bool
AgglomerationHandler<dim, spacedim>::at_boundary(
  const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
  const unsigned int face_index) const
{
  Assert(!is_slave_cell(cell),
         ExcMessage("This function should not be called for a slave cell."));
 
  return polytope_cache.cell_face_at_boundary
    .at({master2polygon.at(cell->active_cell_index()), face_index})
    .first;
}
 
 
template <int dim, int spacedim>
inline unsigned int
AgglomerationHandler<dim, spacedim>::n_dofs_per_cell() const noexcept
{
  return fe->n_dofs_per_cell();
}
 
 
 
template <int dim, int spacedim>
inline types::global_dof_index
AgglomerationHandler<dim, spacedim>::n_dofs() const noexcept
{
  return agglo_dh.n_dofs();
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<typename Triangulation<dim>::active_face_iterator> &
AgglomerationHandler<dim, spacedim>::polytope_boundary(
  const typename Triangulation<dim>::active_cell_iterator &cell)
{
  return polygon_boundary[cell];
}
 
 
 
template <int dim, int spacedim>
inline typename Triangulation<dim, spacedim>::active_cell_iterator &
AgglomerationHandler<dim, spacedim>::is_slave_cell_of(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
{
  return master_slave_relationships_iterators.at(cell->active_cell_index());
}
 
 
 
template <int dim, int spacedim>
inline types::global_cell_index
AgglomerationHandler<dim, spacedim>::get_master_idx_of_cell(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell) const
{
  auto idx = master_slave_relationships[cell->global_active_cell_index()];
  if (idx == -1)
    return cell->global_active_cell_index();
  else
    return static_cast<types::global_cell_index>(idx);
}
 
 
 
template <int dim, int spacedim>
inline bool
AgglomerationHandler<dim, spacedim>::are_cells_agglomerated(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
  const typename Triangulation<dim, spacedim>::active_cell_iterator &other_cell)
  const
{
  // if different subdomain, then **by construction** they will not be together
  // if (cell->subdomain_id() != other_cell->subdomain_id())
  //   return false;
  // else
  return (get_master_idx_of_cell(cell) == get_master_idx_of_cell(other_cell));
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
AgglomerationHandler<dim, spacedim>::n_agglomerates() const
{
  return n_agglomerations;
}
 
 
 
template <int dim, int spacedim>
inline const typename DoFHandler<dim, spacedim>::active_cell_iterator
AgglomerationHandler<dim, spacedim>::polytope_to_dh_iterator(
  const types::global_cell_index polytope_index) const
{
  return master_cells_container[polytope_index]->as_dof_handler_iterator(
    agglo_dh);
}
 
 
 
template <int dim, int spacedim>
AgglomerationIterator<dim, spacedim>
AgglomerationHandler<dim, spacedim>::begin() const
{
  Assert(n_agglomerations > 0,
         ExcMessage("No agglomeration has been performed."));
  return {*master_cells_container.begin(), this};
}
 
 
 
template <int dim, int spacedim>
AgglomerationIterator<dim, spacedim>
AgglomerationHandler<dim, spacedim>::begin()
{
  Assert(n_agglomerations > 0,
         ExcMessage("No agglomeration has been performed."));
  return {*master_cells_container.begin(), this};
}
 
 
 
template <int dim, int spacedim>
AgglomerationIterator<dim, spacedim>
AgglomerationHandler<dim, spacedim>::end() const
{
  Assert(n_agglomerations > 0,
         ExcMessage("No agglomeration has been performed."));
  return {*master_cells_container.end(), this};
}
 
 
 
template <int dim, int spacedim>
AgglomerationIterator<dim, spacedim>
AgglomerationHandler<dim, spacedim>::end()
{
  Assert(n_agglomerations > 0,
         ExcMessage("No agglomeration has been performed."));
  return {*master_cells_container.end(), this};
}
 
 
 
template <int dim, int spacedim>
AgglomerationIterator<dim, spacedim>
AgglomerationHandler<dim, spacedim>::last()
{
  Assert(n_agglomerations > 0,
         ExcMessage("No agglomeration has been performed."));
  return {master_cells_container.back(), this};
}
 
 
 
template <int dim, int spacedim>
IteratorRange<
  typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator>
AgglomerationHandler<dim, spacedim>::polytope_iterators() const
{
  return IteratorRange<
    typename AgglomerationHandler<dim, spacedim>::agglomeration_iterator>(
    begin(), end());
}
 
template <int dim, int spacedim>
template <typename RtreeType>
void
AgglomerationHandler<dim, spacedim>::connect_hierarchy(
  const CellsAgglomerator<dim, RtreeType> &agglomerator)
{
  parent_child_info        = agglomerator.parent_node_to_children_nodes;
  present_extraction_level = agglomerator.extraction_level;
}
 
template <int dim, int spacedim>
inline const hp::FECollection<dim, spacedim> &
AgglomerationHandler<dim, spacedim>::get_fe_collection() const
{
  return *hp_fe_collection;
}
 
template <int dim, int spacedim>
inline bool
AgglomerationHandler<dim, spacedim>::used_fe_collection() const
{
  return is_hp_collection;
}
 
 
#endif