dof_handler.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2020 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_hp_dof_handler_h
#define dealii_hp_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/iterator_range.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/template_constraints.h>

#include <deal.II/dofs/deprecated_function_map.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/dofs/dof_iterator_selector.h>
#include <deal.II/dofs/number_cache.h>

#include <deal.II/hp/dof_faces.h>
#include <deal.II/hp/dof_level.h>
#include <deal.II/hp/fe_collection.h>

#include <map>
#include <memory>
#include <set>
#include <vector>

DEAL_II_NAMESPACE_OPEN

// Forward declarations
#ifndef DOXYGEN
template <int dim, int spacedim>
class Triangulation;

namespace parallel
{
  namespace distributed
  {
    template <int dim, int spacedim, typename VectorType>
    class CellDataTransfer;
  }
} // namespace parallel

namespace internal
{
  namespace DoFHandlerImplementation
  {
    struct Implementation;

    namespace Policy
    {
      template <int dim, int spacedim>
      class PolicyBase;
      struct Implementation;
    } // namespace Policy
  }   // namespace DoFHandlerImplementation

  namespace hp
  {
    class DoFLevel;

    namespace DoFHandlerImplementation
    {
      struct Implementation;
    }
  } // namespace hp
} // namespace internal

namespace internal
{
  namespace DoFAccessorImplementation
  {
    struct Implementation;
  }

  namespace DoFCellAccessorImplementation
  {
    struct Implementation;
  }
} // namespace internal
#endif


namespace hp
{
  template <int dim, int spacedim = dim>
  class DoFHandler : public Subscriptor
  {
    using ActiveSelector = ::internal::DoFHandlerImplementation::
      Iterators<DoFHandler<dim, spacedim>, false>;
    using LevelSelector = ::internal::DoFHandlerImplementation::
      Iterators<DoFHandler<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;

#ifndef _MSC_VER
    using active_cell_iterator = typename ActiveSelector::active_cell_iterator;
#else
    using active_cell_iterator = TriaActiveIterator<
      ::DoFCellAccessor<DoFHandler<dim, spacedim>, false>>;
#endif

    using level_cell_iterator = typename LevelSelector::cell_iterator;

#ifndef _MSC_VER
    using cell_iterator = typename ActiveSelector::cell_iterator;
#else
    using cell_iterator =
      TriaIterator<::DoFCellAccessor<DoFHandler<dim, spacedim>, false>>;
#endif

    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_face_iterator = typename LevelSelector::face_iterator;

    static const unsigned int dimension = dim;

    static const unsigned int space_dimension = spacedim;

    static const bool is_hp_dof_handler = true;

    static const unsigned int default_fe_index = numbers::invalid_unsigned_int;


    DoFHandler();

    DoFHandler(const Triangulation<dim, spacedim> &tria);

    DoFHandler(const DoFHandler &) = delete;

    virtual ~DoFHandler() override;

    DoFHandler &
    operator=(const DoFHandler &) = delete;

    void
    initialize(const Triangulation<dim, spacedim> &   tria,
               const hp::FECollection<dim, spacedim> &fe);

    virtual void
    set_fe(const hp::FECollection<dim, spacedim> &fe);

    virtual void
    distribute_dofs(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;

    virtual void
    clear();

    void
    renumber_dofs(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;

    IteratorRange<cell_iterator>
    cell_iterators() const;

    IteratorRange<active_cell_iterator>
    active_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;

    /*
     * @}
     */

    /*---------------------------------------*/


    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;

    types::global_dof_index
    n_locally_owned_dofs() const;

    const IndexSet &
    locally_owned_dofs() const;

    DEAL_II_DEPRECATED const std::vector<IndexSet> &
                             locally_owned_dofs_per_processor() const;

    DEAL_II_DEPRECATED const std::vector<types::global_dof_index> &
                             n_locally_owned_dofs_per_processor() const;

    const IndexSet &
    locally_owned_mg_dofs(const unsigned int level) const;

    DEAL_II_DEPRECATED const std::vector<IndexSet> &
                             locally_owned_mg_dofs_per_processor(const unsigned int level) const;

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

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

    const Triangulation<dim, spacedim> &
    get_triangulation() const;

    virtual std::size_t
    memory_consumption() const;

    void
    prepare_for_serialization_of_active_fe_indices();

    void
    deserialize_active_fe_indices();

    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(ExcGridsDoNotMatch);
    DeclException0(ExcInvalidBoundaryIndicator);
    DeclException1(ExcMatrixHasWrongSize,
                   int,
                   << "The matrix has the wrong dimension " << arg1);
    DeclException0(ExcFunctionNotUseful);
    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");
    DeclException1(ExcInvalidLevel,
                   int,
                   << "The given level " << arg1
                   << " is not in the valid range!");
    DeclException0(ExcFacesHaveNoLevel);
    DeclException1(ExcEmptyLevel,
                   int,
                   << "You tried to do something on level " << arg1
                   << ", but this level is empty.");

  private:
    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;


    void
    setup_policy_and_listeners();

    void
    clear_space();

    template <int structdim>
    types::global_dof_index
    get_dof_index(const unsigned int obj_level,
                  const unsigned int obj_index,
                  const unsigned int fe_index,
                  const unsigned int local_index) const;

    template <int structdim>
    void
    set_dof_index(const unsigned int            obj_level,
                  const unsigned int            obj_index,
                  const unsigned int            fe_index,
                  const unsigned int            local_index,
                  const types::global_dof_index global_index) const;

    void
    create_active_fe_table();

    void
    pre_refinement_action();

    void
    post_refinement_action();

    void
    pre_active_fe_index_transfer();

    void
    pre_distributed_active_fe_index_transfer();

    void
    post_active_fe_index_transfer();

    void
    post_distributed_active_fe_index_transfer();

    void
    post_distributed_serialization_of_active_fe_indices();

    std::vector<std::unique_ptr<::internal::hp::DoFLevel>> levels;

    std::unique_ptr<::internal::hp::DoFIndicesOnFaces<dim>> faces;

    ::internal::DoFHandlerImplementation::NumberCache number_cache;

    std::vector<::internal::DoFHandlerImplementation::NumberCache>
      mg_number_cache;

    std::vector<types::global_dof_index> vertex_dofs;

    std::vector<unsigned int> vertex_dof_offsets;

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

    std::unique_ptr<ActiveFEIndexTransfer> active_fe_index_transfer;

    std::vector<boost::signals2::connection> tria_listeners;

    // Make accessor objects friends.
    template <int, class, bool>
    friend class ::DoFAccessor;
    template <class, 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.
    template <int>
    friend class ::internal::hp::DoFIndicesOnFacesOrEdges;
    friend struct ::internal::hp::DoFHandlerImplementation::
      Implementation;
    friend struct ::internal::DoFHandlerImplementation::Policy::
      Implementation;
  };



#ifndef DOXYGEN


  /* ----------------------- Inline functions ----------------------------------
   */


  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
  {
    // 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);
  }



  template <>
  inline types::global_dof_index
  DoFHandler<2, 3>::n_boundary_dofs() const
  {
    Assert(false, ExcNotImplemented());
    return 0;
  }



  template <>
  template <typename number>
  inline types::global_dof_index
  DoFHandler<2, 3>::n_boundary_dofs(
    const std::map<types::boundary_id, const Function<3, number> *> &) const
  {
    Assert(false, ExcNotImplemented());
    return 0;
  }



  template <>
  inline types::global_dof_index
  DoFHandler<2, 3>::n_boundary_dofs(const std::set<types::boundary_id> &) const
  {
    Assert(false, ExcNotImplemented());
    return 0;
  }
}


namespace internal
{
  template <int dim, int spacedim>
  std::string
  policy_to_string(const ::internal::DoFHandlerImplementation::Policy::
                     PolicyBase<dim, spacedim> &policy);
} // namespace internal


namespace hp
{
  template <int dim, int spacedim>
  template <class Archive>
  void
  DoFHandler<dim, spacedim>::save(Archive &ar, unsigned int) const
  {
    ar &vertex_dofs;
    ar &vertex_dof_offsets;
    ar &number_cache;
    ar &mg_number_cache;

    // some versions of gcc have trouble with loading vectors of
    // std::unique_ptr objects because std::unique_ptr does not
    // have a copy constructor. do it one level at a time
    const unsigned int n_levels = levels.size();
    ar &               n_levels;
    for (unsigned int i = 0; i < n_levels; ++i)
      ar &levels[i];

    // boost dereferences a nullptr when serializing a nullptr
    // at least up to 1.65.1. This causes problems with clang-5.
    // Therefore, work around it.
    bool faces_is_nullptr = (faces.get() == nullptr);
    ar & faces_is_nullptr;
    if (!faces_is_nullptr)
      ar &faces;

    // 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 = tria->n_cells();
    std::string policy_name    = ::internal::policy_to_string(*policy);

    ar &n_cells &policy_name;
  }



  template <int dim, int spacedim>
  template <class Archive>
  void
  DoFHandler<dim, spacedim>::load(Archive &ar, unsigned int)
  {
    ar &vertex_dofs;
    ar &vertex_dof_offsets;
    ar &number_cache;
    ar &mg_number_cache;

    // boost::serialization can restore pointers just fine, but if the
    // pointer object still points to something useful, that object is not
    // destroyed and we end up with a memory leak. consequently, first delete
    // previous content before re-loading stuff
    levels.clear();
    faces.reset();

    // some versions of gcc have trouble with loading vectors of
    // std::unique_ptr objects because std::unique_ptr does not
    // have a copy constructor. do it one level at a time
    unsigned int size;
    ar &         size;
    levels.resize(size);
    for (unsigned int i = 0; i < size; ++i)
      {
        std::unique_ptr<::internal::hp::DoFLevel> level;
        ar &                                            level;
        levels[i] = std::move(level);
      }

    // Workaround for nullptr, see in save().
    bool faces_is_nullptr = true;
    ar & faces_is_nullptr;
    if (!faces_is_nullptr)
      ar &faces;

    // 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 == 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(*policy),
                ExcMessage(
                  "The policy currently associated with this DoFHandler (" +
                  ::internal::policy_to_string(*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>::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) const
  {
    Assert(false, ExcNotImplemented());
    return numbers::invalid_dof_index;
  }



  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 std::vector<types::global_dof_index> &
  DoFHandler<dim, spacedim>::n_locally_owned_dofs_per_processor() const
  {
    if (number_cache.n_locally_owned_dofs_per_processor.empty() &&
        number_cache.n_global_dofs > 0)
      {
        MPI_Comm comm;

        const parallel::TriangulationBase<dim, spacedim> *tr =
          (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
            &this->get_triangulation()));
        if (tr != nullptr)
          comm = tr->get_communicator();
        else
          comm = MPI_COMM_SELF;

        const_cast<::internal::DoFHandlerImplementation::NumberCache &>(
          number_cache)
          .n_locally_owned_dofs_per_processor =
          number_cache.get_n_locally_owned_dofs_per_processor(comm);
      }
    return number_cache.n_locally_owned_dofs_per_processor;
  }



  template <int dim, int spacedim>
  const std::vector<IndexSet> &
  DoFHandler<dim, spacedim>::locally_owned_dofs_per_processor() const
  {
    if (number_cache.locally_owned_dofs_per_processor.empty() &&
        number_cache.n_global_dofs > 0)
      {
        MPI_Comm comm;

        const parallel::TriangulationBase<dim, spacedim> *tr =
          (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
            &this->get_triangulation()));
        if (tr != nullptr)
          comm = tr->get_communicator();
        else
          comm = MPI_COMM_SELF;

        const_cast<::internal::DoFHandlerImplementation::NumberCache &>(
          number_cache)
          .locally_owned_dofs_per_processor =
          number_cache.get_locally_owned_dofs_per_processor(comm);
      }
    return number_cache.locally_owned_dofs_per_processor;
  }



  template <int dim, int spacedim>
  const IndexSet &
  DoFHandler<dim, spacedim>::locally_owned_mg_dofs(
    const unsigned int level) const
  {
    Assert(false, ExcNotImplemented());
    (void)level;
    Assert(level < this->get_triangulation().n_global_levels(),
           ExcMessage("The given level index exceeds the number of levels "
                      "present in the triangulation"));
    return mg_number_cache[0].locally_owned_dofs;
  }



  template <int dim, int spacedim>
  const std::vector<IndexSet> &
  DoFHandler<dim, spacedim>::locally_owned_mg_dofs_per_processor(
    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()?"));
    if (mg_number_cache[level].locally_owned_dofs_per_processor.empty() &&
        mg_number_cache[level].n_global_dofs > 0)
      {
        MPI_Comm comm;

        const parallel::TriangulationBase<dim, spacedim> *tr =
          (dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
            &this->get_triangulation()));
        if (tr != nullptr)
          comm = tr->get_communicator();
        else
          comm = MPI_COMM_SELF;

        const_cast<::internal::DoFHandlerImplementation::NumberCache &>(
          mg_number_cache[level])
          .locally_owned_dofs_per_processor =
          mg_number_cache[level].get_locally_owned_dofs_per_processor(comm);
      }
    return mg_number_cache[level].locally_owned_dofs_per_processor;
  }



  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
  {
    Assert(fe_collection.size() > 0,
           ExcMessage("No finite element collection is associated with "
                      "this DoFHandler"));
    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;
  }

#endif

} // namespace hp

DEAL_II_NAMESPACE_CLOSE

#endif