libdeal.ii-dev 8.1.0-6ubuntu1 / usr / include / deal.II / dofs / dof_accessor.h

// ---------------------------------------------------------------------
// $Id: dof_accessor.h 30040 2013-07-18 17:06:48Z maier $
//
// Copyright (C) 1998 - 2013 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 __deal2__dof_accessor_h
#define __deal2__dof_accessor_h


#include <deal.II/base/config.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/hp/dof_handler.h>

#include <vector>

DEAL_II_NAMESPACE_OPEN

template <typename number> class FullMatrix;
template <typename number> class SparseMatrix;
template <typename number> class Vector;
class ConstraintMatrix;

template <typename Accessor> class TriaRawIterator;

template <int, int> class FiniteElement;


namespace internal
{
  namespace DoFCellAccessor
  {
    struct Implementation;
  }

  namespace DoFHandler
  {
    struct Implementation;
    namespace Policy
    {
      struct Implementation;
    }
  }

  namespace hp
  {
    namespace DoFHandler
    {
      struct Implementation;
    }
  }
}

// note: the file dof_accessor.templates.h is included at the end of
// this file.  this includes a lot of templates and thus makes
// compilation slower, but at the same time allows for more aggressive
// inlining and thus faster code.


namespace internal
{
  namespace DoFAccessor
  {
    /**
     * This is a switch class which only declares a @p typedef. It is meant to
     * determine which class a DoFAccessor class is to be derived from. By
     * default, <tt>DoFAccessor@<structdim,dim,spacedim@></tt> derives from the
     * typedef in the general
     * <tt>Inheritance@<structdim,dim,spacedim@></tt> class, which is
     * <tt>TriaAccessor@<structdim,dim,spacedim@></tt>, but if
     * <tt>structdim==dim</tt>, then the specialization
     * <tt>Inheritance@<dim,dim,spacedim@></tt> is used which declares
     * its local type to be <tt>CellAccessor@<dim,spacedim@></tt>. Therefore, the
     * inheritance is automatically chosen to be from CellAccessor if the object
     * under consideration has full dimension, i.e. constitutes a cell.
     *
     * @ingroup dofs
     * @ingroup Accessors
     * @author Wolfgang Bangerth, 1999
     */
    template <int structdim, int dim, int spacedim>
    struct Inheritance
    {
      /**
       * Declaration of the @p typedef.
       * See the full documentation for
       * more information.
       */
      typedef dealii::TriaAccessor<structdim,dim,spacedim> BaseClass;
    };


    /**
     * This is the specialization of the general template used for the case where
     * an object has full dimension, i.e. is a cell. See the general template for
     * more details.
     */
    template <int dim, int spacedim>
    struct Inheritance<dim,dim,spacedim>
    {
      /**
       * Declaration of the @p typedef.
       * See the full documentation for
       * more information.
       */
      typedef dealii::CellAccessor<dim,spacedim> BaseClass;
    };
  }
}


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



/**
 * A class that gives access to the degrees of freedom stored in a DoFHandler,
 * MGSoDHandler, or hp::DoFHandler object. Accessors are used to, access the data that
 * pertains to edges, faces, and cells of a triangulation. The concept is
 * explained in more detail in connection to @ref Iterators.
 *
 * This class follows mainly the route laid out by the accessor
 * library declared in the triangulation library (TriaAccessor). It
 * enables the user to access the degrees of freedom on lines, quads,
 * or hexes. The first template argument of this class determines the
 * dimensionality of the object under consideration: 1 for lines, 2
 * for quads, and 3 for hexes. The second argument denotes the type of
 * DoF handler we should work on. It can either be ::DoFHandler or
 * hp::DoFHandler.  From the second template argument we also deduce
 * the dimension of the Triangulation this object refers to as well as
 * the dimension of the space into which it is embedded. Finally, the
 * template argument <code>level_dof_access</code> governs the
 * behavior of the function get_active_or_mg_dof_indices(). See the section on Generic loops below.
 *
 * <h3>Typedefs</h3>
 *
 * Usage is best to happen through the typedefs to the various kinds
 * of iterators provided by the DoFHandler and hp::DoFHandler classes,
 * since they are more secure to changes in the class naming and
 * template interface as well as providing easier typing (much less
 * complicated names!).
 *
 * <h3>Generic loops and the third template argument</h3>
 *
 * Many loops look very similar, whether they operate on the active
 * dofs of the active cells of the Triangulation or on the level dodfs
 * of a single level or the whole grid hierarchy. In order to use
 * polymorphism in such loops, they access degrees of freedom through
 * the function get_active_or_mg_dof_indices(), which changes behavior according to the
 * third template argument.  If the argument is false, then the active
 * dofs of active cells are accessed. If it is true, the level dofs
 * are used. DoFHandler has functions, for instance begin() and
 * begin_mg(), which return either type or the other. Additionally,
 * they can be cast into each other, in case this is needed, since
 * they access the same data.
 *
 * It is highly recommended to use the function get_active_or_mg_dof_indices() in
 * generic loops in lieu of get_dof_indices() or get_mg_dof_indices().
 *
 * <h3>Inheritance</h3>
 *
 * If the structural dimension given by the first template argument
 * equals the dimension of the DoFHandler (given as the second
 * template argument), then we are obviously dealing with cells,
 * rather than lower-dimensional objects. In that case, inheritance is
 * from CellAccessor, to provide access to all the cell specific
 * information afforded by that class. Otherwise, i.e. for
 * lower-dimensional objects, inheritance is from TriaAccessor.
 *
 * There is a DoFCellAccessor class that provides the
 * equivalent to the CellAccessor class.
 *
 * @ingroup dofs
 * @ingroup Accessors
 * @author Wolfgang Bangerth, 1998, 2006, 2008, Timo Heister, Guido Kanschat, 2012, 2013
 */
template <int structdim, class DH, bool level_dof_access>
class DoFAccessor : public dealii::internal::DoFAccessor::Inheritance<structdim, DH::dimension, DH::space_dimension>::BaseClass
{
public:

  /**
   * A static variable that allows users of
   * this class to discover the value of
   * the second template argument.
   */
  static const unsigned int dimension=DH::dimension;

  /**
   * A static variable that allows users of
   * this class to discover the value of
   * the third template argument.
   */
  static const unsigned int space_dimension=DH::space_dimension;

  /**
   * Declare a typedef to the base
   * class to make accessing some
   * of the exception classes
   * simpler.
   */
  typedef
  typename dealii::internal::DoFAccessor::Inheritance<structdim, dimension, space_dimension>::BaseClass
  BaseClass;

  /**
   * Data type passed by the iterator class.
   */
  typedef DH AccessorData;

  /**
   * @name Constructors
   */
  /**
   * @{
   */

  /**
   * Default constructor. Provides
   * an accessor that can't be
   * used.
   */
  DoFAccessor ();

  /**
   * Constructor
   */
  DoFAccessor (const Triangulation<DH::dimension,DH::space_dimension> *tria,
               const int                 level,
               const int                 index,
               const DH                 *local_data);

  /**
   * Conversion constructor. This
   * constructor exists to make certain
   * constructs simpler to write in
   * dimension independent code. For
   * example, it allows assigning a face
   * iterator to a line iterator, an
   * operation that is useful in 2d but
   * doesn't make any sense in 3d. The
   * constructor here exists for the
   * purpose of making the code conform to
   * C++ but it will unconditionally abort;
   * in other words, assigning a face
   * iterator to a line iterator is better
   * put into an if-statement that checks
   * that the dimension is two, and assign
   * to a quad iterator in 3d (an operator
   * that, without this constructor would
   * be illegal if we happen to compile for
   * 2d).
   */
  template <int structdim2, int dim2, int spacedim2>
  DoFAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  /**
   * Another conversion operator
   * between objects that don't
   * make sense, just like the
   * previous one.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  DoFAccessor (const DoFAccessor<dim2, DH2, level_dof_access2> &);

  /**
   * Copy constructor allowing to
   * switch level access and active
   * access.
   */
  template <bool level_dof_access2>
  DoFAccessor(const DoFAccessor<structdim, DH, level_dof_access2> &);
  /**
   * @}
   */

  /**
   * Return a handle on the
   * DoFHandler object which we
   * are using.
   */
  const DH &
  get_dof_handler () const;

  /**
   * Implement the copy operator needed
   * for the iterator classes.
   */
  template <bool level_dof_access2>
  void copy_from (const DoFAccessor<structdim, DH, level_dof_access2> &a);

  /**
   * Copy operator used by the
   * iterator class. Keeps the
   * previously set dof handler,
   * but sets the object
   * coordinates of the TriaAccessor.
   */
  void copy_from (const TriaAccessorBase<structdim, DH::dimension, DH::space_dimension> &da);

  /**
   * Tell the caller whether
   * get_active_or_mg_dof_indices() accesses active
   * or level dofs.
   */
  static bool is_level_cell();

  /**
   * Return an iterator pointing to
   * the the parent.
   */
  TriaIterator<DoFAccessor<structdim,DH, level_dof_access> >
  parent () const;

  /**
   *  @name Accessing sub-objects
   */
  /**
   * @{
   */

  /**
   * Return an iterator pointing to
   * the the @p c-th child.
   */
  TriaIterator<DoFAccessor<structdim,DH, level_dof_access> >
  child (const unsigned int c) const;

  /**
   * Pointer to the @p ith line
   * bounding this object. If the
  * current object is a line itself,
  * then the only valid index is
  * @p i equals to zero, and the
  * function returns an iterator
  * to itself.
   */
  typename dealii::internal::DoFHandler::Iterators<DH, level_dof_access>::line_iterator
  line (const unsigned int i) const;

  /**
   * Pointer to the @p ith quad
   * bounding this object. If the
  * current object is a quad itself,
  * then the only valid index is
  * @p i equals to zero, and the
  * function returns an iterator
  * to itself.
   */
  typename dealii::internal::DoFHandler::Iterators<DH, level_dof_access>::quad_iterator
  quad (const unsigned int i) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the DoF indices of this object
   */
  /**
   * @{
   */

  /**
   * Return the <i>global</i> indices of the degrees of freedom located on
   * this object in the standard ordering defined by the finite element (i.e.,
   * dofs on vertex 0, dofs on vertex 1, etc, dofs on line 0, dofs on line 1,
   * etc, dofs on quad 0, etc.) This function is only available on
   * <i>active</i> objects (see @ref GlossActive "this glossary entry").
   *
   * The cells needs to be an active cell (and not artificial in a
   * parallel distributed computation).
   *
   * The vector has to have the right size before being passed to this
   * function.
   *
   * The last argument denotes the finite element index. For the
   * standard ::DoFHandler class, this value must be equal to its
   * default value since that class only supports the same finite
   * element on all cells anyway.
   *
   * However, for hp objects (i.e. the hp::DoFHandler class),
   * different finite element objects may be used on different
   * cells. On faces between two cells, as well as vertices, there may
   * therefore be two sets of degrees of freedom, one for each of the
   * finite elements used on the adjacent cells. In order to specify
   * which set of degrees of freedom to work on, the last argument is
   * used to disambiguate. Finally, if this function is called for a
   * cell object, there can only be a single set of degrees of
   * freedom, and fe_index has to match the result of
   * active_fe_index().
   *
   * For cells, there is only a single possible finite element index
   * (namely the one for that cell, returned by
   * <code>cell-@>active_fe_index</code>. Consequently, the derived
   * DoFCellAccessor class has an overloaded version of this function
   * that calls the present function with
   * <code>cell-@>active_fe_index</code> as last argument.
   *
   */
  void get_dof_indices (std::vector<types::global_dof_index> &dof_indices,
                        const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Return the global multilevel indices of the degrees of freedom that live
   * on the current object with respect to the given level within the
   * multigrid hierarchy. The indices refer to the local numbering for the
   * level this line lives on.
   */
  void get_mg_dof_indices (const int level,
                           std::vector<types::global_dof_index> &dof_indices,
                           const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Sets the level DoF indices that are returned by get_mg_dof_indices.
   */
  void set_mg_dof_indices (const int level,
                           const std::vector<types::global_dof_index> &dof_indices,
                           const unsigned int fe_index = DH::default_fe_index);

  /**
   * Global DoF index of the <i>i</i>
   * degree associated with the @p vertexth
   * vertex of the present cell.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  types::global_dof_index vertex_dof_index (const unsigned int vertex,
                                            const unsigned int i,
                                            const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Returns the global DoF index of the <code>i</code>th degree of
   * freedom associated with the <code>vertex</code>th vertex on
   * level @p level. Also see vertex_dof_index().
   */
  types::global_dof_index mg_vertex_dof_index (const int level,
                                               const unsigned int vertex,
                                               const unsigned int i,
                                               const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Index of the <i>i</i>th degree
   * of freedom of this object.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   *
   * @note While the get_dof_indices()
   * function returns an array that
   * contains the indices of all degrees of
   * freedom that somehow live on this
   * object (i.e. on the vertices, edges or
   * interior of this object), the current
   * dof_index() function only considers
   * the DoFs that really belong to this
   * particular object's interior. In other
   * words, as an example, if the current
   * object refers to a quad (a cell in 2d,
   * a face in 3d) and the finite element
   * associated with it is a bilinear one,
   * then the get_dof_indices() will return
   * an array of size 4 while dof_index()
   * will produce an exception because no
   * degrees are defined in the interior of
   * the face.
   */
  types::global_dof_index dof_index (const unsigned int i,
                                     const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Returns the dof_index on the given level. Also see dof_index.
   */
  types::global_dof_index mg_dof_index (const int level, const unsigned int i) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the finite element associated with this object
   */
  /**
   * @{
   */

  /**
   * Return the number of finite
   * elements that are active on a
   * given object.
   *
   * For non-hp DoFHandler objects,
   * the answer is of course always
   * one. However, for
   * hp::DoFHandler objects, this
   * isn't the case: If this is a
   * cell, the answer is of course
   * one. If it is a face, the
   * answer may be one or two,
   * depending on whether the two
   * adjacent cells use the same
   * finite element or not. If it
   * is an edge in 3d, the possible
   * return value may be one or any
   * other value larger than that.
   */
  unsigned int
  n_active_fe_indices () const;

  /**
   * Return the @p n-th active fe
   * index on this object. For
   * cells and all non-hp objects,
   * there is only a single active
   * fe index, so the argument must
   * be equal to zero. For
   * lower-dimensional hp objects,
   * there are
   * n_active_fe_indices() active
   * finite elements, and this
   * function can be queried for
   * their indices.
   */
  unsigned int
  nth_active_fe_index (const unsigned int n) const;

  /**
   * Return true if the finite
   * element with given index is
   * active on the present
   * object. For non-hp DoF
   * accessors, this is of course
   * the case only if @p fe_index
   * equals zero. For cells, it is
   * the case if @p fe_index equals
   * active_fe_index() of this
   * cell. For faces and other
   * lower-dimensional objects,
   * there may be more than one @p
   * fe_index that are active on
   * any given object (see
   * n_active_fe_indices()).
   */
  bool
  fe_index_is_active (const unsigned int fe_index) const;

  /**
   * Return a reference to the finite
   * element used on this object with the
   * given @p fe_index. @p fe_index must be
   * used on this object,
   * i.e. <code>fe_index_is_active(fe_index)</code>
   * must return true.
   */
  const FiniteElement<DH::dimension,DH::space_dimension> &
  get_fe (const unsigned int fe_index) const;

  /**
   * @}
   */

  /**
   * Exceptions for child classes
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcInvalidObject);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcVectorNotEmpty);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcVectorDoesNotMatch);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcMatrixDoesNotMatch);
  /**
   * A function has been called for
   * a cell which should be active,
   * but is refined. @ref GlossActive
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcNotActive);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcCantCompareIterators);

protected:

  /**
   * Store the address of the DoFHandler object
   * to be accessed.
   */
  DH *dof_handler;
public:
  /**
   * Compare for equality. Return
   * <tt>true</tt> if the two
   * accessors refer to the same object.
   *
   * The template parameters of this
   * function allow for a comparison
   * of very different
   * objects. Therefore, some of them
   * are disabled.
   * Namely, if the dimension, or the dof
   * handler of the two objects
   * differ, an exception is
   * generated. It can be
   * expected that this is an
   * unwanted comparison.
   *
   * The template parameter
   * <tt>level_dof_access2</tt> is ignored, such
   * that an iterator with level
   * access can be equal to one with
   * access to the active degrees of
   * freedom.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  bool operator == (const DoFAccessor<dim2,DH2,level_dof_access2> &) const;

  /**
   * Compare for inequality. The
   * boolean not of operator==().
   */
  template <int dim2, class DH2, bool level_dof_access2>
  bool operator != (const DoFAccessor<dim2,DH2,level_dof_access2> &) const;
protected:
  /**
   * Reset the DoF handler pointer.
   */
  void set_dof_handler (DH *dh);

  /**
   * Set the index of the
   * <i>i</i>th degree of freedom
   * of this object to @p index.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  void set_dof_index (const unsigned int i,
                      const types::global_dof_index index,
                      const unsigned int fe_index = DH::default_fe_index) const;

  void set_mg_dof_index (const int level, const unsigned int i, const types::global_dof_index index) const;

  /**
   * Set the global index of the <i>i</i>
   * degree on the @p vertex-th vertex of
   * the present cell to @p index.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  void set_vertex_dof_index (const unsigned int vertex,
                             const unsigned int i,
                             const types::global_dof_index index,
                             const unsigned int fe_index = DH::default_fe_index) const;

  void set_mg_vertex_dof_index (const int level, const unsigned int vertex, const unsigned int i, const types::global_dof_index index, const unsigned int fe_index = DH::default_fe_index) const;

  /**
   * Iterator classes need to be friends
   * because they need to access operator==
   * and operator!=.
   */
  template <typename> friend class TriaRawIterator;
  template <int, class, bool> friend class DoFAccessor;

private:
  /**
   *  Copy operator. This is normally used
   *  in a context like <tt>iterator a,b;
   *  *a=*b;</tt>. Presumably, the intent
   *  here is to copy the object pointed to
   *  by @p b to the object pointed to by
   *  @p a. However, the result of
   *  dereferencing an iterator is not an
   *  object but an accessor; consequently,
   *  this operation is not useful for
   *  iterators on triangulations. We
   *  declare this function here private,
   *  thus it may not be used from outside.
   *  Furthermore it is not implemented and
   *  will give a linker error if used
   *  anyway.
   */
  DoFAccessor<structdim,DH, level_dof_access> &
  operator = (const DoFAccessor<structdim,DH, level_dof_access> &da);

  /**
   * Make the DoFHandler class a friend so
   * that it can call the set_xxx()
   * functions.
   */
  template <int dim, int spacedim> friend class DoFHandler;
  template <int dim, int spacedim> friend class hp::DoFHandler;

  friend struct dealii::internal::DoFHandler::Policy::Implementation;
  friend struct dealii::internal::DoFHandler::Implementation;
  friend struct dealii::internal::hp::DoFHandler::Implementation;
  friend struct dealii::internal::DoFCellAccessor::Implementation;
  friend struct dealii::internal::DoFAccessor::Implementation;
};



/**
 * Specialization of the general DoFAccessor class template for the
 * case of zero-dimensional objects (a vertex) that are the face of a
 * one-dimensional cell in spacedim space dimensions. Since vertices
 * function differently than general faces, this class does a few
 * things differently than the general template, but the interface
 * should look the same.
 *
 * @author Wolfgang Bangerth, 2010
 */
template <template <int, int> class DH, int spacedim, bool level_dof_access>
class DoFAccessor<0,DH<1,spacedim>, level_dof_access> : public TriaAccessor<0,1,spacedim>
{
public:

  /**
   * A static variable that allows users of
   * this class to discover the value of
   * the second template argument.
   */
  static const unsigned int dimension=1;

  /**
   * A static variable that allows users of
   * this class to discover the value of
   * the third template argument.
   */
  static const unsigned int space_dimension=spacedim;

  /**
   * Declare a typedef to the base
   * class to make accessing some
   * of the exception classes
   * simpler.
   */
  typedef TriaAccessor<0,1,spacedim> BaseClass;

  /**
   * Data type passed by the iterator class.
   */
  typedef DH<1,spacedim> AccessorData;

  /**
   * @name Constructors
   */
  /**
   * @{
   */

  /**
   * Default constructor. Provides
   * an accessor that can't be
   * used.
   */
  DoFAccessor ();

  /**
   * Constructor to be used if the
   * object here refers to a vertex
   * of a one-dimensional
   * triangulation, i.e. a face of
   * the triangulation.
   *
   * Since there is no mapping from
   * vertices to cells, an accessor
   * object for a point has no way
   * to figure out whether it is at
   * the boundary of the domain or
   * not. Consequently, the second
   * argument must be passed by the
   * object that generates this
   * accessor -- e.g. a 1d cell
   * that can figure out whether
   * its left or right vertex are
   * at the boundary.
   *
   * The third argument is the
   * global index of the vertex we
   * point to.
   *
   * The fourth argument is a
   * pointer to the DoFHandler
   * object.
   *
   * This iterator can only be
   * called for one-dimensional
   * triangulations.
   */
  DoFAccessor (const Triangulation<1,spacedim> *tria,
               const typename TriaAccessor<0,1,spacedim>::VertexKind vertex_kind,
               const unsigned int    vertex_index,
               const DH<1,spacedim> *dof_handler);

  /**
   * Constructor. This constructor
   * exists in order to maintain
   * interface compatibility with
   * the other accessor
   * classes. However, it doesn't
   * do anything useful here and so
   * may not actually be called.
   */
  DoFAccessor (const Triangulation<1,spacedim> *,
               const int = 0,
               const int = 0,
               const DH<1,spacedim> *dof_handler = 0);

  /**
   * Conversion constructor. This
   * constructor exists to make certain
   * constructs simpler to write in
   * dimension independent code. For
   * example, it allows assigning a face
   * iterator to a line iterator, an
   * operation that is useful in 2d but
   * doesn't make any sense in 3d. The
   * constructor here exists for the
   * purpose of making the code conform to
   * C++ but it will unconditionally abort;
   * in other words, assigning a face
   * iterator to a line iterator is better
   * put into an if-statement that checks
   * that the dimension is two, and assign
   * to a quad iterator in 3d (an operator
   * that, without this constructor would
   * be illegal if we happen to compile for
   * 2d).
   */
  template <int structdim2, int dim2, int spacedim2>
  DoFAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  /**
   * Another conversion operator
   * between objects that don't
   * make sense, just like the
   * previous one.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  DoFAccessor (const DoFAccessor<dim2, DH2, level_dof_access2> &);

  /**
   * @}
   */

  /**
   * Return a handle on the
   * DoFHandler object which we
   * are using.
   */
  const DH<1,spacedim> &
  get_dof_handler () const;

  /**
   * Copy operator.
   */
  DoFAccessor<0,DH<1,spacedim>, level_dof_access> &
  operator = (const DoFAccessor<0,DH<1,spacedim>, level_dof_access> &da);

  /**
   * Implement the copy operator needed
   * for the iterator classes.
   */
  template <bool level_dof_access2>
  void copy_from (const DoFAccessor<0, DH<1,spacedim>, level_dof_access2> &a);

  /**
   * Copy operator used by the
   * iterator class. Keeps the
   * previously set dof handler,
   * but sets the object
   * coordinates of the TriaAccessor.
   */
  void copy_from (const TriaAccessorBase<0, 1, spacedim> &da);

  /**
   * Return an iterator pointing to
   * the the parent.
   */
  TriaIterator<DoFAccessor<0,DH<1,spacedim>, level_dof_access> >
  parent () const;

  /**
   *  @name Accessing sub-objects
   */
  /**
   * @{
   */

  /**
   * Return an iterator pointing to
   * the the @p c-th child.
   */
  TriaIterator<DoFAccessor<0,DH<1,spacedim>, level_dof_access > >
  child (const unsigned int c) const;

  /**
   * Pointer to the @p ith line
   * bounding this object. If the
  * current object is a line itself,
  * then the only valid index is
  * @p i equals to zero, and the
  * function returns an iterator
  * to itself.
   */
  typename dealii::internal::DoFHandler::Iterators<DH<1,spacedim>, level_dof_access>::line_iterator
  line (const unsigned int i) const;

  /**
   * Pointer to the @p ith quad
   * bounding this object. If the
  * current object is a quad itself,
  * then the only valid index is
  * @p i equals to zero, and the
  * function returns an iterator
  * to itself.
   */
  typename dealii::internal::DoFHandler::Iterators<DH<1,spacedim>, level_dof_access>::quad_iterator
  quad (const unsigned int i) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the DoF indices of this object
   */
  /**
   * @{
   */

  /**
   * Return the <i>global</i> indices of the degrees of freedom located on
   * this object in the standard ordering defined by the finite element (i.e.,
   * dofs on vertex 0, dofs on vertex 1, etc, dofs on line 0, dofs on line 1,
   * etc, dofs on quad 0, etc.) This function is only available on
   * <i>active</i> objects (see @ref GlossActive "this glossary entry").
   *
   * The cells needs to be an active cell (and not artificial in a
   * parallel distributed computation).
   *
   * The vector has to have the
   * right size before being passed
   * to this function.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   *
   * For cells, there is only a
   * single possible finite element
   * index (namely the one for that
   * cell, returned by
   * <code>cell-@>active_fe_index</code>. Consequently,
   * the derived DoFCellAccessor
   * class has an overloaded
   * version of this function that
   * calls the present function
   * with
   * <code>cell-@>active_fe_index</code>
   * as last argument.
   */
  void get_dof_indices (std::vector<types::global_dof_index> &dof_indices,
                        const unsigned int fe_index = AccessorData::default_fe_index) const;

  /**
   * Global DoF index of the <i>i</i>
   * degree associated with the @p vertexth
   * vertex of the present cell.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  types::global_dof_index vertex_dof_index (const unsigned int vertex,
                                            const unsigned int i,
                                            const unsigned int fe_index = AccessorData::default_fe_index) const;

  /**
   * Index of the <i>i</i>th degree
   * of freedom of this object.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   *
   * @note While the get_dof_indices()
   * function returns an array that
   * contains the indices of all degrees of
   * freedom that somehow live on this
   * object (i.e. on the vertices, edges or
   * interior of this object), the current
   * dof_index() function only considers
   * the DoFs that really belong to this
   * particular object's interior. In other
   * words, as an example, if the current
   * object refers to a quad (a cell in 2d,
   * a face in 3d) and the finite element
   * associated with it is a bilinear one,
   * then the get_dof_indices() will return
   * an array of size 4 while dof_index()
   * will produce an exception because no
   * degrees are defined in the interior of
   * the face.
   */
  types::global_dof_index dof_index (const unsigned int i,
                                     const unsigned int fe_index = AccessorData::default_fe_index) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the finite element associated with this object
   */
  /**
   * @{
   */

  /**
   * Return the number of finite
   * elements that are active on a
   * given object.
   *
   * For non-hp DoFHandler objects,
   * the answer is of course always
   * one. However, for
   * hp::DoFHandler objects, this
   * isn't the case: If this is a
   * cell, the answer is of course
   * one. If it is a face, the
   * answer may be one or two,
   * depending on whether the two
   * adjacent cells use the same
   * finite element or not. If it
   * is an edge in 3d, the possible
   * return value may be one or any
   * other value larger than that.
   */
  unsigned int
  n_active_fe_indices () const;

  /**
   * Return the @p n-th active fe
   * index on this object. For
   * cells and all non-hp objects,
   * there is only a single active
   * fe index, so the argument must
   * be equal to zero. For
   * lower-dimensional hp objects,
   * there are
   * n_active_fe_indices() active
   * finite elements, and this
   * function can be queried for
   * their indices.
   */
  unsigned int
  nth_active_fe_index (const unsigned int n) const;

  /**
   * Return true if the finite
   * element with given index is
   * active on the present
   * object. For non-hp DoF
   * accessors, this is of course
   * the case only if @p fe_index
   * equals zero. For cells, it is
   * the case if @p fe_index equals
   * active_fe_index() of this
   * cell. For faces and other
   * lower-dimensional objects,
   * there may be more than one @p
   * fe_index that are active on
   * any given object (see
   * n_active_fe_indices()).
   */
  bool
  fe_index_is_active (const unsigned int fe_index) const;

  /**
   * Return a reference to the finite
   * element used on this object with the
   * given @p fe_index. @p fe_index must be
   * used on this object,
   * i.e. <code>fe_index_is_active(fe_index)</code>
   * must return true.
   */
  const FiniteElement<DH<1,spacedim>::dimension,DH<1,spacedim>::space_dimension> &
  get_fe (const unsigned int fe_index) const;

  /**
   * @}
   */

  /**
   * Exceptions for child classes
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcInvalidObject);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcVectorNotEmpty);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcVectorDoesNotMatch);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcMatrixDoesNotMatch);
  /**
   * A function has been called for
   * a cell which should be active,
   * but is refined. @ref GlossActive
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcNotActive);
  /**
   * Exception
   *
   * @ingroup Exceptions
   */
  DeclException0 (ExcCantCompareIterators);

protected:

  /**
   * Store the address of the DoFHandler object
   * to be accessed.
   */
  DH<1,spacedim> *dof_handler;

  /**
   *  Compare for equality.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  bool operator == (const DoFAccessor<dim2,DH2,level_dof_access2> &) const;

  /**
   * Compare for inequality.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  bool operator != (const DoFAccessor<dim2,DH2,level_dof_access2> &) const;

  /**
   * Reset the DoF handler pointer.
   */
  void set_dof_handler (DH<1,spacedim> *dh);

  /**
   * Set the index of the
   * <i>i</i>th degree of freedom
   * of this object to @p index.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  void set_dof_index (const unsigned int i,
                      const types::global_dof_index index,
                      const unsigned int fe_index = AccessorData::default_fe_index) const;

  /**
   * Set the global index of the <i>i</i>
   * degree on the @p vertex-th vertex of
   * the present cell to @p index.
   *
   * The last argument denotes the
   * finite element index. For the
   * standard ::DoFHandler class,
   * this value must be equal to
   * its default value since that
   * class only supports the same
   * finite element on all cells
   * anyway.
   *
   * However, for hp objects
   * (i.e. the hp::DoFHandler
   * class), different finite
   * element objects may be used on
   * different cells. On faces
   * between two cells, as well as
   * vertices, there may therefore
   * be two sets of degrees of
   * freedom, one for each of the
   * finite elements used on the
   * adjacent cells. In order to
   * specify which set of degrees
   * of freedom to work on, the
   * last argument is used to
   * disambiguate. Finally, if this
   * function is called for a cell
   * object, there can only be a
   * single set of degrees of
   * freedom, and fe_index has to
   * match the result of
   * active_fe_index().
   */
  void set_vertex_dof_index (const unsigned int vertex,
                             const unsigned int i,
                             const types::global_dof_index index,
                             const unsigned int fe_index = AccessorData::default_fe_index) const;

  /**
   * Iterator classes need to be friends
   * because they need to access operator==
   * and operator!=.
   */
  template <typename> friend class TriaRawIterator;


  /**
   * Make the DoFHandler class a friend so
   * that it can call the set_xxx()
   * functions.
   */
  template <int, int> friend class DoFHandler;
  template <int, int> friend class hp::DoFHandler;

  friend struct dealii::internal::DoFHandler::Policy::Implementation;
  friend struct dealii::internal::DoFHandler::Implementation;
  friend struct dealii::internal::hp::DoFHandler::Implementation;
  friend struct dealii::internal::DoFCellAccessor::Implementation;
};


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


/**
 * Grant access to the degrees of freedom on a cell.
 *
 * Note that since for the class we derive from, i.e. <tt>DoFAccessor<dim></tt>,
 * the two template parameters are equal, the base class is actually derived from
 * CellAccessor, which makes the functions of this class available to the
 * DoFCellAccessor class as well.
 *
 * @ingroup dofs
 * @ingroup Accessors
 * @author Wolfgang Bangerth, 1998, Timo Heister, Guido Kanschat, 2012
 */
template <class DH, bool level_dof_access>
class DoFCellAccessor :  public DoFAccessor<DH::dimension,DH, level_dof_access>
{
public:
  /**
   * Extract dimension from DH.
   */
  static const unsigned int dim = DH::dimension;

  /**
   * Extract space dimension from DH.
   */
  static const unsigned int spacedim = DH::space_dimension;


  /**
   * Data type passed by the iterator class.
   */
  typedef DH AccessorData;

  /**
   * Declare a typedef to the base
   * class to make accessing some
   * of the exception classes
   * simpler.
   */
  typedef DoFAccessor<DH::dimension,DH, level_dof_access> BaseClass;

  /**
   * Define the type of the
   * container this is part of.
   */
  typedef DH Container;

  /**
   * @name Constructors and initialization
   */
  /**
   * @{
   */

  /**
   * Constructor
   */
  DoFCellAccessor (const Triangulation<DH::dimension,DH::space_dimension> *tria,
                   const int                 level,
                   const int                 index,
                   const AccessorData       *local_data);

  /**
   * Conversion constructor. This
   * constructor exists to make certain
   * constructs simpler to write in
   * dimension independent code. For
   * example, it allows assigning a face
   * iterator to a line iterator, an
   * operation that is useful in 2d but
   * doesn't make any sense in 3d. The
   * constructor here exists for the
   * purpose of making the code conform to
   * C++ but it will unconditionally abort;
   * in other words, assigning a face
   * iterator to a line iterator is better
   * put into an if-statement that checks
   * that the dimension is two, and assign
   * to a quad iterator in 3d (an operator
   * that, without this constructor would
   * be illegal if we happen to compile for
   * 2d).
   */
  template <int structdim2, int dim2, int spacedim2>
  DoFCellAccessor (const InvalidAccessor<structdim2,dim2,spacedim2> &);

  /**
   * Another conversion operator
   * between objects that don't
   * make sense, just like the
   * previous one.
   */
  template <int dim2, class DH2, bool level_dof_access2>
  explicit
  DoFCellAccessor (const DoFAccessor<dim2, DH2, level_dof_access2> &);

  /**
   * @}
   */

  /**
   * Return the parent as a DoF
   * cell iterator. This
   * function is needed since the
   * parent function of the base
   * class returns a cell accessor
   * without access to the DoF
   * data.
   */
  TriaIterator<DoFCellAccessor<DH, level_dof_access> >
  parent () const;

  /**
   *  @name Accessing sub-objects and neighbors
   */
  /**
   * @{
   */

  /**
   * Return the @p ith neighbor as
   * a DoF cell iterator. This
   * function is needed since the
   * neighbor function of the base
   * class returns a cell accessor
   * without access to the DoF
   * data.
   */
  TriaIterator<DoFCellAccessor<DH, level_dof_access> >
  neighbor (const unsigned int) const;

  /**
   * Return the @p ith child as a
   * DoF cell iterator. This
   * function is needed since the
   * child function of the base
   * class returns a cell accessor
   * without access to the DoF
   * data.
   */
  TriaIterator<DoFCellAccessor<DH, level_dof_access> >
  child (const unsigned int) const;

  /**
   * Return an iterator to the @p ith face
   * of this cell.
   *
   * This function is not implemented in
   * 1D, and maps to DoFAccessor::line
   * in 2D.
   */
  TriaIterator<DoFAccessor<DH::dimension-1,DH, level_dof_access> >
  face (const unsigned int i) const;

  /**
   * Return the result of the
   * @p neighbor_child_on_subface
   * function of the base class,
   * but convert it so that one can
   * also access the DoF data (the
   * function in the base class
   * only returns an iterator with
   * access to the triangulation
   * data).
   */
  TriaIterator<DoFCellAccessor<DH, level_dof_access> >
  neighbor_child_on_subface (const unsigned int face_no,
                             const unsigned int subface_no) const;

  /**
   * @}
   */

  /**
   *  @name Extracting values from global vectors
   */
  /**
   * @{
   */

  /**
   * Return the values of the given vector
   * restricted to the dofs of this
   * cell in the standard ordering: dofs
   * on vertex 0, dofs on vertex 1, etc,
   * dofs on line 0, dofs on line 1, etc,
   * dofs on quad 0, etc.
   *
   * The vector has to have the
   * right size before being passed
   * to this function. This
   * function is only callable for
   * active cells.
   *
   * The input vector may be either
   * a <tt>Vector<float></tt>,
   * Vector<double>, or a
   * BlockVector<double>, or a
   * PETSc or Trilinos vector if
   * deal.II is compiled to support
   * these libraries. It is in the
   * responsibility of the caller
   * to assure that the types of
   * the numbers stored in input
   * and output vectors are
   * compatible and with similar
   * accuracy.
   */
  template <class InputVector, typename number>
  void get_dof_values (const InputVector &values,
                       Vector<number>    &local_values) const;

  /**
   * Return the values of the given vector
   * restricted to the dofs of this
   * cell in the standard ordering: dofs
   * on vertex 0, dofs on vertex 1, etc,
   * dofs on line 0, dofs on line 1, etc,
   * dofs on quad 0, etc.
   *
   * The vector has to have the
   * right size before being passed
   * to this function. This
   * function is only callable for
   * active cells.
   *
   * The input vector may be either
   * a <tt>Vector<float></tt>,
   * Vector<double>, or a
   * BlockVector<double>, or a
   * PETSc or Trilinos vector if
   * deal.II is compiled to support
   * these libraries. It is in the
   * responsibility of the caller
   * to assure that the types of
   * the numbers stored in input
   * and output vectors are
   * compatible and with similar
   * accuracy.
   */
  template <class InputVector, typename ForwardIterator>
  void get_dof_values (const InputVector &values,
                       ForwardIterator    local_values_begin,
                       ForwardIterator    local_values_end) const;

  /**
   * Return the values of the given vector
   * restricted to the dofs of this
   * cell in the standard ordering: dofs
   * on vertex 0, dofs on vertex 1, etc,
   * dofs on line 0, dofs on line 1, etc,
   * dofs on quad 0, etc.
   *
   * The vector has to have the
   * right size before being passed
   * to this function. This
   * function is only callable for
   * active cells.
   *
   * The input vector may be either a
   * <tt>Vector<float></tt>,
   * Vector<double>, or a
   * BlockVector<double>, or a PETSc or
   * Trilinos vector if deal.II is
   * compiled to support these
   * libraries. It is in the
   * responsibility of the caller to
   * assure that the types of the numbers
   * stored in input and output vectors
   * are compatible and with similar
   * accuracy. The ConstraintMatrix
   * passed as an argument to this
   * function makes sure that constraints
   * are correctly distributed when the
   * dof values are calculated.
   */
  template <class InputVector, typename ForwardIterator>
  void get_dof_values (const ConstraintMatrix &constraints,
                       const InputVector      &values,
                       ForwardIterator         local_values_begin,
                       ForwardIterator         local_values_end) const;

  /**
   * This function is the counterpart to
   * get_dof_values(): it takes a vector
   * of values for the degrees of freedom
   * of the cell pointed to by this iterator
   * and writes these values into the global
   * data vector @p values. This function
   * is only callable for active cells.
   *
   * Note that for continuous finite
   * elements, calling this function affects
   * the dof values on neighboring cells as
   * well. It may also violate continuity
   * requirements for hanging nodes, if
   * neighboring cells are less refined than
   * the present one. These requirements
   * are not taken care of and must be
   * enforced by the user afterwards.
   *
   * The vector has to have the
   * right size before being passed
   * to this function.
   *
   * The output vector may be either a
   * Vector<float>,
   * Vector<double>, or a
   * BlockVector<double>, or a
   * PETSc vector if deal.II is compiled to
   * support these libraries. It is in the
   * responsibility of the caller to assure
   * that the types of the numbers stored
   * in input and output vectors are
   * compatible and with similar accuracy.
   */
  template <class OutputVector, typename number>
  void set_dof_values (const Vector<number> &local_values,
                       OutputVector         &values) const;

  /**
   * Return the interpolation of
   * the given finite element
   * function to the present
   * cell. In the simplest case,
   * the cell is a terminal one,
   * i.e. has no children; then,
   * the returned value is the
   * vector of nodal values on that
   * cell. You could then as well
   * get the desired values through
   * the @p get_dof_values
   * function. In the other case,
   * when the cell has children, we
   * use the restriction matrices
   * provided by the finite element
   * class to compute the
   * interpolation from the
   * children to the present cell.
   *
   * It is assumed that both
   * vectors already have the right
   * size beforehand.
   *
   * Unlike the get_dof_values()
   * function, this function works
   * on cells rather than to lines,
   * quads, and hexes, since
   * interpolation is presently
   * only provided for cells by the
   * finite element classes.
   */
  template <class InputVector, typename number>
  void get_interpolated_dof_values (const InputVector &values,
                                    Vector<number>    &interpolated_values) const;

  /**
   * This, again, is the
   * counterpart to
   * get_interpolated_dof_values():
   * you specify the dof values on
   * a cell and these are
   * interpolated to the children
   * of the present cell and set on
   * the terminal cells.
   *
   * In principle, it works as
   * follows: if the cell pointed
   * to by this object is terminal,
   * then the dof values are set in
   * the global data vector by
   * calling the set_dof_values()
   * function; otherwise, the
   * values are prolonged to each
   * of the children and this
   * function is called for each of
   * them.
   *
   * Using the
   * get_interpolated_dof_values()
   * and this function, you can
   * compute the interpolation of a
   * finite element function to a
   * coarser grid by first getting
   * the interpolated solution on a
   * cell of the coarse grid and
   * afterwards redistributing it
   * using this function.
   *
   * Note that for continuous
   * finite elements, calling this
   * function affects the dof
   * values on neighboring cells as
   * well. It may also violate
   * continuity requirements for
   * hanging nodes, if neighboring
   * cells are less refined than
   * the present one, or if their
   * children are less refined than
   * the children of this
   * cell. These requirements are
   * not taken care of and must be
   * enforced by the user
   * afterward.
   *
   * It is assumed that both
   * vectors already have the right
   * size beforehand. This function
   * relies on the existence of a
   * natural interpolation property
   * of finite element spaces of a
   * cell to its children, denoted
   * by the prolongation matrices
   * of finite element classes. For
   * some elements, the spaces on
   * coarse and fine grids are not
   * nested, in which case the
   * interpolation to a child is
   * not the identity; refer to the
   * documentation of the
   * respective finite element
   * class for a description of
   * what the prolongation matrices
   * represent in this case.
   *
   * Unlike the set_dof_values()
   * function, this function is
   * associated to cells rather
   * than to lines, quads, and
   * hexes, since interpolation is
   * presently only provided for
   * cells by the finite element
   * objects.
   *
   * The output vector may be either a
   * Vector<float>,
   * Vector<double>, or a
   * BlockVector<double>, or a
   * PETSc vector if deal.II is compiled to
   * support these libraries. It is in the
   * responsibility of the caller to assure
   * that the types of the numbers stored
   * in input and output vectors are
   * compatible and with similar accuracy.
   */
  template <class OutputVector, typename number>
  void set_dof_values_by_interpolation (const Vector<number> &local_values,
                                        OutputVector         &values) const;

  /**
   * Distribute a local (cell
   * based) vector to a global one
   * by mapping the local numbering
   * of the degrees of freedom to
   * the global one and entering
   * the local values into the
   * global vector.
   *
   * The elements are
   * <em>added</em> up to the
   * elements in the global vector,
   * rather than just set, since
   * this is usually what one
   * wants.
   */
  template <typename number, typename OutputVector>
  void
  distribute_local_to_global (const Vector<number> &local_source,
                              OutputVector         &global_destination) const;

  /**
   * Distribute a local (cell based)
   * vector in iterator format to a
   * global one by mapping the local
   * numbering of the degrees of freedom
   * to the global one and entering the
   * local values into the global vector.
   *
   * The elements are <em>added</em> up
   * to the elements in the global
   * vector, rather than just set, since
   * this is usually what one wants.
   */
  template <typename ForwardIterator, typename OutputVector>
  void
  distribute_local_to_global (ForwardIterator   local_source_begin,
                              ForwardIterator   local_source_end,
                              OutputVector     &global_destination) const;

  /**
   * Distribute a local (cell based)
   * vector in iterator format to a
   * global one by mapping the local
   * numbering of the degrees of freedom
   * to the global one and entering the
   * local values into the global vector.
   *
   * The elements are <em>added</em> up
   * to the elements in the global
   * vector, rather than just set, since
   * this is usually what one
   * wants. Moreover, the
   * ConstraintMatrix passed to this
   * function makes sure that also
   * constraints are eliminated in this
   * process.
   */
  template <typename ForwardIterator, typename OutputVector>
  void
  distribute_local_to_global (const ConstraintMatrix &constraints,
                              ForwardIterator         local_source_begin,
                              ForwardIterator         local_source_end,
                              OutputVector           &global_destination) const;

  /**
   * This function does much the
   * same as the
   * <tt>distribute_local_to_global(Vector,Vector)</tt>
   * function, but operates on
   * matrices instead of
   * vectors. If the matrix type is
   * a sparse matrix then it is
   * supposed to have non-zero
   * entry slots where required.
   */
  template <typename number, typename OutputMatrix>
  void
  distribute_local_to_global (const FullMatrix<number> &local_source,
                              OutputMatrix             &global_destination) const;

  /**
   * This function does what the two
   * <tt>distribute_local_to_global</tt>
   * functions with vector and matrix
   * argument do, but all at once.
   */
  template <typename number, typename OutputMatrix, typename OutputVector>
  void
  distribute_local_to_global (const FullMatrix<number> &local_matrix,
                              const Vector<number>     &local_vector,
                              OutputMatrix             &global_matrix,
                              OutputVector             &global_vector) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the DoF indices of this object
   */
  /**
   * @{
   */

  /**
   * Obtain the global indices of the local degrees of freedom on this cell.
   *
   * If this object accesses a level
   * cell (indicated by the third
   * template argument or
   * #is_level_cell), then return the
   * result of get_mg_dof_indices(),
   * else return get_dof_indices().
   *
   * You will get a
   * level_cell_iterator when calling
   * begin_mg() and a normal one
   * otherwise.
   *
   * Examples for this use are in the implementation of DoFRenumbering.
   */
  void get_active_or_mg_dof_indices (std::vector<types::global_dof_index> &dof_indices) const;

  /**
   * Return the <i>global</i> indices of the degrees of freedom located on
   * this object in the standard ordering defined by the finite element (i.e.,
   * dofs on vertex 0, dofs on vertex 1, etc, dofs on line 0, dofs on line 1,
   * etc, dofs on quad 0, etc.) This function is only available on
   * <i>active</i> objects (see @ref GlossActive "this glossary entry").
   *
   * @param[out] dof_indices The vector into which the indices will be
   *     written. It has to have the right size (namely,
   *     <code>fe.dofs_per_cell</code>, <code>fe.dofs_per_face</code>,
   *     or <code>fe.dofs_per_line</code>, depending on which kind of
   *     object this function is called) before being passed to this
   *     function.
   *
   * This function reimplements the same function in the base class.
   * In contrast to the function in the base class, we do not need the
   * <code>fe_index</code> here because there is always a unique finite
   * element index on cells.
   *
   * This is a function which requires that the cell is active.
   *
   * Also see get_active_or_mg_dof_indices().
   *
   * @note In many places in the tutorial and elsewhere in the library,
   *   the argument to this function is called <code>local_dof_indices</code>
   *   by convention. The name is not meant to indicate the <i>local</i>
   *   numbers of degrees of freedom (which are always between zero and
   *   <code>fe.dofs_per_cell</code>) but instead that the returned values
   *   are the <i>global</i> indices of those degrees of freedom that
   *   are located locally on the current cell.
   *
   * @deprecated Currently, this function can also be called for non-active cells, if all degrees of freedom of the FiniteElement are located in vertices. This functionality will vanish in a future release.
   */
  void get_dof_indices (std::vector<types::global_dof_index> &dof_indices) const;

  /**
   * @deprecated Use get_active_or_mg_dof_indices() with level_cell_iterator returned from begin_mg().
   *
   * Retrieve the global indices of the degrees of freedom on this cell in the
   * level vector associated to the level of the cell.
   */
  void get_mg_dof_indices (std::vector<types::global_dof_index> &dof_indices) const;

  /**
   * @}
   */

  /**
   *  @name Accessing the finite element associated with this object
   */
  /**
   * @{
   */

  /**
   * Return the finite element that
   * is used on the cell pointed to
   * by this iterator. For non-hp
   * DoF handlers, this is of
   * course always the same
   * element, independent of the
   * cell we are presently on, but
   * for hp DoF handlers, this may
   * change from cell to cell.
   */
  const FiniteElement<DH::dimension,DH::space_dimension> &
  get_fe () const;

  /**
   *  Returns the index inside the
   *  hp::FECollection of the FiniteElement
   *  used for this cell.
   */
  unsigned int active_fe_index () const;

  /**
   *  Sets the index of the FiniteElement used for
   *  this cell.
   */
  void set_active_fe_index (const unsigned int i);
  /**
   * @}
   */

  /**
   * Set the DoF indices of this
   * cell to the given values. This
   * function bypasses the DoF
   * cache, if one exists for the
   * given DoF handler class.
   */
  void set_dof_indices (const std::vector<types::global_dof_index> &dof_indices);

  /**
   * Set the Level DoF indices of this
   * cell to the given values.
   */
  void set_mg_dof_indices (const std::vector<types::global_dof_index> &dof_indices);

  /**
   * Update the cache in which we
   * store the dof indices of this
   * cell.
   */
  void update_cell_dof_indices_cache () const;

private:
  /**
   *  Copy operator. This is normally used
   *  in a context like <tt>iterator a,b;
   *  *a=*b;</tt>. Presumably, the intent
   *  here is to copy the object pointed to
   *  by @p b to the object pointed to by
   *  @p a. However, the result of
   *  dereferencing an iterator is not an
   *  object but an accessor; consequently,
   *  this operation is not useful for
   *  iterators on triangulations. We
   *  declare this function here private,
   *  thus it may not be used from outside.
   *  Furthermore it is not implemented and
   *  will give a linker error if used
   *  anyway.
   */
  DoFCellAccessor<DH, level_dof_access> &
  operator = (const DoFCellAccessor<DH, level_dof_access> &da);

  /**
   * Make the DoFHandler class a
   * friend so that it can call the
   * update_cell_dof_indices_cache()
   * function
   */
  template <int dim, int spacedim> friend class DoFHandler;
  friend struct dealii::internal::DoFCellAccessor::Implementation;
};


template <int sd, class DH, bool lda>
inline
bool
DoFAccessor<sd, DH, lda>::is_level_cell()
{
  return lda;
}



DEAL_II_NAMESPACE_CLOSE

// include more templates
#include "dof_accessor.templates.h"


#endif