dof_handler.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2005 - 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__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 hp
  {
    class DoFLevel;

    namespace DoFHandler
    {
      struct Implementation;
    }
  }
}

namespace internal
{
  namespace DoFAccessor
  {
    struct Implementation;
  }

  namespace DoFCellAccessor
  {
    struct Implementation;
  }
}



namespace hp
{

  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;

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

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

    virtual ~DoFHandler ();

    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 hp::FECollection<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 (ExcInvalidTriangulation);
    DeclException0 (ExcNoFESelected);
    DeclException0 (ExcRenumberingIncomplete);
    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.");

  protected:

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

    SmartPointer<const hp::FECollection<dim,spacedim>,hp::DoFHandler<dim,spacedim> > finite_elements;

  private:

    DoFHandler (const DoFHandler &);

    DoFHandler &operator = (const DoFHandler &);

    class MGVertexDoFs
    {
    public:
      MGVertexDoFs ();
      ~MGVertexDoFs ();
      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);
    };

    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
    compute_vertex_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const;

    void
    compute_line_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const;

    void
    compute_quad_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const;

    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
                                 ::internal::int2type<0>);

    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
                                 ::internal::int2type<1>);

    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
                                 ::internal::int2type<2>);

    void renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
                                 ::internal::int2type<3>);

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

    ::internal::hp::DoFIndicesOnFaces<dim> *faces;

    ::internal::DoFHandler::NumberCache number_cache;

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

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

    std::vector<MGVertexDoFs> mg_vertex_dofs;  // we should really remove this field!

    std::vector<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::DoFAccessor::Implementation;
    friend struct ::internal::DoFCellAccessor::Implementation;

    template <int> friend class ::internal::hp::DoFIndicesOnFacesOrEdges;
    friend struct ::internal::hp::DoFHandler::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 <int dim, int spacedim>
  template <class Archive>
  void DoFHandler<dim, spacedim>::save(Archive &ar, unsigned int) const
  {
    ar &vertex_dofs;
    ar &vertex_dofs_offsets;
    ar &number_cache;
    ar &levels;
    ar &faces;
    ar &has_children;

    // write out the number of triangulation cells and later check during
    // loading that this number is indeed correct;
    unsigned int n_cells = tria->n_cells();

    ar &n_cells;
  }

  template <int dim, int spacedim>
  template <class Archive>
  void DoFHandler<dim, spacedim>::load(Archive &ar, unsigned int)
  {
    ar &vertex_dofs;
    ar &vertex_dofs_offsets;
    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];
    for (unsigned int i = 0; i<has_children.size(); ++i)
      delete has_children[i];
    levels.resize(0);
    has_children.resize(0);
    delete faces;
    faces = 0;

    ar &levels;
    ar &faces;
    ar &has_children;

    // these are the checks that correspond to the last block in the save()
    // function
    unsigned int n_cells;

    ar &n_cells;

    AssertThrow(n_cells == tria->n_cells(),
                ExcMessage("The object being loaded into does not match the triangulation "
                           "that has been stored previously."));
  }

  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
  {
    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>
  inline
  const hp::FECollection<dim,spacedim> &
  DoFHandler<dim,spacedim>::get_fe () const
  {
    Assert (finite_elements != 0,
            ExcMessage ("No finite element collection is associated with "
                        "this DoFHandler"));
    return *finite_elements;
  }



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



  template<int dim, int spacedim>
  inline
  const Triangulation<dim,spacedim> &
  DoFHandler<dim,spacedim>::get_triangulation () const
  {
    return *tria;
  }



  template<int dim, int spacedim>
  inline
  DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs()
  {
    Assert (false, ExcNotImplemented ());
  }

  template<int dim, int spacedim>
  inline
  DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs()
  {
    Assert (false, ExcNotImplemented ());
  }

  template<int dim, int spacedim>
  inline
  types::global_dof_index DoFHandler<dim, spacedim>::MGVertexDoFs::get_index (const unsigned int,
      const unsigned int) const
  {
    Assert (false, ExcNotImplemented ());
    return invalid_dof_index;
  }

  template<int dim, int spacedim>
  inline
  void DoFHandler<dim, spacedim>::MGVertexDoFs::set_index (const unsigned int,
                                                           const unsigned int,
                                                           types::global_dof_index)
  {
    Assert (false, ExcNotImplemented ());
  }


#endif

}

DEAL_II_NAMESPACE_CLOSE

#endif