agglomeration_accessor.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_accessor_h
#define agglomeration_accessor_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/iterator_range.h>
 
#include <deal.II/grid/filtered_iterator.h>
 
#include <vector>
 
using namespace dealii;
 
 
// Forward declarations
#ifndef DOXYGEN
template <int, int>
class AgglomerationHandler;
template <int, int>
class AgglomerationIterator;
#endif
 

template <int dim, int spacedim = dim>
class AgglomerationAccessor
{
public:
  using AgglomerationContainer =
    std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator>;
 

  void
  get_dof_indices(std::vector<types::global_dof_index> &) const;

  unsigned int
  n_faces() const;

  unsigned int
  n_agglomerated_faces() const;

  const AgglomerationIterator<dim, spacedim>
  neighbor(const unsigned int f) const;

  unsigned int
  neighbor_of_agglomerated_neighbor(const unsigned int f) const;

  bool
  at_boundary(const unsigned int f) const;

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

  double
  volume() const;

  double
  diameter() const;

  AgglomerationContainer
  get_agglomerate() const;

  const BoundingBox<dim> &
  get_bounding_box() const;

  types::global_cell_index
  index() const;

  typename DoFHandler<dim, spacedim>::active_cell_iterator
  as_dof_handler_iterator(const DoFHandler<dim, spacedim> &dof_handler) const;

  unsigned int
  n_background_cells() const;
 
  /* Returns true if this polygon is owned by the current processor. On a serial
   * Triangulation this returs always true, but may yield false for a
   * parallel::distributed::Triangulation.
   */
  bool
  is_locally_owned() const;

  CellId
  id() const;

  types::subdomain_id
  subdomain_id() const;

  inline const std::vector<types::global_cell_index> &
  children() const;

  const FiniteElement<dim, spacedim> &
  get_fe() const;

  void
  set_active_fe_index(const types::fe_index index) const;

  types::fe_index
  active_fe_index() const;
 
private:
  AgglomerationAccessor();

  AgglomerationAccessor(
    const typename Triangulation<dim, spacedim>::active_cell_iterator
                                              &master_cell,
    const AgglomerationHandler<dim, spacedim> *ah);

  AgglomerationAccessor(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
    const CellId                                                      &cell_id,
    const AgglomerationHandler<dim, spacedim>                         *ah);

  ~AgglomerationAccessor() = default;
 

  typename Triangulation<dim, spacedim>::active_cell_iterator master_cell;

  types::global_cell_index present_index;

  CellId present_id;

  types::subdomain_id present_subdomain_id;

  AgglomerationHandler<dim, spacedim> *handler;

  bool
  operator==(const AgglomerationAccessor<dim, spacedim> &other) const;

  bool
  operator!=(const AgglomerationAccessor<dim, spacedim> &other) const;

  void
  next();

  void
  prev();

  const AgglomerationContainer &
  get_slaves() const;
 
  unsigned int
  n_agglomerated_faces_per_cell(
    const typename Triangulation<dim, spacedim>::active_cell_iterator &cell)
    const;
 
  template <int, int>
  friend class AgglomerationIterator;
};
 
 
 
template <int dim, int spacedim>
unsigned int
AgglomerationAccessor<dim, spacedim>::n_agglomerated_faces_per_cell(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell) const
{
  unsigned int n_neighbors = 0;
  for (const auto &f : cell->face_indices())
    {
      const auto &neighboring_cell = cell->neighbor(f);
      if ((cell->face(f)->at_boundary()) ||
          (neighboring_cell->is_active() &&
           !handler->are_cells_agglomerated(cell, neighboring_cell)))
        {
          ++n_neighbors;
        }
    }
  return n_neighbors;
}
 
 
 
template <int dim, int spacedim>
unsigned int
AgglomerationAccessor<dim, spacedim>::n_faces() const
{
  Assert(!handler->is_slave_cell(master_cell),
         ExcMessage("You cannot pass a slave cell."));
  return handler->number_of_agglomerated_faces[present_index];
}
 
 
 
template <int dim, int spacedim>
const AgglomerationIterator<dim, spacedim>
AgglomerationAccessor<dim, spacedim>::neighbor(const unsigned int f) const
{
  if (!at_boundary(f))
    {
      if (master_cell->is_ghost())
        {
          // The following path is needed when the present function is called
          // from neighbor_of_neighbor()
 
          const unsigned int sender_rank = master_cell->subdomain_id();
 
          const CellId &master_id_ghosted_neighbor =
            handler->recv_ghosted_master_id.at(sender_rank)
              .at(present_id)
              .at(f);
 
          // Use the id of the master cell to uniquely identify the neighboring
          // agglomerate
 
          return {master_cell,
                  master_id_ghosted_neighbor,
                  handler}; // dummy master?
        }
 
      const types::global_cell_index polytope_index =
        handler->master2polygon.at(master_cell->active_cell_index());
 
      // #ifdef AGGLO_DEBUG
      //       std::cout << "Poly index: " << polytope_index << std::endl;
      //       std::cout << "Face index " << f << std::endl;
 
      //       if (Utilities::MPI::this_mpi_process(handler->communicator) == 1)
      //         for (const auto &[key, value] :
      //              handler->polytope_cache.cell_face_at_boundary)
      //           {
      //             std::cout << "p.f = " << key.first << std::endl;
      //             std::cout << "p.s = " << key.second << std::endl;
      //           }
      // #endif
 
      const auto &neigh =
        handler->polytope_cache.cell_face_at_boundary.at({polytope_index, f})
          .second;
 
 
      if (neigh->is_locally_owned())
        {
          typename DoFHandler<dim, spacedim>::active_cell_iterator cell_dh(
            *neigh, &(handler->agglo_dh));
          return {cell_dh, handler};
        }
      else
        {
          // #ifdef AGGLO_DEBUG
          //           std::cout << "before calling ghosted_masteer_id" <<
          //           std::endl; std::cout << "present_id = " << present_id <<
          //           std::endl; std::cout << "che master cell? " <<
          //           master_cell << std::endl;
 
          //           if
          //           (Utilities::MPI::this_mpi_process(handler->communicator)
          //           == 1)
          //             for (const auto &[key, value] :
          //                  handler->polytope_cache.ghosted_master_id)
          //               {
          //                 std::cout << "poly_id = " << key.first <<
          //                 std::endl; std::cout << "face_idx = " << key.second
          //                 << std::endl;
          //               }
 
          // #endif
          // Get master_id from the neighboring ghost polytope. This uniquely
          // identifies the neighboring polytope among all processors.
          const CellId &master_id_neighbor =
            handler->polytope_cache.ghosted_master_id.at({present_id, f});
 
          // Use the id of the master cell to uniquely identify the neighboring
          // agglomerate
          return {neigh, master_id_neighbor, handler};
        }
    }
  else
    {
      return {};
    }
}
 
 
 
template <int dim, int spacedim>
unsigned int
AgglomerationAccessor<dim, spacedim>::neighbor_of_agglomerated_neighbor(
  const unsigned int f) const
{
  // First, make sure it's not a boundary face.
  if (!at_boundary(f))
    {
      const auto &neigh_polytope =
        neighbor(f); // returns the neighboring master and id
 
      AssertThrow(neigh_polytope.state() == IteratorState::valid,
                  ExcInternalError());
 
      unsigned int n_faces_agglomerated_neighbor;
 
      // if it is locally owned, retrieve the number of faces
      if (neigh_polytope->is_locally_owned())
        {
          n_faces_agglomerated_neighbor = neigh_polytope->n_faces();
        }
      else
        {
          // The neighboring polytope is not locally owned. We need to get the
          // number of its faces from the neighboring rank.
 
          // First, retrieve the CellId of the neighboring polytope.
          const CellId &master_id_neighbor = neigh_polytope->id();
 
          // Then, get the neighboring rank
          const unsigned int sender_rank = neigh_polytope->subdomain_id();
 
          // From the neighboring rank, use the CellId of the neighboring
          // polytope to get the number of its faces.
          n_faces_agglomerated_neighbor =
            handler->recv_n_faces.at(sender_rank).at(master_id_neighbor);
        }
 
 
      // Loop over all faces of neighboring agglomerate
      for (unsigned int f_out = 0; f_out < n_faces_agglomerated_neighbor;
           ++f_out)
        {
          // Check if same CellId
          if (neigh_polytope->neighbor(f_out).state() == IteratorState::valid)
            if (neigh_polytope->neighbor(f_out)->id() == present_id)
              return f_out;
        }
      return numbers::invalid_unsigned_int;
    }
  else
    {
      // Face is at boundary
      return numbers::invalid_unsigned_int;
    }
}
 
// ------------------------------ inline functions -------------------------
 
template <int dim, int spacedim>
inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor()
{}
 
 
 
template <int dim, int spacedim>
inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &cell,
  const AgglomerationHandler<dim, spacedim>                         *ah)
{
  handler = const_cast<AgglomerationHandler<dim, spacedim> *>(ah);
  if (&(*handler->master_cells_container.end()) == std::addressof(cell))
    {
      present_index        = handler->master_cells_container.size();
      master_cell          = *handler->master_cells_container.end();
      present_id           = CellId(); // invalid id (TODO)
      present_subdomain_id = numbers::invalid_subdomain_id;
    }
  else
    {
      present_index = handler->master2polygon.at(cell->active_cell_index());
      master_cell   = cell;
      present_id    = master_cell->id();
      present_subdomain_id = master_cell->subdomain_id();
    }
}
 
 
 
template <int dim, int spacedim>
inline AgglomerationAccessor<dim, spacedim>::AgglomerationAccessor(
  const typename Triangulation<dim, spacedim>::active_cell_iterator &neigh_cell,
  const CellId                              &master_cell_id,
  const AgglomerationHandler<dim, spacedim> *ah)
{
  Assert(neigh_cell->is_ghost(), ExcInternalError());
  // neigh_cell is ghosted
 
  handler       = const_cast<AgglomerationHandler<dim, spacedim> *>(ah);
  master_cell   = neigh_cell;
  present_index = numbers::invalid_unsigned_int;
  // neigh_cell is ghosted, use the CellId of that agglomerate
  present_id           = master_cell_id;
  present_subdomain_id = master_cell->subdomain_id();
}
 
 
 
template <int dim, int spacedim>
inline void
AgglomerationAccessor<dim, spacedim>::get_dof_indices(
  std::vector<types::global_dof_index> &dof_indices) const
{
  Assert(dof_indices.size() > 0,
         ExcMessage(
           "The vector of DoFs indices must be already properly resized."));
  if (is_locally_owned())
    {
      // Forward the call to the master cell
      typename DoFHandler<dim, spacedim>::cell_iterator master_cell_dh(
        *master_cell, &(handler->agglo_dh));
      master_cell_dh->get_dof_indices(dof_indices);
    }
  else
    {
      const std::vector<types::global_dof_index> &recv_dof_indices =
        handler->recv_ghost_dofs.at(present_subdomain_id).at(present_id);
 
      std::copy(recv_dof_indices.cbegin(),
                recv_dof_indices.cend(),
                dof_indices.begin());
    }
}
 
 
 
template <int dim, int spacedim>
inline typename AgglomerationAccessor<dim, spacedim>::AgglomerationContainer
AgglomerationAccessor<dim, spacedim>::get_agglomerate() const
{
  auto agglomeration = get_slaves();
  agglomeration.push_back(master_cell);
  return agglomeration;
}
 
 
 
template <int dim, int spacedim>
inline const std::vector<typename Triangulation<dim>::active_face_iterator> &
AgglomerationAccessor<dim, spacedim>::polytope_boundary() const
{
  return handler->polygon_boundary[master_cell];
}
 
 
 
template <int dim, int spacedim>
inline double
AgglomerationAccessor<dim, spacedim>::diameter() const
{
  Assert(!handler->is_slave_cell(master_cell),
         ExcMessage("The present function cannot be called for slave cells."));
 
  if (handler->is_master_cell(master_cell))
    {
      // Get the bounding box associated with the master cell
      const auto &bdary_pts =
        handler->bboxes[present_index].get_boundary_points();
      return (bdary_pts.second - bdary_pts.first).norm();
    }
  else
    {
      // Standard deal.II way to get the measure of a cell.
      return master_cell->diameter();
    }
}
 
 
 
template <int dim, int spacedim>
inline const BoundingBox<dim> &
AgglomerationAccessor<dim, spacedim>::get_bounding_box() const
{
  if (is_locally_owned())
    return handler->bboxes[present_index];
  else
    return handler->recv_ghosted_bbox.at(present_subdomain_id).at(present_id);
}
 
 
 
template <int dim, int spacedim>
inline double
AgglomerationAccessor<dim, spacedim>::volume() const
{
  Assert(!handler->is_slave_cell(master_cell),
         ExcMessage("The present function cannot be called for slave cells."));
 
  if (handler->is_master_cell(master_cell))
    {
      return handler->bboxes[present_index].volume();
    }
  else
    {
      return master_cell->measure();
    }
}
 
 
 
template <int dim, int spacedim>
inline void
AgglomerationAccessor<dim, spacedim>::next()
{
  // Increment the present index and update the polytope
  ++present_index;
 
  // Make sure not to query the CellId if it's past the last
  if (present_index < handler->master_cells_container.size())
    {
      master_cell          = handler->master_cells_container[present_index];
      present_id           = master_cell->id();
      present_subdomain_id = master_cell->subdomain_id();
    }
}
 
 
 
template <int dim, int spacedim>
inline void
AgglomerationAccessor<dim, spacedim>::prev()
{
  // Decrement the present index and update the polytope
  --present_index;
  master_cell = handler->master_cells_container[present_index];
  present_id  = master_cell->id();
}
 
 
template <int dim, int spacedim>
inline bool
AgglomerationAccessor<dim, spacedim>::operator==(
  const AgglomerationAccessor<dim, spacedim> &other) const
{
  return present_index == other.present_index;
}
 
template <int dim, int spacedim>
inline bool
AgglomerationAccessor<dim, spacedim>::operator!=(
  const AgglomerationAccessor<dim, spacedim> &other) const
{
  return !(*this == other);
}
 
 
 
template <int dim, int spacedim>
inline types::global_cell_index
AgglomerationAccessor<dim, spacedim>::index() const
{
  return present_index;
}
 
 
 
template <int dim, int spacedim>
typename DoFHandler<dim, spacedim>::active_cell_iterator
AgglomerationAccessor<dim, spacedim>::as_dof_handler_iterator(
  const DoFHandler<dim, spacedim> &dof_handler) const
{
  // Forward the call to the master cell using the right DoFHandler.
  return master_cell->as_dof_handler_iterator(dof_handler);
}
 
 
 
template <int dim, int spacedim>
inline const typename AgglomerationAccessor<dim,
                                            spacedim>::AgglomerationContainer &
AgglomerationAccessor<dim, spacedim>::get_slaves() const
{
  return handler->master2slaves.at(master_cell->active_cell_index());
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
AgglomerationAccessor<dim, spacedim>::n_background_cells() const
{
  AssertThrow(get_agglomerate().size() > 0, ExcMessage("Empty agglomeration."));
  return get_agglomerate().size();
}
 
 
 
template <int dim, int spacedim>
unsigned int
AgglomerationAccessor<dim, spacedim>::n_agglomerated_faces() const
{
  const auto  &agglomeration = get_agglomerate();
  unsigned int n_neighbors   = 0;
  for (const auto &cell : agglomeration)
    n_neighbors += n_agglomerated_faces_per_cell(cell);
  return n_neighbors;
}
 
 
 
template <int dim, int spacedim>
inline bool
AgglomerationAccessor<dim, spacedim>::at_boundary(const unsigned int f) const
{
  if (master_cell->is_ghost())
    {
      const unsigned int sender_rank = master_cell->subdomain_id();
      // #ifdef AGGLO_DEBUG
      // const auto &       map_recv =
      //   handler->recv_bdary_info.at(sender_rank).at(present_id);
      //       std::cout << "size of map " << map_recv.size() << std::endl;
      //       for (const auto &[face_idx, boundary_or_not] : map_recv)
      //         {
      //           // std::cout << "id: " << face_idx << std::endl;
      //           std::cout << "face: " << face_idx << std::endl;
      //           std::cout << "boundary or not: " << boundary_or_not <<
      //           std::endl;
      //         }
 
      //       std::cout << "il present_id e': " << present_id << "and size "
      //                 << map_recv.size() << std::endl;
      // #endif
 
      return handler->recv_bdary_info.at(sender_rank).at(present_id).at(f);
    }
  else
    {
      Assert(!handler->is_slave_cell(master_cell),
             ExcMessage(
               "This function should not be called for a slave cell."));
 
 
      typename DoFHandler<dim, spacedim>::active_cell_iterator cell_dh(
        *master_cell, &(handler->agglo_dh));
      return handler->at_boundary(cell_dh, f);
    }
}
 
 
 
template <int dim, int spacedim>
inline bool
AgglomerationAccessor<dim, spacedim>::is_locally_owned() const
{
  return master_cell->is_locally_owned();
}
 
 
 
template <int dim, int spacedim>
inline CellId
AgglomerationAccessor<dim, spacedim>::id() const
{
  return present_id;
}
 
 
 
template <int dim, int spacedim>
inline types::subdomain_id
AgglomerationAccessor<dim, spacedim>::subdomain_id() const
{
  return present_subdomain_id;
}
 
template <int dim, int spacedim>
inline const std::vector<types::global_cell_index> &
AgglomerationAccessor<dim, spacedim>::children() const
{
  Assert(!handler->parent_child_info.empty(), ExcInternalError());
  return handler->parent_child_info.at(
    {present_index, handler->present_extraction_level});
}
 
template <int dim, int spacedim>
inline const FiniteElement<dim, spacedim> &
AgglomerationAccessor<dim, spacedim>::get_fe() const
{
  typename DoFHandler<dim, spacedim>::active_cell_iterator
    master_cell_as_dof_handler_iterator =
      master_cell->as_dof_handler_iterator(handler->agglo_dh);
  return master_cell_as_dof_handler_iterator->get_fe();
}
 
template <int dim, int spacedim>
inline void
AgglomerationAccessor<dim, spacedim>::set_active_fe_index(
  const types::fe_index index) const
{
  Assert(!handler->is_slave_cell(master_cell),
         ExcMessage("The present function cannot be called for slave cells."));
  typename DoFHandler<dim, spacedim>::active_cell_iterator
    master_cell_as_dof_handler_iterator =
      master_cell->as_dof_handler_iterator(handler->agglo_dh);
  master_cell_as_dof_handler_iterator->set_active_fe_index(index);
}
 
template <int dim, int spacedim>
inline types::fe_index
AgglomerationAccessor<dim, spacedim>::active_fe_index() const
{
  typename DoFHandler<dim, spacedim>::active_cell_iterator
    master_cell_as_dof_handler_iterator =
      master_cell->as_dof_handler_iterator(handler->agglo_dh);
  return master_cell_as_dof_handler_iterator->active_fe_index();
}
 
#endif