dof_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 2018 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 at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#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/template_constraints.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/iterator_range.h>
#include <deal.II/dofs/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/fe_collection.h>
#include <deal.II/hp/dof_faces.h>
#include <deal.II/hp/dof_level.h>

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

DEAL_II_NAMESPACE_OPEN

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 hp
  {
    class DoFLevel;

    namespace DoFHandlerImplementation
    {
      struct Implementation;
    }
  }
}

namespace internal
{
  namespace DoFAccessorImplementation
  {
    struct Implementation;
  }

  namespace DoFCellAccessorImplementation
  {
    struct Implementation;
  }
}



namespace hp
{

  template <int dim, int spacedim=dim>
  class DoFHandler : public Subscriptor
  {
    typedef ::internal::DoFHandlerImplementation::Iterators<DoFHandler<dim,spacedim>, false> ActiveSelector;
    typedef ::internal::DoFHandlerImplementation::Iterators<DoFHandler<dim,spacedim>, true> LevelSelector;
  public:
    typedef typename ActiveSelector::CellAccessor         cell_accessor;
    typedef typename ActiveSelector::FaceAccessor         face_accessor;

    typedef typename ActiveSelector::line_iterator        line_iterator;
    typedef typename ActiveSelector::active_line_iterator active_line_iterator;

    typedef typename ActiveSelector::quad_iterator        quad_iterator;
    typedef typename ActiveSelector::active_quad_iterator active_quad_iterator;

    typedef typename ActiveSelector::hex_iterator         hex_iterator;
    typedef typename ActiveSelector::active_hex_iterator  active_hex_iterator;

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

    typedef typename LevelSelector::cell_iterator         level_cell_iterator;

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

    typedef typename ActiveSelector::face_iterator        face_iterator;

    typedef typename ActiveSelector::active_face_iterator active_face_iterator;

    typedef typename LevelSelector::CellAccessor          level_cell_accessor;
    typedef typename LevelSelector::FaceAccessor          level_face_accessor;

    typedef typename LevelSelector::face_iterator         level_face_iterator;

    static const unsigned int dimension = dim;

    static const unsigned int space_dimension = spacedim;

    DEAL_II_DEPRECATED
    static const types::global_dof_index invalid_dof_index = numbers::invalid_dof_index;

    static const unsigned int default_fe_index = numbers::invalid_unsigned_int;


    DoFHandler ();

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

    DoFHandler (const DoFHandler &) = delete;

    virtual ~DoFHandler ();

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

    void initialize(const Triangulation<dim,spacedim> &tria,
                    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;

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

    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;

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

    DEAL_II_DEPRECATED
    const hp::FECollection<dim,spacedim> &get_fe () 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;

    template <class Archive>
    void save(Archive &ar, const unsigned int version) const;

    template <class Archive>
    void load(Archive &ar, const unsigned int version);

    BOOST_SERIALIZATION_SPLIT_MEMBER()

    
    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 ();

    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;

    std::vector<std::unique_ptr<std::vector<bool> > > has_children;

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

    template <int, class, bool> friend class ::DoFAccessor;
    template <class, bool> friend class ::DoFCellAccessor;
    friend struct ::internal::DoFAccessorImplementation::Implementation;
    friend struct ::internal::DoFCellAccessorImplementation::Implementation;

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


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


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;

    // the same issue as above
    const unsigned int n_has_children = has_children.size();
    ar &n_has_children;
    for (unsigned int i = 0; i < n_has_children; ++i)
      ar &has_children[i];

    // 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 ();
    has_children.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;

    // the same issue as above
    ar &size;
    has_children.resize(size);
    for (unsigned int i = 0; i < size; ++i)
      {
        std::unique_ptr<std::vector<bool> > has_children_on_level;
        ar &has_children_on_level;
        has_children[i] = std::move(has_children_on_level);
      }

    // 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
  {
    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
  {
    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("invalid level in locally_owned_mg_dofs"));
    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(false, ExcNotImplemented());
    (void)level;
    Assert(level < this->get_triangulation().n_global_levels(),
           ExcMessage("invalid level in locally_owned_mg_dofs_per_processor"));
    return mg_number_cache[0].locally_owned_dofs_per_processor;
  }



  template <int dim, int spacedim>
  inline
  const hp::FECollection<dim,spacedim> &
  DoFHandler<dim,spacedim>::get_fe () 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 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

}

DEAL_II_NAMESPACE_CLOSE

#endif