dof_handler.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2022 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_dof_handler_h
#define dealii_dof_handler_h
 
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/exceptions.h>
#include <deal.II/base/function.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/iterator_range.h>
#include <deal.II/base/smartpointer.h>
 
#include <deal.II/distributed/tria_base.h>
 
#include <deal.II/dofs/block_info.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_faces.h>
#include <deal.II/dofs/dof_iterator_selector.h>
#include <deal.II/dofs/dof_levels.h>
#include <deal.II/dofs/number_cache.h>
 
#include <deal.II/hp/fe_collection.h>
 
#include <boost/serialization/split_member.hpp>
 
#include <map>
#include <memory>
#include <set>
#include <vector>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <int dim, int spacedim>
class FiniteElement;
template <int dim, int spacedim>
class Triangulation;
 
namespace internal
{
  namespace DoFHandlerImplementation
  {
    struct Implementation;
 
    namespace Policy
    {
      template <int dim, int spacedim>
      class PolicyBase;
      struct Implementation;
    } // namespace Policy
  }   // namespace DoFHandlerImplementation
 
  namespace DoFAccessorImplementation
  {
    struct Implementation;
  }
 
  namespace DoFCellAccessorImplementation
  {
    struct Implementation;
  }
 
  namespace hp
  {
    namespace DoFHandlerImplementation
    {
      struct Implementation;
    }
  } // namespace hp
} // namespace internal
 
namespace parallel
{
  namespace distributed
  {
    template <int dim, int spacedim, typename VectorType>
    class CellDataTransfer;
  }
} // namespace parallel
#endif
 
template <int dim, int spacedim = dim>
class DoFHandler : public Subscriptor
{
  using ActiveSelector =
    ::internal::DoFHandlerImplementation::Iterators<dim, spacedim, false>;
  using LevelSelector =
    ::internal::DoFHandlerImplementation::Iterators<dim, spacedim, true>;
 
public:
  using cell_accessor = typename ActiveSelector::CellAccessor;
 
  using face_accessor = typename ActiveSelector::FaceAccessor;
 
  using line_iterator = typename ActiveSelector::line_iterator;
 
  using active_line_iterator = typename ActiveSelector::active_line_iterator;
 
  using quad_iterator = typename ActiveSelector::quad_iterator;
 
  using active_quad_iterator = typename ActiveSelector::active_quad_iterator;
 
  using hex_iterator = typename ActiveSelector::hex_iterator;
 
  using active_hex_iterator = typename ActiveSelector::active_hex_iterator;
 
  using active_cell_iterator = typename ActiveSelector::active_cell_iterator;
 
  using cell_iterator = typename ActiveSelector::cell_iterator;
 
  using face_iterator = typename ActiveSelector::face_iterator;
 
  using active_face_iterator = typename ActiveSelector::active_face_iterator;
 
  using level_cell_accessor = typename LevelSelector::CellAccessor;
  using level_face_accessor = typename LevelSelector::FaceAccessor;
 
  using level_cell_iterator = typename LevelSelector::cell_iterator;
  using level_face_iterator = typename LevelSelector::face_iterator;
 
 
  static constexpr unsigned int dimension = dim;
 
  static constexpr unsigned int space_dimension = spacedim;
 
  static const unsigned int default_fe_index = 0;
 
  static const unsigned int invalid_fe_index = numbers::invalid_unsigned_int;
 
  using active_fe_index_type = unsigned short int;
 
  using offset_type = unsigned int;
 
  static const active_fe_index_type invalid_active_fe_index =
    static_cast<active_fe_index_type>(-1);
 
  DoFHandler();
 
  explicit DoFHandler(const Triangulation<dim, spacedim> &tria);
 
  DoFHandler(const DoFHandler &) = delete;
 
  virtual ~DoFHandler() override;
 
  DoFHandler &
  operator=(const DoFHandler &) = delete;
 
  DEAL_II_DEPRECATED
  void
  initialize(const Triangulation<dim, spacedim> &tria,
             const FiniteElement<dim, spacedim> &fe);
 
  DEAL_II_DEPRECATED
  void
  initialize(const Triangulation<dim, spacedim> &   tria,
             const hp::FECollection<dim, spacedim> &fe);
 
  DEAL_II_DEPRECATED
  void
  set_fe(const FiniteElement<dim, spacedim> &fe);
 
  DEAL_II_DEPRECATED
  void
  set_fe(const hp::FECollection<dim, spacedim> &fe);
 
  void
  set_active_fe_indices(const std::vector<unsigned int> &active_fe_indices);
 
  void
  get_active_fe_indices(std::vector<unsigned int> &active_fe_indices) const;
 
  void
  reinit(const Triangulation<dim, spacedim> &tria);
 
  void
  distribute_dofs(const FiniteElement<dim, spacedim> &fe);
 
  void
  distribute_dofs(const hp::FECollection<dim, spacedim> &fe);
 
  void
  distribute_mg_dofs();
 
  bool
  has_hp_capabilities() const;
 
  bool
  has_level_dofs() const;
 
  bool
  has_active_dofs() const;
 
  void
  initialize_local_block_info();
 
  void
  clear();
 
  void
  renumber_dofs(const std::vector<types::global_dof_index> &new_numbers);
 
  void
  renumber_dofs(const unsigned int                          level,
                const std::vector<types::global_dof_index> &new_numbers);
 
  unsigned int
  max_couplings_between_dofs() const;
 
  unsigned int
  max_couplings_between_boundary_dofs() const;
 
  /*--------------------------------------*/
 
  /*
   * @{
   */
 
  cell_iterator
  begin(const unsigned int level = 0) const;
 
  active_cell_iterator
  begin_active(const unsigned int level = 0) const;
 
  cell_iterator
  end() const;
 
  cell_iterator
  end(const unsigned int level) const;
 
  active_cell_iterator
  end_active(const unsigned int level) const;
 
  level_cell_iterator
  begin_mg(const unsigned int level = 0) const;
 
  level_cell_iterator
  end_mg(const unsigned int level) const;
 
  level_cell_iterator
  end_mg() const;
 
  IteratorRange<cell_iterator>
  cell_iterators() const;
 
  IteratorRange<active_cell_iterator>
  active_cell_iterators() const;
 
  IteratorRange<level_cell_iterator>
  mg_cell_iterators() const;
 
  IteratorRange<cell_iterator>
  cell_iterators_on_level(const unsigned int level) const;
 
  IteratorRange<active_cell_iterator>
  active_cell_iterators_on_level(const unsigned int level) const;
 
  IteratorRange<level_cell_iterator>
  mg_cell_iterators_on_level(const unsigned int level) const;
 
  /*
   * @}
   */
 
 
  /*---------------------------------------*/
 
 
  types::global_dof_index
  n_dofs() const;
 
  types::global_dof_index
  n_dofs(const unsigned int level) const;
 
  types::global_dof_index
  n_boundary_dofs() const;
 
  template <typename number>
  types::global_dof_index
  n_boundary_dofs(
    const std::map<types::boundary_id, const Function<spacedim, number> *>
      &boundary_ids) const;
 
  types::global_dof_index
  n_boundary_dofs(const std::set<types::boundary_id> &boundary_ids) const;
 
  const BlockInfo &
  block_info() const;
 
  types::global_dof_index
  n_locally_owned_dofs() const;
 
  const IndexSet &
  locally_owned_dofs() const;
 
  const IndexSet &
  locally_owned_mg_dofs(const unsigned int level) const;
 
  const FiniteElement<dim, spacedim> &
  get_fe(const unsigned int index = 0) const;
 
  const hp::FECollection<dim, spacedim> &
  get_fe_collection() const;
 
  const Triangulation<dim, spacedim> &
  get_triangulation() const;
 
  MPI_Comm
  get_communicator() const;
 
  void
  prepare_for_serialization_of_active_fe_indices();
 
  void
  deserialize_active_fe_indices();
 
  virtual std::size_t
  memory_consumption() const;
 
  template <class Archive>
  void
  save(Archive &ar, const unsigned int version) const;
 
  template <class Archive>
  void
  load(Archive &ar, const unsigned int version);
 
#ifdef DOXYGEN
  template <class Archive>
  void
  serialize(Archive &archive, const unsigned int version);
#else
  // This macro defines the serialize() method that is compatible with
  // the templated save() and load() method that have been implemented.
  BOOST_SERIALIZATION_SPLIT_MEMBER()
#endif
 
  DeclException0(ExcNoFESelected);
  DeclException0(ExcInvalidBoundaryIndicator);
  DeclException1(ExcInvalidLevel,
                 int,
                 << "The given level " << arg1
                 << " is not in the valid range!");
  DeclException1(ExcNewNumbersNotConsecutive,
                 types::global_dof_index,
                 << "The given list of new dof indices is not consecutive: "
                 << "the index " << arg1 << " does not exist.");
  DeclException2(ExcInvalidFEIndex,
                 int,
                 int,
                 << "The mesh contains a cell with an active FE index of "
                 << arg1 << ", but the finite element collection only has "
                 << arg2 << " elements");
 
  DeclExceptionMsg(ExcOnlyAvailableWithHP,
                   "The current function doesn't make sense when used with a "
                   "DoFHandler without hp-capabilities.");
 
  DeclExceptionMsg(ExcNotImplementedWithHP,
                   "The current function has not yet been implemented for a "
                   "DoFHandler with hp-capabilities.");
 
private:
  class MGVertexDoFs
  {
  public:
    MGVertexDoFs();
 
    void
    init(const unsigned int coarsest_level,
         const unsigned int finest_level,
         const unsigned int dofs_per_vertex);
 
    unsigned int
    get_coarsest_level() const;
 
    unsigned int
    get_finest_level() const;
 
    types::global_dof_index &
    access_index(const unsigned int level,
                 const unsigned int dof_number,
                 const unsigned int dofs_per_vertex);
 
  private:
    unsigned int coarsest_level;
 
    unsigned int finest_level;
 
    std::unique_ptr<types::global_dof_index[]> indices;
  };
 
  struct ActiveFEIndexTransfer
  {
    std::map<const cell_iterator, const unsigned int> persisting_cells_fe_index;
 
    std::map<const cell_iterator, const unsigned int> refined_cells_fe_index;
 
    std::map<const cell_iterator, const unsigned int> coarsened_cells_fe_index;
 
    std::vector<unsigned int> active_fe_indices;
 
    std::unique_ptr<
      parallel::distributed::
        CellDataTransfer<dim, spacedim, std::vector<unsigned int>>>
      cell_data_transfer;
  };
 
  BlockInfo block_info_object;
 
  bool hp_capability_enabled;
 
  SmartPointer<const Triangulation<dim, spacedim>, DoFHandler<dim, spacedim>>
    tria;
 
  hp::FECollection<dim, spacedim> fe_collection;
 
  std::unique_ptr<::internal::DoFHandlerImplementation::Policy::
                    PolicyBase<dim, spacedim>>
    policy;
 
  ::internal::DoFHandlerImplementation::NumberCache number_cache;
 
  std::vector<::internal::DoFHandlerImplementation::NumberCache>
    mg_number_cache;
 
  mutable std::vector<std::vector<types::global_dof_index>>
    cell_dof_cache_indices;
 
  mutable std::vector<std::vector<offset_type>> cell_dof_cache_ptr;
 
  mutable std::vector<std::array<std::vector<types::global_dof_index>, dim + 1>>
    object_dof_indices;
 
  mutable std::vector<std::array<std::vector<offset_type>, dim + 1>>
    object_dof_ptr;
 
  mutable std::array<std::vector<active_fe_index_type>, dim + 1>
    hp_object_fe_indices;
 
  mutable std::array<std::vector<offset_type>, dim + 1> hp_object_fe_ptr;
 
  mutable std::vector<std::vector<active_fe_index_type>>
    hp_cell_active_fe_indices;
 
  mutable std::vector<std::vector<active_fe_index_type>>
    hp_cell_future_fe_indices;
 
  std::vector<MGVertexDoFs> mg_vertex_dofs;
 
  std::vector<
    std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel<dim>>>
    mg_levels;
 
  std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces<dim>>
    mg_faces;
 
  std::unique_ptr<ActiveFEIndexTransfer> active_fe_index_transfer;
 
  std::vector<boost::signals2::connection> tria_listeners;
 
  std::vector<boost::signals2::connection> tria_listeners_for_transfer;
 
  void
  clear_space();
 
  void
  clear_mg_space();
 
  void
  setup_policy();
 
  void
  connect_to_triangulation_signals();
 
  void
  create_active_fe_table();
 
  void
  update_active_fe_table();
 
  void
  pre_transfer_action();
 
  void
  post_transfer_action();
 
  void
  pre_distributed_transfer_action();
 
  void
  post_distributed_transfer_action();
 
 
  // Make accessor objects friends.
  template <int, int, int, bool>
  friend class ::DoFAccessor;
  template <int, int, bool>
  friend class ::DoFCellAccessor;
  friend struct ::internal::DoFAccessorImplementation::Implementation;
  friend struct ::internal::DoFCellAccessorImplementation::Implementation;
 
  // Likewise for DoFLevel objects since they need to access the vertex dofs
  // in the functions that set and retrieve vertex dof indices.
  friend struct ::internal::DoFHandlerImplementation::Implementation;
  friend struct ::internal::hp::DoFHandlerImplementation::Implementation;
  friend struct ::internal::DoFHandlerImplementation::Policy::
    Implementation;
 
  // explicitly check for sensible template arguments, but not on windows
  // because MSVC creates bogus warnings during normal compilation
#ifndef DEAL_II_MSVC
  static_assert(dim <= spacedim,
                "The dimension <dim> of a DoFHandler must be less than or "
                "equal to the space dimension <spacedim> in which it lives.");
#endif
};
 
namespace internal
{
  namespace hp
  {
    namespace DoFHandlerImplementation
    {
      template <int dim, int spacedim>
      void
      communicate_future_fe_indices(DoFHandler<dim, spacedim> &dof_handler);
 
      template <int dim, int spacedim = dim>
      unsigned int
      dominated_future_fe_on_children(
        const typename DoFHandler<dim, spacedim>::cell_iterator &parent);
 
      DeclExceptionMsg(
        ExcNoDominatedFiniteElementOnChildren,
        "No FiniteElement has been found in your FECollection that is "
        "dominated by all children of a cell you are trying to coarsen!");
    } // namespace DoFHandlerImplementation
  }   // namespace hp
} // namespace internal
 
#ifndef DOXYGEN
 
/* ----------------------- Inline functions ----------------------------------
 */
 
 
template <int dim, int spacedim>
inline bool
DoFHandler<dim, spacedim>::has_hp_capabilities() const
{
  return hp_capability_enabled;
}
 
 
 
template <int dim, int spacedim>
inline bool
DoFHandler<dim, spacedim>::has_level_dofs() const
{
  return mg_number_cache.size() > 0;
}
 
 
 
template <int dim, int spacedim>
inline bool
DoFHandler<dim, spacedim>::has_active_dofs() const
{
  return number_cache.n_global_dofs > 0;
}
 
 
 
template <int dim, int spacedim>
inline types::global_dof_index
DoFHandler<dim, spacedim>::n_dofs() const
{
  return number_cache.n_global_dofs;
}
 
 
 
template <int dim, int spacedim>
inline types::global_dof_index
DoFHandler<dim, spacedim>::n_dofs(const unsigned int level) const
{
  Assert(has_level_dofs(),
         ExcMessage(
           "n_dofs(level) can only be called after distribute_mg_dofs()"));
  Assert(level < mg_number_cache.size(), ExcInvalidLevel(level));
  return mg_number_cache[level].n_global_dofs;
}
 
 
 
template <int dim, int spacedim>
types::global_dof_index
DoFHandler<dim, spacedim>::n_locally_owned_dofs() const
{
  return number_cache.n_locally_owned_dofs;
}
 
 
 
template <int dim, int spacedim>
const IndexSet &
DoFHandler<dim, spacedim>::locally_owned_dofs() const
{
  return number_cache.locally_owned_dofs;
}
 
 
 
template <int dim, int spacedim>
const IndexSet &
DoFHandler<dim, spacedim>::locally_owned_mg_dofs(const unsigned int level) const
{
  Assert(level < this->get_triangulation().n_global_levels(),
         ExcMessage("The given level index exceeds the number of levels "
                    "present in the triangulation"));
  Assert(
    mg_number_cache.size() == this->get_triangulation().n_global_levels(),
    ExcMessage(
      "The level dofs are not set up properly! Did you call distribute_mg_dofs()?"));
  return mg_number_cache[level].locally_owned_dofs;
}
 
 
 
template <int dim, int spacedim>
inline const FiniteElement<dim, spacedim> &
DoFHandler<dim, spacedim>::get_fe(const unsigned int number) const
{
  Assert(fe_collection.size() > 0,
         ExcMessage("No finite element collection is associated with "
                    "this DoFHandler"));
  return fe_collection[number];
}
 
 
 
template <int dim, int spacedim>
inline const hp::FECollection<dim, spacedim> &
DoFHandler<dim, spacedim>::get_fe_collection() const
{
  return fe_collection;
}
 
 
 
template <int dim, int spacedim>
inline const Triangulation<dim, spacedim> &
DoFHandler<dim, spacedim>::get_triangulation() const
{
  Assert(tria != nullptr,
         ExcMessage("This DoFHandler object has not been associated "
                    "with a triangulation."));
  return *tria;
}
 
 
 
template <int dim, int spacedim>
inline MPI_Comm
DoFHandler<dim, spacedim>::get_communicator() const
{
  Assert(tria != nullptr,
         ExcMessage("This DoFHandler object has not been associated "
                    "with a triangulation."));
  return tria->get_communicator();
}
 
 
 
template <int dim, int spacedim>
inline const BlockInfo &
DoFHandler<dim, spacedim>::block_info() const
{
  Assert(this->hp_capability_enabled == false, ExcNotImplementedWithHP());
 
  return block_info_object;
}
 
 
 
template <int dim, int spacedim>
template <typename number>
types::global_dof_index
DoFHandler<dim, spacedim>::n_boundary_dofs(
  const std::map<types::boundary_id, const Function<spacedim, number> *>
    &boundary_ids) const
{
  Assert(!(dim == 2 && spacedim == 3) || this->hp_capability_enabled == false,
         ExcNotImplementedWithHP());
 
  // extract the set of boundary ids and forget about the function object
  // pointers
  std::set<types::boundary_id> boundary_ids_only;
  for (typename std::map<types::boundary_id,
                         const Function<spacedim, number> *>::const_iterator p =
         boundary_ids.begin();
       p != boundary_ids.end();
       ++p)
    boundary_ids_only.insert(p->first);
 
  // then just hand everything over to the other function that does the work
  return n_boundary_dofs(boundary_ids_only);
}
 
 
 
namespace internal
{
  template <int dim, int spacedim>
  std::string
  policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
                     PolicyBase<dim, spacedim> &policy);
} // namespace internal
 
 
 
template <int dim, int spacedim>
template <class Archive>
void
DoFHandler<dim, spacedim>::save(Archive &ar, const unsigned int) const
{
  if (this->hp_capability_enabled)
    {
      ar &this->object_dof_indices;
      ar &this->object_dof_ptr;
 
      ar &this->cell_dof_cache_indices;
      ar &this->cell_dof_cache_ptr;
 
      ar &this->hp_cell_active_fe_indices;
      ar &this->hp_cell_future_fe_indices;
 
      ar &hp_object_fe_ptr;
      ar &hp_object_fe_indices;
 
      ar &number_cache;
 
      ar &mg_number_cache;
 
      // write out the number of triangulation cells and later check during
      // loading that this number is indeed correct; same with something that
      // identifies the policy
      const unsigned int n_cells = this->tria->n_cells();
      std::string        policy_name =
        ::internal::policy_to_string(*this->policy);
 
      ar &n_cells &policy_name;
    }
  else
    {
      ar &this->block_info_object;
      ar &number_cache;
 
      ar &this->object_dof_indices;
      ar &this->object_dof_ptr;
 
      ar &this->cell_dof_cache_indices;
      ar &this->cell_dof_cache_ptr;
 
      // write out the number of triangulation cells and later check during
      // loading that this number is indeed correct; same with something that
      // identifies the FE and the policy
      unsigned int n_cells     = this->tria->n_cells();
      std::string  fe_name     = this->get_fe(0).get_name();
      std::string  policy_name = internal::policy_to_string(*this->policy);
 
      ar &n_cells &fe_name &policy_name;
    }
}
 
 
 
template <int dim, int spacedim>
template <class Archive>
void
DoFHandler<dim, spacedim>::load(Archive &ar, const unsigned int)
{
  if (this->hp_capability_enabled)
    {
      ar &this->object_dof_indices;
      ar &this->object_dof_ptr;
 
      ar &this->cell_dof_cache_indices;
      ar &this->cell_dof_cache_ptr;
 
      ar &this->hp_cell_active_fe_indices;
      ar &this->hp_cell_future_fe_indices;
 
      ar &hp_object_fe_ptr;
      ar &hp_object_fe_indices;
 
      ar &number_cache;
 
      ar &mg_number_cache;
 
      // these are the checks that correspond to the last block in the save()
      // function
      unsigned int n_cells;
      std::string  policy_name;
 
      ar &n_cells &policy_name;
 
      AssertThrow(
        n_cells == this->tria->n_cells(),
        ExcMessage(
          "The object being loaded into does not match the triangulation "
          "that has been stored previously."));
      AssertThrow(
        policy_name == ::internal::policy_to_string(*this->policy),
        ExcMessage("The policy currently associated with this DoFHandler (" +
                   ::internal::policy_to_string(*this->policy) +
                   ") does not match the one that was associated with the "
                   "DoFHandler previously stored (" +
                   policy_name + ")."));
    }
  else
    {
      ar &this->block_info_object;
      ar &number_cache;
 
      object_dof_indices.clear();
 
      object_dof_ptr.clear();
 
      ar &this->object_dof_indices;
      ar &this->object_dof_ptr;
 
      ar &this->cell_dof_cache_indices;
      ar &this->cell_dof_cache_ptr;
 
      // these are the checks that correspond to the last block in the save()
      // function
      unsigned int n_cells;
      std::string  fe_name;
      std::string  policy_name;
 
      ar &n_cells &fe_name &policy_name;
 
      AssertThrow(
        n_cells == this->tria->n_cells(),
        ExcMessage(
          "The object being loaded into does not match the triangulation "
          "that has been stored previously."));
      AssertThrow(
        fe_name == this->get_fe(0).get_name(),
        ExcMessage(
          "The finite element associated with this DoFHandler does not match "
          "the one that was associated with the DoFHandler previously stored."));
      AssertThrow(policy_name == internal::policy_to_string(*this->policy),
                  ExcMessage(
                    "The policy currently associated with this DoFHandler (" +
                    internal::policy_to_string(*this->policy) +
                    ") does not match the one that was associated with the "
                    "DoFHandler previously stored (" +
                    policy_name + ")."));
    }
}
 
 
 
template <int dim, int spacedim>
inline types::global_dof_index &
DoFHandler<dim, spacedim>::MGVertexDoFs::access_index(
  const unsigned int level,
  const unsigned int dof_number,
  const unsigned int dofs_per_vertex)
{
  Assert((level >= coarsest_level) && (level <= finest_level),
         ExcInvalidLevel(level));
  return indices[dofs_per_vertex * (level - coarsest_level) + dof_number];
}
 
 
 
extern template class DoFHandler<1, 1>;
extern template class DoFHandler<1, 2>;
extern template class DoFHandler<1, 3>;
extern template class DoFHandler<2, 2>;
extern template class DoFHandler<2, 3>;
extern template class DoFHandler<3, 3>;
 
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif