]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Next round of multigrid implementations.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 26 Aug 1998 08:27:02 +0000 (08:27 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 26 Aug 1998 08:27:02 +0000 (08:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@522 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/mg_dof_accessor.h
deal.II/deal.II/include/dofs/mg_dof_handler.h
deal.II/deal.II/include/multigrid/mg_dof_accessor.h
deal.II/deal.II/include/multigrid/mg_dof_handler.h
deal.II/deal.II/source/dofs/mg_dof_accessor.cc
deal.II/deal.II/source/dofs/mg_dof_handler.cc
deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
deal.II/deal.II/source/multigrid/mg_dof_handler.cc

index 03274755dc4537b30b651fe83c76c715d3e7442a..75ba6a5a746c2e31a992e15bf5b76f8e3e2d1429 100644 (file)
@@ -51,15 +51,22 @@ class MGDoFAccessor {
 
                                     /**
                                      * This should be the default constructor.
+                                     * We cast away the #const#ness of the
+                                     * pointer which clearly is EVIL but
+                                     * we can't help without making all
+                                     * functions which could somehow use
+                                     * iterators (directly or indirectly) make
+                                     * non-const, even if they preserve
+                                     * constness.
                                      */
-    MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
-                   mg_dof_handler(mg_dof_handler) {};
+    MGDoFAccessor (const MGDoFHandler<dim> *mg_dof_handler) :
+                   mg_dof_handler(const_cast<MGDoFHandler<dim>*>(mg_dof_handler)) {};
 
                                     /**
                                      * Reset the DoF handler pointer.
                                      */
     void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
-      Assert (dh != 0, ExcInvalidObject());
+      Assert (dh != 0, DoFAccessor<dim>::ExcInvalidObject());
       mg_dof_handler = dh;
     };
 
@@ -133,9 +140,15 @@ class MGDoFAccessor {
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim, typename BaseClass>
-class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
-                          public DoFLineAccessor<dim, BaseClass> {
+class MGDoFLineAccessor :  public MGDoFAccessor<dim>, public DoFLineAccessor<dim, BaseClass> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFHandler<dim> AccessorData;
+    
                                     /**
                                      * Default constructor, unused thus
                                      * not implemented.
@@ -150,7 +163,7 @@ class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
     MGDoFLineAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data);
+                      const AccessorData *local_data);
     
                                     /**
                                      * Return the index of the #i#th degree
@@ -223,9 +236,15 @@ class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
  * @see DoFLineAccessor
  */
 template <int dim, typename BaseClass>
-class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
-                          public DoFQuadAccessor<dim, BaseClass> {
+class MGDoFQuadAccessor :  public MGDoFAccessor<dim>, public DoFQuadAccessor<dim, BaseClass> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFHandler<dim> AccessorData;
+
                                     /**
                                      * Default constructor, unused thus
                                      * not implemented.
@@ -240,7 +259,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
     MGDoFQuadAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data);
+                      const AccessorData *local_data);
     
                                     /**
                                      * Return the index of the #i#th degree
@@ -307,7 +326,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
                                      * Implement the copy operator needed
                                      * for the iterator classes.
                                      */
-    void copy_from (const DoFQuadAccessor<dim,BaseClass> &a);
+    void copy_from (const MGDoFQuadAccessor<dim,BaseClass> &a);
 };
 
 
@@ -329,8 +348,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
  * the compiler) happens to use this class.
  */
 template <int dim>
-class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
-                              public TriaAccessor<dim> {
+class MGDoFSubstructAccessor : public MGDoFLineAccessor<1,CellAccessor<1> > {
   public:
     DoFSubstructAccessor () {
       Assert (false, ExcInternalError());
@@ -351,8 +369,8 @@ class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
  * with template parameters. This class and #DoFSubstructAccessor<2>#
  * wrap the following names:
  * \begin{verbatim}
- *   DoFSubstructAccessor<1> := DoFLineAccessor<1,CellAccessor<1> >;
- *   DoFSubstructAccessor<2> := DoFQuadAccessor<2,CellAccessor<2> >;
+ *   MGDoFSubstructAccessor<1> := MGDoFLineAccessor<1,CellAccessor<1> >;
+ *   MGDoFSubstructAccessor<2> := MGDoFQuadAccessor<2,CellAccessor<2> >;
  * \end{verbatim}
  * We do this rather complex (and needless, provided C++ the needed constructs!)
  * class hierarchy manipulation, since this way we can declare and implement
@@ -366,13 +384,20 @@ class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
 template <>
 class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> > {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef typename MGDoFLineAccessor<1,CellAccessor<1> >::AccessorData AccessorData;
+
                                     /**
                                      * Constructor
                                      */
     MGDoFSubstructAccessor (Triangulation<1> *tria,
                            const int         level,
                            const int         index,
-                           const void       *local_data) :
+                           const AccessorData *local_data) :
                    MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
@@ -393,13 +418,20 @@ class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> >
 template <>
 class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFQuadAccessor<2,CellAccessor<2> >::AccessorData AccessorData;
+    
                                     /**
                                      * Constructor
                                      */
     MGDoFSubstructAccessor (Triangulation<2> *tria,
                            const int         level,
                            const int         index,
-                           const void       *local_data) :
+                           const AccessorData *local_data) :
                    MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
@@ -431,13 +463,20 @@ class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
 template <int dim>
 class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef typename MGDoFSubstructAccessor<dim>::AccessorData AccessorData;
+    
                                     /**
                                      * Constructor
                                      */
     MGDoFCellAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data) :
+                      const AccessorData *local_data) :
                    MGDoFSubstructAccessor<dim> (tria,level,index,local_data) {};
 
                                     /**
index 21ad59db2066bef73d069211d0fe36965501e5e7..0e52401f3750e54df73d3c6c39b902a8c4c35693 100644 (file)
 
 // forward declarations
 template <int dim, typename BaseClass> class MGDoFLineAccessor;
-template <int dim, typename BaseClass> class MGDoFLineAccessor;
-template <int dim, typename BaseClass> class MGDoFQuadAccessor;
 template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+template <int dim> class MGDoFCellAccessor;
+
+
+
+
+/**
+ * Define some types which differ between the dimensions. This class
+ * is analogous to the \Ref{TriaDimensionInfo} class hierarchy.
+ * 
+ * @see MGDoFDimensionInfo<1>
+ * @see MGDoFDimensionInfo<2>
+ */
+template <int dim>
+class MGDoFDimensionInfo;
+
+
+
+
+
+/**
+ * Define some types for the DoF handling in one dimension.
+ *
+ * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
+ */
+class MGDoFDimensionInfo<1> {
+  public:
+    typedef TriaRawIterator<1,MGDoFCellAccessor<1> >    raw_line_iterator;
+    typedef TriaIterator<1,MGDoFCellAccessor<1> >       line_iterator;
+    typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator;
+
+    typedef void * raw_quad_iterator;
+    typedef void * quad_iterator;
+    typedef void * active_quad_iterator;
+
+    typedef raw_line_iterator    raw_cell_iterator;
+    typedef line_iterator        cell_iterator;
+    typedef active_line_iterator active_cell_iterator;
+
+    typedef void * raw_face_iterator;
+    typedef void * face_iterator;
+    typedef void * active_face_iterator;
+};
+
+
+
+
+
+
+/**
+ * Define some types for the DoF handling in two dimensions.
+ *
+ * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
+ */
+class MGDoFDimensionInfo<2> {
+  public:
+    typedef TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >    raw_line_iterator;
+    typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >       line_iterator;
+    typedef TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > active_line_iterator;
+    
+    typedef TriaRawIterator<2,MGDoFCellAccessor<2> >               raw_quad_iterator;
+    typedef TriaIterator<2,MGDoFCellAccessor<2> >                  quad_iterator;
+    typedef TriaActiveIterator<2,MGDoFCellAccessor<2> >            active_quad_iterator;
+
+    typedef raw_quad_iterator    raw_cell_iterator;
+    typedef quad_iterator        cell_iterator;
+    typedef active_quad_iterator active_cell_iterator;
+
+    typedef raw_line_iterator    raw_face_iterator;
+    typedef line_iterator        face_iterator;
+    typedef active_line_iterator active_face_iterator;    
+};
 
 
 
@@ -22,7 +91,26 @@ template <int dim, typename BaseClass> class MGDoFQuadAccessor;
 template <int dim>
 class MGDoFHandler : public DoFHandler<dim> {
   public:
-                                    /**
+                                    // insert these definitions for gcc2.8,
+                                    // since it can't inherit typedefs (I
+                                    // believe it should, but it can't)
+    typedef typename MGDoFDimensionInfo<dim>::raw_line_iterator raw_line_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::line_iterator line_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_line_iterator active_line_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_quad_iterator raw_quad_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::face_iterator face_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
+
+                                    /**
                                      * Constructor. Take #tria# as the
                                      * triangulation to work on.
                                      */
@@ -69,6 +157,473 @@ class MGDoFHandler : public DoFHandler<dim> {
                                const vector<int> &starting_points = vector<int>());
 
 
+                                    /*--------------------------------------*/
+    
+                                    /**
+                                     *  @name Cell iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first cell, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no cells, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 1D and #begin_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    begin_raw   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used cell
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 1D and #begin_quad# in 2D.
+                                     */
+    cell_iterator        begin       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  cell on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 1D and #begin_active_quad# in 2D.
+                                     */
+    active_cell_iterator begin_active(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 1D and #end_quad# in 2D.
+                                     */
+    raw_cell_iterator    end () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    cell_iterator        end (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_cell_iterator    end_raw (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_cell_iterator end_active (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last cell, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 1D and #last_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    last_raw () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  cell of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 1D and #last_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    last_raw (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used cell.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 1D and #last_quad# in 2D.
+                                     */
+    cell_iterator        last () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used cell on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 1D and #last_quad# in 2D.
+                                     */
+    cell_iterator        last (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active cell.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 1D and #last_active_quad# in 2D.
+                                     */
+    active_cell_iterator last_active () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active cell on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 1D and #last_active_quad# in 2D.
+                                     */
+    active_cell_iterator last_active (const unsigned int level) const;
+                                    //@}
+
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Face iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first face, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no faces, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 2D and #begin_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    begin_raw_face   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used face
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 2D and #begin_quad# in 3D.
+                                     */
+    face_iterator        begin_face       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  face on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 2D and #begin_active_quad# in 3D.
+                                     */
+    active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 2D and #end_quad# in 3D.
+                                     */
+    raw_face_iterator    end_face () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    face_iterator        end_face (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_face_iterator    end_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_face_iterator end_active_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last face, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  face of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face (const unsigned int level) const;
+                                    //@}
+
+    
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Line iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first line, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no lines, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_line_iterator
+    begin_raw_line   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used line
+                                     *  on level #level#.
+                                     */
+    line_iterator
+    begin_line       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  line on level #level#.
+                                     */
+    active_line_iterator
+    begin_active_line(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_line_iterator
+    end_line () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    line_iterator        end_line (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_line_iterator    end_raw_line (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_line_iterator end_active_line (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last line, used or not.
+                                     */
+    raw_line_iterator
+    last_raw_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  line of the level #level#, used or not.
+
+                                     */
+    raw_line_iterator
+    last_raw_line (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used line.
+                                     */
+    line_iterator
+    last_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used line on level #level#.
+                                     */
+    line_iterator
+    last_line (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active line.
+                                     */
+    active_line_iterator
+    last_active_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active line on level #level#.
+                                     */
+    active_line_iterator
+    last_active_line (const unsigned int level) const;
+                                    /*@}*/       
+
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Quad iterator functions*/
+                                    /*@{
+                                     */
+                                    /**
+                                     *  Return iterator to the first quad, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no quads, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_quad_iterator
+    begin_raw_quad   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used quad
+                                     *  on level #level#.
+                                     */
+    quad_iterator
+    begin_quad       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  quad on level #level#.
+                                     */
+    active_quad_iterator
+    begin_active_quad(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_quad_iterator
+    end_quad () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    quad_iterator        end_quad (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_quad_iterator    end_raw_quad (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_quad_iterator end_active_quad (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last quad, used or not.
+                                     */
+    raw_quad_iterator
+    last_raw_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  quad of the level #level#, used or not.
+
+                                     */
+    raw_quad_iterator
+    last_raw_quad (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used quad.
+                                     */
+    quad_iterator
+    last_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used quad on level #level#.
+                                     */
+    quad_iterator
+    last_quad (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active quad.
+                                     */
+    active_quad_iterator
+    last_active_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active quad on level #level#.
+                                     */
+    active_quad_iterator
+    last_active_quad (const unsigned int level) const;
+                                    /*@}*/
+
+                                    /*---------------------------------------*/
+    
   private:
 
                                     /** We need for each vertex a list of the
@@ -93,12 +648,23 @@ class MGDoFHandler : public DoFHandler<dim> {
       public:
        
                                         /**
-                                         * Constructor. Allocates memory and
-                                         * sets all indices to #-1#.
+                                         * Constructor. This one is empty
+                                         * because it is difficult to make it
+                                         * efficient to use vector<>'s and
+                                         * still construct the object using
+                                         * the constructor. Use the #init#
+                                         * function to really allocate
+                                         * memory.
                                          */
-       MGVertexDoFs (const unsigned int coarsest_level,
-                     const unsigned int n_levels,
-                     const unsigned int dofs_per_vertex);
+       MGVertexDoFs ();
+
+                                        /**
+                                         * Allocate memory and
+                                         * set all indices to #-1#.
+                                         */
+       void init (const unsigned int coarsest_level,
+                  const unsigned int n_levels,
+                  const unsigned int dofs_per_vertex);
 
                                         /**
                                          * Destructor
@@ -148,6 +714,7 @@ class MGDoFHandler : public DoFHandler<dim> {
        DeclException0 (ExcNoMemory);
                                         /**
                                          * Exception.
+
                                          */
        DeclException0 (ExcInvalidIndex);
                                         /**
@@ -157,6 +724,10 @@ class MGDoFHandler : public DoFHandler<dim> {
                        int,
                        << "The given level number " << arg1 << " is below the "
                        << "coarsest level this vertex lives on.");
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcInternalError);
        
       private:
                                         /**
@@ -165,12 +736,12 @@ class MGDoFHandler : public DoFHandler<dim> {
                                          * as an offset when accessing the
                                          * dofs of a specific level.
                                          */
-       const unsigned int coarsest_level;
+       unsigned int coarsest_level;
 
                                         /**
                                          * Array holding the indices.
                                          */
-       int *const         indices;
+       int         *indices;
     };
 
     
@@ -180,6 +751,23 @@ class MGDoFHandler : public DoFHandler<dim> {
                                      */
     void reserve_space ();
 
+                                    /**
+                                     * Distribute dofs on the given cell,
+                                     * with new dofs starting with index
+                                     * #next_free_dof#. Return the next
+                                     * unused index number. The finite
+                                     * element used is the one given to
+                                     * #distribute_dofs#, which is copied
+                                     * to #selected_fe#.
+                                     *
+                                     * This function is excluded from the
+                                     * #distribute_dofs# function since
+                                     * it can not be implemented dimension
+                                     * independent.
+                                     */
+    unsigned int distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int next_free_dof);
+
                                     /**
                                      * Space to store the DoF numbers for the
                                      * different levels. Unlike the #levels#
index 03274755dc4537b30b651fe83c76c715d3e7442a..75ba6a5a746c2e31a992e15bf5b76f8e3e2d1429 100644 (file)
@@ -51,15 +51,22 @@ class MGDoFAccessor {
 
                                     /**
                                      * This should be the default constructor.
+                                     * We cast away the #const#ness of the
+                                     * pointer which clearly is EVIL but
+                                     * we can't help without making all
+                                     * functions which could somehow use
+                                     * iterators (directly or indirectly) make
+                                     * non-const, even if they preserve
+                                     * constness.
                                      */
-    MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
-                   mg_dof_handler(mg_dof_handler) {};
+    MGDoFAccessor (const MGDoFHandler<dim> *mg_dof_handler) :
+                   mg_dof_handler(const_cast<MGDoFHandler<dim>*>(mg_dof_handler)) {};
 
                                     /**
                                      * Reset the DoF handler pointer.
                                      */
     void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
-      Assert (dh != 0, ExcInvalidObject());
+      Assert (dh != 0, DoFAccessor<dim>::ExcInvalidObject());
       mg_dof_handler = dh;
     };
 
@@ -133,9 +140,15 @@ class MGDoFAccessor {
  * @author Wolfgang Bangerth, 1998
  */
 template <int dim, typename BaseClass>
-class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
-                          public DoFLineAccessor<dim, BaseClass> {
+class MGDoFLineAccessor :  public MGDoFAccessor<dim>, public DoFLineAccessor<dim, BaseClass> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFHandler<dim> AccessorData;
+    
                                     /**
                                      * Default constructor, unused thus
                                      * not implemented.
@@ -150,7 +163,7 @@ class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
     MGDoFLineAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data);
+                      const AccessorData *local_data);
     
                                     /**
                                      * Return the index of the #i#th degree
@@ -223,9 +236,15 @@ class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
  * @see DoFLineAccessor
  */
 template <int dim, typename BaseClass>
-class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
-                          public DoFQuadAccessor<dim, BaseClass> {
+class MGDoFQuadAccessor :  public MGDoFAccessor<dim>, public DoFQuadAccessor<dim, BaseClass> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFHandler<dim> AccessorData;
+
                                     /**
                                      * Default constructor, unused thus
                                      * not implemented.
@@ -240,7 +259,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
     MGDoFQuadAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data);
+                      const AccessorData *local_data);
     
                                     /**
                                      * Return the index of the #i#th degree
@@ -307,7 +326,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
                                      * Implement the copy operator needed
                                      * for the iterator classes.
                                      */
-    void copy_from (const DoFQuadAccessor<dim,BaseClass> &a);
+    void copy_from (const MGDoFQuadAccessor<dim,BaseClass> &a);
 };
 
 
@@ -329,8 +348,7 @@ class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
  * the compiler) happens to use this class.
  */
 template <int dim>
-class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
-                              public TriaAccessor<dim> {
+class MGDoFSubstructAccessor : public MGDoFLineAccessor<1,CellAccessor<1> > {
   public:
     DoFSubstructAccessor () {
       Assert (false, ExcInternalError());
@@ -351,8 +369,8 @@ class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
  * with template parameters. This class and #DoFSubstructAccessor<2>#
  * wrap the following names:
  * \begin{verbatim}
- *   DoFSubstructAccessor<1> := DoFLineAccessor<1,CellAccessor<1> >;
- *   DoFSubstructAccessor<2> := DoFQuadAccessor<2,CellAccessor<2> >;
+ *   MGDoFSubstructAccessor<1> := MGDoFLineAccessor<1,CellAccessor<1> >;
+ *   MGDoFSubstructAccessor<2> := MGDoFQuadAccessor<2,CellAccessor<2> >;
  * \end{verbatim}
  * We do this rather complex (and needless, provided C++ the needed constructs!)
  * class hierarchy manipulation, since this way we can declare and implement
@@ -366,13 +384,20 @@ class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
 template <>
 class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> > {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef typename MGDoFLineAccessor<1,CellAccessor<1> >::AccessorData AccessorData;
+
                                     /**
                                      * Constructor
                                      */
     MGDoFSubstructAccessor (Triangulation<1> *tria,
                            const int         level,
                            const int         index,
-                           const void       *local_data) :
+                           const AccessorData *local_data) :
                    MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
@@ -393,13 +418,20 @@ class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> >
 template <>
 class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef MGDoFQuadAccessor<2,CellAccessor<2> >::AccessorData AccessorData;
+    
                                     /**
                                      * Constructor
                                      */
     MGDoFSubstructAccessor (Triangulation<2> *tria,
                            const int         level,
                            const int         index,
-                           const void       *local_data) :
+                           const AccessorData *local_data) :
                    MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
                                     // do this here, since this way the
                                     // CellAccessor has the possibility to know
@@ -431,13 +463,20 @@ class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
 template <int dim>
 class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
   public:
+                                    /**
+                                     * Declare the data type that this accessor
+                                     * class expects to get passed from the
+                                     * iterator classes.
+                                     */
+    typedef typename MGDoFSubstructAccessor<dim>::AccessorData AccessorData;
+    
                                     /**
                                      * Constructor
                                      */
     MGDoFCellAccessor (Triangulation<dim> *tria,
                       const int           level,
                       const int           index,
-                      const void         *local_data) :
+                      const AccessorData *local_data) :
                    MGDoFSubstructAccessor<dim> (tria,level,index,local_data) {};
 
                                     /**
index 21ad59db2066bef73d069211d0fe36965501e5e7..0e52401f3750e54df73d3c6c39b902a8c4c35693 100644 (file)
 
 // forward declarations
 template <int dim, typename BaseClass> class MGDoFLineAccessor;
-template <int dim, typename BaseClass> class MGDoFLineAccessor;
-template <int dim, typename BaseClass> class MGDoFQuadAccessor;
 template <int dim, typename BaseClass> class MGDoFQuadAccessor;
+template <int dim> class MGDoFCellAccessor;
+
+
+
+
+/**
+ * Define some types which differ between the dimensions. This class
+ * is analogous to the \Ref{TriaDimensionInfo} class hierarchy.
+ * 
+ * @see MGDoFDimensionInfo<1>
+ * @see MGDoFDimensionInfo<2>
+ */
+template <int dim>
+class MGDoFDimensionInfo;
+
+
+
+
+
+/**
+ * Define some types for the DoF handling in one dimension.
+ *
+ * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
+ */
+class MGDoFDimensionInfo<1> {
+  public:
+    typedef TriaRawIterator<1,MGDoFCellAccessor<1> >    raw_line_iterator;
+    typedef TriaIterator<1,MGDoFCellAccessor<1> >       line_iterator;
+    typedef TriaActiveIterator<1,MGDoFCellAccessor<1> > active_line_iterator;
+
+    typedef void * raw_quad_iterator;
+    typedef void * quad_iterator;
+    typedef void * active_quad_iterator;
+
+    typedef raw_line_iterator    raw_cell_iterator;
+    typedef line_iterator        cell_iterator;
+    typedef active_line_iterator active_cell_iterator;
+
+    typedef void * raw_face_iterator;
+    typedef void * face_iterator;
+    typedef void * active_face_iterator;
+};
+
+
+
+
+
+
+/**
+ * Define some types for the DoF handling in two dimensions.
+ *
+ * The types have the same meaning as those declared in \Ref{TriaDimensionInfo<2>}.
+ */
+class MGDoFDimensionInfo<2> {
+  public:
+    typedef TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >    raw_line_iterator;
+    typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >       line_iterator;
+    typedef TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > active_line_iterator;
+    
+    typedef TriaRawIterator<2,MGDoFCellAccessor<2> >               raw_quad_iterator;
+    typedef TriaIterator<2,MGDoFCellAccessor<2> >                  quad_iterator;
+    typedef TriaActiveIterator<2,MGDoFCellAccessor<2> >            active_quad_iterator;
+
+    typedef raw_quad_iterator    raw_cell_iterator;
+    typedef quad_iterator        cell_iterator;
+    typedef active_quad_iterator active_cell_iterator;
+
+    typedef raw_line_iterator    raw_face_iterator;
+    typedef line_iterator        face_iterator;
+    typedef active_line_iterator active_face_iterator;    
+};
 
 
 
@@ -22,7 +91,26 @@ template <int dim, typename BaseClass> class MGDoFQuadAccessor;
 template <int dim>
 class MGDoFHandler : public DoFHandler<dim> {
   public:
-                                    /**
+                                    // insert these definitions for gcc2.8,
+                                    // since it can't inherit typedefs (I
+                                    // believe it should, but it can't)
+    typedef typename MGDoFDimensionInfo<dim>::raw_line_iterator raw_line_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::line_iterator line_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_line_iterator active_line_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_quad_iterator raw_quad_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
+
+    typedef typename MGDoFDimensionInfo<dim>::raw_face_iterator raw_face_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::face_iterator face_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_face_iterator active_face_iterator;
+
+                                    /**
                                      * Constructor. Take #tria# as the
                                      * triangulation to work on.
                                      */
@@ -69,6 +157,473 @@ class MGDoFHandler : public DoFHandler<dim> {
                                const vector<int> &starting_points = vector<int>());
 
 
+                                    /*--------------------------------------*/
+    
+                                    /**
+                                     *  @name Cell iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first cell, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no cells, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 1D and #begin_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    begin_raw   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used cell
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 1D and #begin_quad# in 2D.
+                                     */
+    cell_iterator        begin       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  cell on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 1D and #begin_active_quad# in 2D.
+                                     */
+    active_cell_iterator begin_active(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 1D and #end_quad# in 2D.
+                                     */
+    raw_cell_iterator    end () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    cell_iterator        end (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_cell_iterator    end_raw (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_cell_iterator end_active (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last cell, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 1D and #last_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    last_raw () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  cell of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 1D and #last_raw_quad# in 2D.
+                                     */
+    raw_cell_iterator    last_raw (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used cell.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 1D and #last_quad# in 2D.
+                                     */
+    cell_iterator        last () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used cell on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 1D and #last_quad# in 2D.
+                                     */
+    cell_iterator        last (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active cell.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 1D and #last_active_quad# in 2D.
+                                     */
+    active_cell_iterator last_active () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active cell on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 1D and #last_active_quad# in 2D.
+                                     */
+    active_cell_iterator last_active (const unsigned int level) const;
+                                    //@}
+
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Face iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first face, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no faces, a past-the-end iterator
+                                     *  is returned.
+                                     *
+                                     *  This function calls #begin_raw_line#
+                                     *  in 2D and #begin_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    begin_raw_face   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used face
+                                     *  on level #level#.
+                                     *
+                                     *  This function calls #begin_line#
+                                     *  in 2D and #begin_quad# in 3D.
+                                     */
+    face_iterator        begin_face       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  face on level #level#.
+                                     *
+                                     *  This function calls #begin_active_line#
+                                     *  in 2D and #begin_active_quad# in 3D.
+                                     */
+    active_face_iterator begin_active_face(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     *
+                                     *  This function calls #end_line#
+                                     *  in 2D and #end_quad# in 3D.
+                                     */
+    raw_face_iterator    end_face () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    face_iterator        end_face (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_face_iterator    end_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_face_iterator end_active_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last face, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  face of the level #level#, used or not.
+                                     *
+                                     *  This function calls #last_raw_line#
+                                     *  in 2D and #last_raw_quad# in 3D.
+                                     */
+    raw_face_iterator    last_raw_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used face on level #level#.
+                                     *
+                                     *  This function calls #last_line#
+                                     *  in 2D and #last_quad# in 3D.
+                                     */
+    face_iterator        last_face (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active face on level #level#.
+                                     *
+                                     *  This function calls #last_active_line#
+                                     *  in 2D and #last_active_quad# in 3D.
+                                     */
+    active_face_iterator last_active_face (const unsigned int level) const;
+                                    //@}
+
+    
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Line iterator functions
+                                     */
+                                    /*@{*/
+                                    /**
+                                     *  Return iterator to the first line, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no lines, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_line_iterator
+    begin_raw_line   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used line
+                                     *  on level #level#.
+                                     */
+    line_iterator
+    begin_line       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  line on level #level#.
+                                     */
+    active_line_iterator
+    begin_active_line(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_line_iterator
+    end_line () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    line_iterator        end_line (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_line_iterator    end_raw_line (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_line_iterator end_active_line (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last line, used or not.
+                                     */
+    raw_line_iterator
+    last_raw_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  line of the level #level#, used or not.
+
+                                     */
+    raw_line_iterator
+    last_raw_line (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used line.
+                                     */
+    line_iterator
+    last_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used line on level #level#.
+                                     */
+    line_iterator
+    last_line (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active line.
+                                     */
+    active_line_iterator
+    last_active_line () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active line on level #level#.
+                                     */
+    active_line_iterator
+    last_active_line (const unsigned int level) const;
+                                    /*@}*/       
+
+                                    /*---------------------------------------*/
+
+                                    /**
+                                     *  @name Quad iterator functions*/
+                                    /*@{
+                                     */
+                                    /**
+                                     *  Return iterator to the first quad, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no quads, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_quad_iterator
+    begin_raw_quad   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used quad
+                                     *  on level #level#.
+                                     */
+    quad_iterator
+    begin_quad       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  quad on level #level#.
+                                     */
+    active_quad_iterator
+    begin_active_quad(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_quad_iterator
+    end_quad () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    quad_iterator        end_quad (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_quad_iterator    end_raw_quad (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_quad_iterator end_active_quad (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last quad, used or not.
+                                     */
+    raw_quad_iterator
+    last_raw_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  quad of the level #level#, used or not.
+
+                                     */
+    raw_quad_iterator
+    last_raw_quad (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used quad.
+                                     */
+    quad_iterator
+    last_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used quad on level #level#.
+                                     */
+    quad_iterator
+    last_quad (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active quad.
+                                     */
+    active_quad_iterator
+    last_active_quad () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active quad on level #level#.
+                                     */
+    active_quad_iterator
+    last_active_quad (const unsigned int level) const;
+                                    /*@}*/
+
+                                    /*---------------------------------------*/
+    
   private:
 
                                     /** We need for each vertex a list of the
@@ -93,12 +648,23 @@ class MGDoFHandler : public DoFHandler<dim> {
       public:
        
                                         /**
-                                         * Constructor. Allocates memory and
-                                         * sets all indices to #-1#.
+                                         * Constructor. This one is empty
+                                         * because it is difficult to make it
+                                         * efficient to use vector<>'s and
+                                         * still construct the object using
+                                         * the constructor. Use the #init#
+                                         * function to really allocate
+                                         * memory.
                                          */
-       MGVertexDoFs (const unsigned int coarsest_level,
-                     const unsigned int n_levels,
-                     const unsigned int dofs_per_vertex);
+       MGVertexDoFs ();
+
+                                        /**
+                                         * Allocate memory and
+                                         * set all indices to #-1#.
+                                         */
+       void init (const unsigned int coarsest_level,
+                  const unsigned int n_levels,
+                  const unsigned int dofs_per_vertex);
 
                                         /**
                                          * Destructor
@@ -148,6 +714,7 @@ class MGDoFHandler : public DoFHandler<dim> {
        DeclException0 (ExcNoMemory);
                                         /**
                                          * Exception.
+
                                          */
        DeclException0 (ExcInvalidIndex);
                                         /**
@@ -157,6 +724,10 @@ class MGDoFHandler : public DoFHandler<dim> {
                        int,
                        << "The given level number " << arg1 << " is below the "
                        << "coarsest level this vertex lives on.");
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcInternalError);
        
       private:
                                         /**
@@ -165,12 +736,12 @@ class MGDoFHandler : public DoFHandler<dim> {
                                          * as an offset when accessing the
                                          * dofs of a specific level.
                                          */
-       const unsigned int coarsest_level;
+       unsigned int coarsest_level;
 
                                         /**
                                          * Array holding the indices.
                                          */
-       int *const         indices;
+       int         *indices;
     };
 
     
@@ -180,6 +751,23 @@ class MGDoFHandler : public DoFHandler<dim> {
                                      */
     void reserve_space ();
 
+                                    /**
+                                     * Distribute dofs on the given cell,
+                                     * with new dofs starting with index
+                                     * #next_free_dof#. Return the next
+                                     * unused index number. The finite
+                                     * element used is the one given to
+                                     * #distribute_dofs#, which is copied
+                                     * to #selected_fe#.
+                                     *
+                                     * This function is excluded from the
+                                     * #distribute_dofs# function since
+                                     * it can not be implemented dimension
+                                     * independent.
+                                     */
+    unsigned int distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int next_free_dof);
+
                                     /**
                                      * Space to store the DoF numbers for the
                                      * different levels. Unlike the #levels#
index 322b395cfe6b9d8e17d3228ce9062dad5d0d637e..d65ad36ca2ef22681a3c7182dc41172142fb6b33 100644 (file)
@@ -23,11 +23,9 @@ template <int dim, typename BaseClass>
 MGDoFLineAccessor<dim,BaseClass>::MGDoFLineAccessor (Triangulation<dim> *tria,
                                                     const int           level,
                                                     const int           index,
-                                                    const void         *local_data) :
-               MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
-               DoFLineAccessor(tria,level,index,
-                               static_cast<DoFHandler<dim>*>
-                               (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+                                                    const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFLineAccessor(tria,level,index,local_data) {};
 
 
 
@@ -139,6 +137,15 @@ MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
 
 
 
+template <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::copy_from (const MGDoFLineAccessor<dim,BaseClass> &a) {
+  DoFLineAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
 
 /* ------------------------ MGDoFQuadAccessor --------------------------- */
 
@@ -146,11 +153,9 @@ template <int dim, typename BaseClass>
 MGDoFQuadAccessor<dim,BaseClass>::MGDoFQuadAccessor (Triangulation<dim> *tria,
                                                     const int           level,
                                                     const int           index,
-                                                    const void         *local_data) :
-               MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
-               DoFQuadAccessor(tria,level,index,
-                               static_cast<DoFHandler<dim>*>
-                               (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+                                                    const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFQuadAccessor(tria,level,index,local_data) {};
 
 
 
@@ -282,6 +287,15 @@ MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
 
 
 
+template <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::copy_from (const MGDoFQuadAccessor<dim,BaseClass> &a) {
+  DoFQuadAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
 
 
 /*------------------------- Functions: MGDoFCellAccessor -----------------------*/
index 14b9a55a01e6676ed2307c1696c3ff14cbd0deca..bb3d7d999fd6107a4c5b22c2aad1a998c99c4252 100644 (file)
@@ -2,7 +2,7 @@
 
 #include <grid/mg_dof.h>
 #include <grid/dof_levels.h>
-#include <grid/dof_accessor.h>
+#include <grid/mg_dof_accessor.h>
 #include <grid/dof_constraints.h>
 #include <grid/tria_levels.h>
 #include <grid/tria_accessor.h>
 /* ------------------------ MGVertexDoFs ----------------------------------- */
 
 template <int dim>
-MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level,
-                                              const unsigned int n_levels,
-                                              const unsigned int dofs_per_vertex) :
-               coarsest_level (coarsest_level),
-               indices (new int[n_levels * dofs_per_vertex])
-{
+MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs () :
+               coarsest_level (1<<30),
+               indices (0)
+{};
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::MGVertexDoFs::init (const unsigned int cl,
+                                           const unsigned int n_levels,
+                                           const unsigned int dofs_per_vertex) {
+  Assert (indices == 0, ExcInternalError());
+
+  coarsest_level = cl;
+  
+  indices = new int[n_levels * dofs_per_vertex];
   Assert (indices != 0, ExcNoMemory ());
 
   for (unsigned int i=0; i<n_levels*dofs_per_vertex; ++i)
@@ -32,6 +42,8 @@ MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level
 
 template <int dim>
 MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+  Assert (indices != 0, ExcInternalError ());
+  
   delete[] indices;
 };
 
@@ -56,6 +68,766 @@ MGDoFHandler<dim>::~MGDoFHandler () {};
 
 
 
+
+
+#if deal_II_dimension == 1
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::begin_raw (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::begin (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::begin_active (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::end () const {
+  return end_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::last_raw () const {
+  return last_raw_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::last_raw (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::last () const {
+  return last_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::last (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::last_active () const {
+  return last_active_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::last_active (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::begin_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::begin_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::begin_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::end_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::last_raw_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::last_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::last_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::last_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::last_active_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::last_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::begin_raw_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::begin_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::begin_active_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::end_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::last_raw_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::last_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::last_active_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::last_raw_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::last_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::last_active_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::begin_raw (const unsigned int level) const {
+  return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::begin (const unsigned int level) const {
+  return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::begin_active (const unsigned int level) const {
+  return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::end () const {
+  return end_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::last_raw () const {
+  return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::last_raw (const unsigned int level) const {
+  return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::last () const {
+  return last_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::last (const unsigned int level) const {
+  return last_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::last_active () const {
+  return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::last_active (const unsigned int level) const {
+  return last_active_quad (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::begin_face (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::begin_active_face (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::end_face () const {
+  return end_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::last_raw_face () const {
+  return last_raw_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::last_raw_face (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::last_face () const {
+  return last_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::last_face (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::last_active_face () const {
+  return last_active_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::last_active_face (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+#endif
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::begin_raw_line (const unsigned int level) const {
+  return raw_line_iterator (tria,
+                           tria->begin_raw_line(level)->level(),
+                           tria->begin_raw_line(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::begin_line (const unsigned int level) const {
+  return line_iterator (tria,
+                       tria->begin_line(level)->level(),
+                       tria->begin_line(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::begin_active_line (const unsigned int level) const {
+  return active_line_iterator (tria,
+                              tria->begin_active_line(level)->level(),
+                              tria->begin_active_line(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::begin_raw_quad (const unsigned int level) const {
+  return raw_quad_iterator (tria,
+                           tria->begin_raw_quad(level)->level(),
+                           tria->begin_raw_quad(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::begin_quad (const unsigned int level) const {
+  return quad_iterator (tria,
+                       tria->begin_quad(level)->level(),
+                       tria->begin_quad(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::begin_active_quad (const unsigned int level) const {
+  return active_quad_iterator (tria,
+                              tria->begin_active_quad(level)->level(),
+                              tria->begin_active_quad(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::end_line () const {
+  return raw_line_iterator (tria, -1, -1, this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_quad () const {
+  return raw_quad_iterator (tria, -1, -1, this);
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_cell_iterator
+MGDoFHandler<dim>::end_raw (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end() :
+         begin_raw (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::cell_iterator
+MGDoFHandler<dim>::end (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         cell_iterator(end()) :
+         begin (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_cell_iterator
+MGDoFHandler<dim>::end_active (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_cell_iterator(end()) :
+         begin_active (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_face_iterator
+MGDoFHandler<dim>::end_raw_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_face() :
+         begin_raw_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::face_iterator
+MGDoFHandler<dim>::end_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         face_iterator(end_face()) :
+         begin_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_face_iterator
+MGDoFHandler<dim>::end_active_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_face_iterator(end_face()) :
+         begin_active_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_line_iterator
+MGDoFHandler<dim>::end_raw_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_line() :
+         begin_raw_line (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::line_iterator
+MGDoFHandler<dim>::end_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         line_iterator(end_line()) :
+         begin_line (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_line_iterator
+MGDoFHandler<dim>::end_active_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_line_iterator(end_line()) :
+         begin_active_line (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_quad() :
+         begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         quad_iterator(end_quad()) :
+         begin_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_quad_iterator(end_quad()) :
+         begin_active_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::last_raw_line (const unsigned int level) const {
+  return raw_line_iterator (tria,
+                           tria->last_raw_line(level)->level(),
+                           tria->last_raw_line(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::last_line (const unsigned int level) const {
+  return line_iterator (tria,
+                       tria->last_line(level)->level(),
+                       tria->last_line(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::last_active_line (const unsigned int level) const {
+  return active_line_iterator (tria,
+                              tria->last_active_line(level)->level(),
+                              tria->last_active_line(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::last_raw_quad (const unsigned int level) const {
+  return raw_quad_iterator (tria,
+                           tria->last_raw_quad(level)->level(),
+                           tria->last_raw_quad(level)->index(),
+                           this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::last_quad (const unsigned int level) const {
+  return quad_iterator (tria,
+                       tria->last_quad(level)->level(),
+                       tria->last_quad(level)->index(),
+                       this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::last_active_quad (const unsigned int level) const {
+  return active_quad_iterator (tria,
+                              tria->last_active_quad(level)->level(),
+                              tria->last_active_quad(level)->index(),
+                              this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::last_raw_line () const {
+  return last_raw_line (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::last_raw_quad () const {
+  return last_raw_quad (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::last_line () const {
+  return last_line (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::last_quad () const {
+  return last_quad (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::last_active_line () const {
+  return last_active_line (mg_levels.size()-1);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::last_active_quad () const {
+  return last_active_quad (mg_levels.size()-1);
+};
+
+
+
+
+
+
+//------------------------------------------------------------------
+
+
+
+
+
+
 template <int dim>
 void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
                                   // first distribute global dofs
@@ -64,9 +836,131 @@ void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
 
                                   // reserve space for the MG dof numbers
   reserve_space ();
+
+                                  // clear user flags because we will
+                                  // need them
+  tria->clear_user_flags ();
+
+                                  // now distribute indices on each level
+                                  // separately
+  for (unsigned int level=0; level<tria->n_levels(); ++level)
+    {
+      unsigned int next_free_dof = 0;   
+      active_cell_iterator cell = begin(level),
+                          endc = end(level);
+
+      for (; cell != endc; ++cell) 
+       next_free_dof = distribute_dofs_on_cell (cell, next_free_dof);
+  
+      mg_used_dofs[level] = next_free_dof;
+    };
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+unsigned int
+MGDoFHandler<1>::distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int          next_free_dof) {
+
+                                  // distribute dofs of vertices
+  if (selected_fe->dofs_per_vertex > 0)
+    for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+      {
+       cell_iterator neighbor = cell->neighbor(v);
+       
+       if (neighbor.state() == valid)
+         {
+                                            // has neighbor already been processed?
+           if (neighbor->user_flag_set() &&
+               (neighbor->level() == cell->level()))
+                                              // copy dofs if the neighbor is on
+                                              // the same level (only then are
+                                              // mg dofs the same)
+             {
+               if (v==0) 
+                 for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+                   cell->set_mg_vertex_dof_index (0, d,
+                                                  neighbor->mg_vertex_dof_index (1, d));
+               else
+                 for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+                   cell->set_mg_vertex_dof_index (1, d,
+                                                  neighbor->mg_vertex_dof_index (0, d));
+               
+                                                // next neighbor
+               continue;
+             };
+         };
+       
+                                        // otherwise: create dofs newly
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (v, d, next_free_dof++);
+      };
+  
+                                  // dofs of line
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
+  
+  return next_free_dof;
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+unsigned int
+MGDoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int          next_free_dof) {
+  if (selected_fe->dofs_per_vertex > 0)
+                                    // number dofs on vertices
+    for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
+                                      // check whether dofs for this
+                                      // vertex have been distributed
+                                      // (only check the first dof)
+      if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+    
+                                  // for the four sides
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
+      {
+       line_iterator line = cell->line(side);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (line->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           line->set_mg_dof_index (d, next_free_dof++);            
+      };
+  
+
+                                      // dofs of quad
+  if (selected_fe->dofs_per_quad > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_quad; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+  
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
   
+  return next_free_dof;
 };
 
+#endif
+
+
+
 
 
 
@@ -74,6 +968,81 @@ void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
 
 template <>
 void MGDoFHandler<1>::reserve_space () {
+  const unsigned int dim = 1;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+
+                                  //////////////////////////
+                                  // DESTRUCTION
+  
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.resize (0);
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.resize (0);
+
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<1>);
+
+      mg_levels.back()->line_dofs = vector<int>(tria->levels[i]->lines.lines.size() *
+                                               selected_fe->dofs_per_line,
+                                               -1);
+    };
+
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex]-min_level[vertex]+1,
+                                  selected_fe->dofs_per_vertex);
+    };
 };
 
 #endif
@@ -83,6 +1052,85 @@ void MGDoFHandler<1>::reserve_space () {
 
 template <>
 void MGDoFHandler<2>::reserve_space () {
+  const unsigned int dim = 2;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+  
+                                  ////////////////////////////
+                                  // DESTRUCTION
+
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.resize (0);
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.resize (0);
+
+
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines and quads on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<2>);
+
+      mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+                                                selected_fe->dofs_per_line,
+                                                -1);
+      mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+                                                selected_fe->dofs_per_quad,
+                                                -1);
+    };
+
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex]-min_level[vertex]+1,
+                                  selected_fe->dofs_per_vertex);
+    };
 };
 
 #endif
index 322b395cfe6b9d8e17d3228ce9062dad5d0d637e..d65ad36ca2ef22681a3c7182dc41172142fb6b33 100644 (file)
@@ -23,11 +23,9 @@ template <int dim, typename BaseClass>
 MGDoFLineAccessor<dim,BaseClass>::MGDoFLineAccessor (Triangulation<dim> *tria,
                                                     const int           level,
                                                     const int           index,
-                                                    const void         *local_data) :
-               MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
-               DoFLineAccessor(tria,level,index,
-                               static_cast<DoFHandler<dim>*>
-                               (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+                                                    const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFLineAccessor(tria,level,index,local_data) {};
 
 
 
@@ -139,6 +137,15 @@ MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
 
 
 
+template <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::copy_from (const MGDoFLineAccessor<dim,BaseClass> &a) {
+  DoFLineAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
 
 /* ------------------------ MGDoFQuadAccessor --------------------------- */
 
@@ -146,11 +153,9 @@ template <int dim, typename BaseClass>
 MGDoFQuadAccessor<dim,BaseClass>::MGDoFQuadAccessor (Triangulation<dim> *tria,
                                                     const int           level,
                                                     const int           index,
-                                                    const void         *local_data) :
-               MGDoFAccessor<dim> ((MGDoFHandler<dim>*)local_data),
-               DoFQuadAccessor(tria,level,index,
-                               static_cast<DoFHandler<dim>*>
-                               (reinterpret_cast<MGDoFHandler<dim>*>(local_data))) {};
+                                                    const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFQuadAccessor(tria,level,index,local_data) {};
 
 
 
@@ -282,6 +287,15 @@ MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
 
 
 
+template <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::copy_from (const MGDoFQuadAccessor<dim,BaseClass> &a) {
+  DoFQuadAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
 
 
 /*------------------------- Functions: MGDoFCellAccessor -----------------------*/
index 14b9a55a01e6676ed2307c1696c3ff14cbd0deca..bb3d7d999fd6107a4c5b22c2aad1a998c99c4252 100644 (file)
@@ -2,7 +2,7 @@
 
 #include <grid/mg_dof.h>
 #include <grid/dof_levels.h>
-#include <grid/dof_accessor.h>
+#include <grid/mg_dof_accessor.h>
 #include <grid/dof_constraints.h>
 #include <grid/tria_levels.h>
 #include <grid/tria_accessor.h>
 /* ------------------------ MGVertexDoFs ----------------------------------- */
 
 template <int dim>
-MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level,
-                                              const unsigned int n_levels,
-                                              const unsigned int dofs_per_vertex) :
-               coarsest_level (coarsest_level),
-               indices (new int[n_levels * dofs_per_vertex])
-{
+MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs () :
+               coarsest_level (1<<30),
+               indices (0)
+{};
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::MGVertexDoFs::init (const unsigned int cl,
+                                           const unsigned int n_levels,
+                                           const unsigned int dofs_per_vertex) {
+  Assert (indices == 0, ExcInternalError());
+
+  coarsest_level = cl;
+  
+  indices = new int[n_levels * dofs_per_vertex];
   Assert (indices != 0, ExcNoMemory ());
 
   for (unsigned int i=0; i<n_levels*dofs_per_vertex; ++i)
@@ -32,6 +42,8 @@ MGDoFHandler<dim>::MGVertexDoFs::MGVertexDoFs (const unsigned int coarsest_level
 
 template <int dim>
 MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+  Assert (indices != 0, ExcInternalError ());
+  
   delete[] indices;
 };
 
@@ -56,6 +68,766 @@ MGDoFHandler<dim>::~MGDoFHandler () {};
 
 
 
+
+
+#if deal_II_dimension == 1
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::begin_raw (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::begin (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::begin_active (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::end () const {
+  return end_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::last_raw () const {
+  return last_raw_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_cell_iterator
+MGDoFHandler<1>::last_raw (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::last () const {
+  return last_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::cell_iterator
+MGDoFHandler<1>::last (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::last_active () const {
+  return last_active_line ();
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_cell_iterator
+MGDoFHandler<1>::last_active (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::begin_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::begin_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::begin_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::end_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::last_raw_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::raw_face_iterator
+MGDoFHandler<1>::last_raw_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::last_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::face_iterator
+MGDoFHandler<1>::last_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::last_active_face () const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+template <>
+MGDoFDimensionInfo<1>::active_face_iterator
+MGDoFHandler<1>::last_active_face (const unsigned int) const {
+  Assert (false, ExcFunctionNotUseful());
+  return 0;
+};
+
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::begin_raw_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::begin_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::begin_active_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::end_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::last_raw_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::last_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::last_active_quad (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_quad_iterator
+MGDoFHandler<1>::last_raw_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::quad_iterator
+MGDoFHandler<1>::last_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_quad_iterator
+MGDoFHandler<1>::last_active_quad () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::begin_raw (const unsigned int level) const {
+  return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::begin (const unsigned int level) const {
+  return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::begin_active (const unsigned int level) const {
+  return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::end () const {
+  return end_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::last_raw () const {
+  return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_cell_iterator
+MGDoFHandler<2>::last_raw (const unsigned int level) const {
+  return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::last () const {
+  return last_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::cell_iterator
+MGDoFHandler<2>::last (const unsigned int level) const {
+  return last_quad (level);
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::last_active () const {
+  return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_cell_iterator
+MGDoFHandler<2>::last_active (const unsigned int level) const {
+  return last_active_quad (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::begin_face (const unsigned int level) const {
+  return begin_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::begin_active_face (const unsigned int level) const {
+  return begin_active_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::end_face () const {
+  return end_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::last_raw_face () const {
+  return last_raw_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::raw_face_iterator
+MGDoFHandler<2>::last_raw_face (const unsigned int level) const {
+  return last_raw_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::last_face () const {
+  return last_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::face_iterator
+MGDoFHandler<2>::last_face (const unsigned int level) const {
+  return last_line (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::last_active_face () const {
+  return last_active_line ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<2>::active_face_iterator
+MGDoFHandler<2>::last_active_face (const unsigned int level) const {
+  return last_active_line (level);
+};
+
+#endif
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::begin_raw_line (const unsigned int level) const {
+  return raw_line_iterator (tria,
+                           tria->begin_raw_line(level)->level(),
+                           tria->begin_raw_line(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::begin_line (const unsigned int level) const {
+  return line_iterator (tria,
+                       tria->begin_line(level)->level(),
+                       tria->begin_line(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::begin_active_line (const unsigned int level) const {
+  return active_line_iterator (tria,
+                              tria->begin_active_line(level)->level(),
+                              tria->begin_active_line(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::begin_raw_quad (const unsigned int level) const {
+  return raw_quad_iterator (tria,
+                           tria->begin_raw_quad(level)->level(),
+                           tria->begin_raw_quad(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::begin_quad (const unsigned int level) const {
+  return quad_iterator (tria,
+                       tria->begin_quad(level)->level(),
+                       tria->begin_quad(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::begin_active_quad (const unsigned int level) const {
+  return active_quad_iterator (tria,
+                              tria->begin_active_quad(level)->level(),
+                              tria->begin_active_quad(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::end_line () const {
+  return raw_line_iterator (tria, -1, -1, this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_quad () const {
+  return raw_quad_iterator (tria, -1, -1, this);
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_cell_iterator
+MGDoFHandler<dim>::end_raw (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end() :
+         begin_raw (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::cell_iterator
+MGDoFHandler<dim>::end (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         cell_iterator(end()) :
+         begin (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_cell_iterator
+MGDoFHandler<dim>::end_active (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_cell_iterator(end()) :
+         begin_active (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_face_iterator
+MGDoFHandler<dim>::end_raw_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_face() :
+         begin_raw_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::face_iterator
+MGDoFHandler<dim>::end_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         face_iterator(end_face()) :
+         begin_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_face_iterator
+MGDoFHandler<dim>::end_active_face (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_face_iterator(end_face()) :
+         begin_active_face (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_line_iterator
+MGDoFHandler<dim>::end_raw_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_line() :
+         begin_raw_line (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::line_iterator
+MGDoFHandler<dim>::end_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         line_iterator(end_line()) :
+         begin_line (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_line_iterator
+MGDoFHandler<dim>::end_active_line (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_line_iterator(end_line()) :
+         begin_active_line (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_quad() :
+         begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         quad_iterator(end_quad()) :
+         begin_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_quad_iterator(end_quad()) :
+         begin_active_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::last_raw_line (const unsigned int level) const {
+  return raw_line_iterator (tria,
+                           tria->last_raw_line(level)->level(),
+                           tria->last_raw_line(level)->index(),
+                           this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::last_line (const unsigned int level) const {
+  return line_iterator (tria,
+                       tria->last_line(level)->level(),
+                       tria->last_line(level)->index(),
+                       this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::last_active_line (const unsigned int level) const {
+  return active_line_iterator (tria,
+                              tria->last_active_line(level)->level(),
+                              tria->last_active_line(level)->index(),
+                              this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::last_raw_quad (const unsigned int level) const {
+  return raw_quad_iterator (tria,
+                           tria->last_raw_quad(level)->level(),
+                           tria->last_raw_quad(level)->index(),
+                           this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::last_quad (const unsigned int level) const {
+  return quad_iterator (tria,
+                       tria->last_quad(level)->level(),
+                       tria->last_quad(level)->index(),
+                       this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::last_active_quad (const unsigned int level) const {
+  return active_quad_iterator (tria,
+                              tria->last_active_quad(level)->level(),
+                              tria->last_active_quad(level)->index(),
+                              this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_line_iterator
+MGDoFHandler<dim>::last_raw_line () const {
+  return last_raw_line (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::raw_quad_iterator
+MGDoFHandler<dim>::last_raw_quad () const {
+  return last_raw_quad (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::line_iterator
+MGDoFHandler<dim>::last_line () const {
+  return last_line (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::quad_iterator
+MGDoFHandler<dim>::last_quad () const {
+  return last_quad (mg_levels.size()-1);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_line_iterator
+MGDoFHandler<dim>::last_active_line () const {
+  return last_active_line (mg_levels.size()-1);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_quad_iterator
+MGDoFHandler<dim>::last_active_quad () const {
+  return last_active_quad (mg_levels.size()-1);
+};
+
+
+
+
+
+
+//------------------------------------------------------------------
+
+
+
+
+
+
 template <int dim>
 void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
                                   // first distribute global dofs
@@ -64,9 +836,131 @@ void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
 
                                   // reserve space for the MG dof numbers
   reserve_space ();
+
+                                  // clear user flags because we will
+                                  // need them
+  tria->clear_user_flags ();
+
+                                  // now distribute indices on each level
+                                  // separately
+  for (unsigned int level=0; level<tria->n_levels(); ++level)
+    {
+      unsigned int next_free_dof = 0;   
+      active_cell_iterator cell = begin(level),
+                          endc = end(level);
+
+      for (; cell != endc; ++cell) 
+       next_free_dof = distribute_dofs_on_cell (cell, next_free_dof);
+  
+      mg_used_dofs[level] = next_free_dof;
+    };
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+unsigned int
+MGDoFHandler<1>::distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int          next_free_dof) {
+
+                                  // distribute dofs of vertices
+  if (selected_fe->dofs_per_vertex > 0)
+    for (unsigned int v=0; v<GeometryInfo<1>::vertices_per_cell; ++v)
+      {
+       cell_iterator neighbor = cell->neighbor(v);
+       
+       if (neighbor.state() == valid)
+         {
+                                            // has neighbor already been processed?
+           if (neighbor->user_flag_set() &&
+               (neighbor->level() == cell->level()))
+                                              // copy dofs if the neighbor is on
+                                              // the same level (only then are
+                                              // mg dofs the same)
+             {
+               if (v==0) 
+                 for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+                   cell->set_mg_vertex_dof_index (0, d,
+                                                  neighbor->mg_vertex_dof_index (1, d));
+               else
+                 for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+                   cell->set_mg_vertex_dof_index (1, d,
+                                                  neighbor->mg_vertex_dof_index (0, d));
+               
+                                                // next neighbor
+               continue;
+             };
+         };
+       
+                                        // otherwise: create dofs newly
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (v, d, next_free_dof++);
+      };
+  
+                                  // dofs of line
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
+  
+  return next_free_dof;
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+unsigned int
+MGDoFHandler<2>::distribute_dofs_on_cell (active_cell_iterator &cell,
+                                         unsigned int          next_free_dof) {
+  if (selected_fe->dofs_per_vertex > 0)
+                                    // number dofs on vertices
+    for (unsigned int vertex=0; vertex<GeometryInfo<2>::vertices_per_cell; ++vertex)
+                                      // check whether dofs for this
+                                      // vertex have been distributed
+                                      // (only check the first dof)
+      if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+    
+                                  // for the four sides
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int side=0; side<GeometryInfo<2>::faces_per_cell; ++side)
+      {
+       line_iterator line = cell->line(side);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (line->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           line->set_mg_dof_index (d, next_free_dof++);            
+      };
+  
+
+                                      // dofs of quad
+  if (selected_fe->dofs_per_quad > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_quad; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+  
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
   
+  return next_free_dof;
 };
 
+#endif
+
+
+
 
 
 
@@ -74,6 +968,81 @@ void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
 
 template <>
 void MGDoFHandler<1>::reserve_space () {
+  const unsigned int dim = 1;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+
+                                  //////////////////////////
+                                  // DESTRUCTION
+  
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.resize (0);
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.resize (0);
+
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<1>);
+
+      mg_levels.back()->line_dofs = vector<int>(tria->levels[i]->lines.lines.size() *
+                                               selected_fe->dofs_per_line,
+                                               -1);
+    };
+
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex]-min_level[vertex]+1,
+                                  selected_fe->dofs_per_vertex);
+    };
 };
 
 #endif
@@ -83,6 +1052,85 @@ void MGDoFHandler<1>::reserve_space () {
 
 template <>
 void MGDoFHandler<2>::reserve_space () {
+  const unsigned int dim = 2;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+  
+                                  ////////////////////////////
+                                  // DESTRUCTION
+
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.resize (0);
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.resize (0);
+
+
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines and quads on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<2>);
+
+      mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+                                                selected_fe->dofs_per_line,
+                                                -1);
+      mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+                                                selected_fe->dofs_per_quad,
+                                                -1);
+    };
+
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex]-min_level[vertex]+1,
+                                  selected_fe->dofs_per_vertex);
+    };
 };
 
 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.