dof_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 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_dof_handler_h
#define dealii_dof_handler_h



#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/smartpointer.h>
#include <deal.II/base/index_set.h>
#include <deal.II/base/iterator_range.h>
#include <deal.II/dofs/block_info.h>
#include <deal.II/dofs/dof_iterator_selector.h>
#include <deal.II/dofs/number_cache.h>
#include <deal.II/dofs/dof_faces.h>
#include <deal.II/dofs/dof_levels.h>
#include <deal.II/dofs/function_map.h>
#include <deal.II/hp/fe_collection.h>

#include <boost/serialization/split_member.hpp>

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

DEAL_II_NAMESPACE_OPEN

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 DoFAccessorImplementation
  {
    struct Implementation;
  }

  namespace DoFCellAccessorImplementation
  {
    struct Implementation;
  }
}


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;

  typedef typename ActiveSelector::active_cell_iterator active_cell_iterator;

  typedef typename ActiveSelector::cell_iterator        cell_iterator;

  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::cell_iterator         level_cell_iterator;
  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 = 0;

  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 FiniteElement<dim,spacedim> &fe);

  virtual void distribute_dofs (const FiniteElement<dim,spacedim> &fe);

  DEAL_II_DEPRECATED
  virtual void distribute_mg_dofs (const FiniteElement<dim, spacedim> &fe);

  virtual void distribute_mg_dofs ();

  bool has_level_dofs() const;

  bool has_active_dofs() const;

  void initialize_local_block_info();

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


  unsigned int n_locally_owned_dofs() const;

  const IndexSet &locally_owned_dofs() const;

  const IndexSet &locally_owned_mg_dofs(const unsigned int level) 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 std::vector<IndexSet> &
  locally_owned_mg_dofs_per_processor (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;

  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 (ExcGridsDoNotMatch);
  DeclException0 (ExcInvalidBoundaryIndicator);
  DeclException1 (ExcNewNumbersNotConsecutive,
                  types::global_dof_index,
                  << "The given list of new dof indices is not consecutive: "
                  << "the index " << arg1 << " does not exist.");
  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:

  BlockInfo block_info_object;

  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;

  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
    get_index (const unsigned int level,
               const unsigned int dof_number,
               const unsigned int dofs_per_vertex) const;

    void set_index (const unsigned int            level,
                    const unsigned int            dof_number,
                    const unsigned int            dofs_per_vertex,
                    const types::global_dof_index index);

  private:
    unsigned int coarsest_level;

    unsigned int finest_level;

    std::unique_ptr<types::global_dof_index[]> indices;
  };

  void clear_space ();

  void clear_mg_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;

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

  std::vector<MGVertexDoFs> mg_vertex_dofs;

  std::vector<std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel<dim> > > levels;

  std::vector<std::unique_ptr<::internal::DoFHandlerImplementation::DoFLevel<dim> > > mg_levels;

  std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces<dim> > faces;

  std::unique_ptr<::internal::DoFHandlerImplementation::DoFFaces<dim> > mg_faces;

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

  friend struct ::internal::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
};




#ifndef DOXYGEN

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


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>
unsigned int
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("invalid level in locally_owned_mg_dofs"));
  return mg_number_cache[level].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 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("invalid level in locally_owned_mg_dofs_per_processor"));
  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 index) const
{
  (void) index;
  Assert(index == 0, ExcMessage("There is only one FiniteElement stored. The index must be zero!"));
  return get_fe_collection()[0];
}



template <int dim, int spacedim>
inline
const hp::FECollection<dim,spacedim> &
DoFHandler<dim,spacedim>::get_fe_collection () const
{
  Assert(fe_collection.size() > 0,
         ExcMessage("You are trying to access the DoFHandler's FECollection object before it has been initialized."));
  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
const BlockInfo &
DoFHandler<dim,spacedim>::block_info () const
{
  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
{
  // 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);
}



template <int dim, int spacedim>
template <class Archive>
void DoFHandler<dim,spacedim>::save (Archive &ar,
                                     const unsigned int) const
{
  ar &block_info_object;
  ar &vertex_dofs;
  ar &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
  unsigned int n_levels = levels.size();
  ar &n_levels;
  for (unsigned int i = 0; i < levels.size(); ++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 FE and the policy
  unsigned int n_cells = tria->n_cells();
  std::string  fe_name = this->get_fe(0).get_name();
  std::string  policy_name = internal::policy_to_string(*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)
{
  ar &block_info_object;
  ar &vertex_dofs;
  ar &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 < levels.size(); ++i)
    {
      std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<dim>> 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  fe_name;
  std::string  policy_name;

  ar &n_cells &fe_name &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 (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(*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>::MGVertexDoFs::get_index (const unsigned int level,
                                                    const unsigned int dof_number,
                                                    const unsigned int dofs_per_vertex) const
{
  Assert ((level >= coarsest_level) && (level <= finest_level), ExcInvalidLevel (level));
  return indices[dofs_per_vertex * (level - coarsest_level) + dof_number];
}



template <int dim, int spacedim>
inline
void
DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int            level,
                                                    const unsigned int            dof_number,
                                                    const unsigned int            dofs_per_vertex,
                                                    const types::global_dof_index index)
{
  Assert ((level >= coarsest_level) && (level <= finest_level), ExcInvalidLevel (level));
  indices[dofs_per_vertex * (level - coarsest_level) + dof_number] = index;
}


#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif