dof_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 1998 - 2017 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/base/std_cxx11/shared_ptr.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 <boost/serialization/split_member.hpp>

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

DEAL_II_NAMESPACE_OPEN

template <int dim, int spacedim> class FiniteElement;
template <int dim, int spacedim> class Triangulation;

namespace internal
{
  namespace DoFHandler
  {
    struct Implementation;

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

  namespace DoFAccessor
  {
    struct Implementation;
  }

  namespace DoFCellAccessor
  {
    struct Implementation;
  }
}


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

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

  virtual ~DoFHandler ();

  void initialize(const Triangulation<dim,spacedim> &tria,
                  const FiniteElement<dim,spacedim> &fe);

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

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

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

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

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

  const Triangulation<dim,spacedim> &get_tria () const DEAL_II_DEPRECATED;

  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 (ExcRenumberingIncomplete);
  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:
  DoFHandler (const DoFHandler &);

  DoFHandler &operator = (const DoFHandler &);


  BlockInfo block_info_object;

  SmartPointer<const Triangulation<dim,spacedim>,DoFHandler<dim,spacedim> >
  tria;

  SmartPointer<const FiniteElement<dim,spacedim>,DoFHandler<dim,spacedim> >
  selected_fe;

  std_cxx11::shared_ptr<::internal::DoFHandler::Policy::PolicyBase<dim,spacedim> > policy;

  ::internal::DoFHandler::NumberCache number_cache;

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

  class MGVertexDoFs
  {
  public:
    MGVertexDoFs ();

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

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

    DeclException0 (ExcNoMemory);

  private:
    unsigned int coarsest_level;

    unsigned int finest_level;

    types::global_dof_index *indices;

    types::global_dof_index *indices_offset;
  };

  void clear_mg_space ();

  void clear_space ();

  void reserve_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<::internal::DoFHandler::DoFLevel<dim>*> levels;

  std::vector<::internal::DoFHandler::DoFLevel<dim>*> mg_levels;

  ::internal::DoFHandler::DoFFaces<dim> *faces;

  ::internal::DoFHandler::DoFFaces<dim> *mg_faces;

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

  friend struct ::internal::DoFHandler::Implementation;
  friend struct ::internal::DoFHandler::Policy::Implementation;

  // explicitly check for sensible template arguments, but not on windows
  // because MSVC creates bogus warnings during normal compilation
#ifdef DEAL_II_WITH_CXX11
#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
#endif
};




/* -------------- declaration of explicit specializations ------------- */

#ifndef DOXYGEN

template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs () const;
template <> types::global_dof_index DoFHandler<1>::n_boundary_dofs (const std::set<types::boundary_id> &) const;

template <> void DoFHandler<1>::renumber_dofs(unsigned int,const std::vector<types::global_dof_index>  &new_numbers);
template <> void DoFHandler<2>::renumber_dofs(unsigned int,const std::vector<types::global_dof_index>  &new_numbers);
template <> void DoFHandler<3>::renumber_dofs(unsigned int,const std::vector<types::global_dof_index>  &new_numbers);


/* ----------------------- 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
{
  Assert(selected_fe!=0, ExcMessage("You are trying to access the DoFHandler's FiniteElement object before it has been initialized."));
  return *selected_fe;
}



template <int dim, int spacedim>
inline
const Triangulation<dim,spacedim> &
DoFHandler<dim,spacedim>::get_tria () const
{
  Assert(tria != 0, ExcNotInitialized());
  return *tria;
}



template <int dim, int spacedim>
inline
const Triangulation<dim,spacedim> &
DoFHandler<dim,spacedim>::get_triangulation () const
{
  Assert(tria != 0, ExcNotInitialized());
  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::DoFHandler::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;
  ar &levels;
  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 = selected_fe->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
  for (unsigned int i=0; i<levels.size(); ++i)
    delete levels[i];
  levels.resize (0);
  delete faces;
  faces = 0;

  ar &levels;
  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 == selected_fe->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 (std::string ("The policy currently associated with this DoFHandler (")
                           + internal::policy_to_string(*policy)
                           +std::string(") 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
{
  Assert ((level >= coarsest_level) && (level <= finest_level), ExcInvalidLevel (level));
  return indices[indices_offset[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 types::global_dof_index index)
{
  Assert ((level >= coarsest_level) && (level <= finest_level), ExcInvalidLevel (level));
  indices[indices_offset[level - coarsest_level] + dof_number] = index;
}

#endif // DOXYGEN

DEAL_II_NAMESPACE_CLOSE

#endif