]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Initial version of the multigrid components.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Aug 1998 16:19:05 +0000 (16:19 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 25 Aug 1998 16:19:05 +0000 (16:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@518 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/mg_dof_accessor.h [new file with mode: 0644]
deal.II/deal.II/include/dofs/mg_dof_handler.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/mg_dof_accessor.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/mg_dof_handler.h [new file with mode: 0644]
deal.II/deal.II/source/dofs/mg_dof_accessor.cc [new file with mode: 0644]
deal.II/deal.II/source/dofs/mg_dof_handler.cc [new file with mode: 0644]
deal.II/deal.II/source/multigrid/mg_dof_accessor.cc [new file with mode: 0644]
deal.II/deal.II/source/multigrid/mg_dof_handler.cc [new file with mode: 0644]

diff --git a/deal.II/deal.II/include/dofs/mg_dof_accessor.h b/deal.II/deal.II/include/dofs/mg_dof_accessor.h
new file mode 100644 (file)
index 0000000..0327475
--- /dev/null
@@ -0,0 +1,505 @@
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __mg_dof_accessor_H
+#define __mg_dof_accessor_H
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+
+
+#include <grid/dof_accessor.h>
+
+
+// forward declaration
+template <int dim> class MGDoFHandler;
+
+
+
+
+/**
+ * Define the basis for accessors to the degrees of freedom for
+ * a multigrid DoF handler object.
+ *
+ * Note that it is allowed to construct an object of which the
+ * #mg_dof_handler# pointer is a Null pointer. Such an object would
+ * result in a strange kind of behaviour, though every reasonable
+ * operating system should disallow access through that pointer.
+ * The reason we do not check for the null pointer in the
+ * constructor which gets passed the #DoFHandler# pointer is that
+ * if we did we could not make dof iterators member of other classes
+ * (like in the #FEValues# class) if we did not know about the
+ * #DoFHandler# object to be used upon construction of that object.
+ * Through the way this class is implemented here, we allow the
+ * creation of a kind of virgin object which only gets useful if
+ * assigned to from another object before first usage.
+ *
+ * Opposite to construction, it is not possible to copy an object
+ * which has an invalid dof handler pointer. This is to guarantee
+ * that every iterator which is once assigned to is a valid
+ * object. However, this assertion only holds in debug mode, when
+ * the #Assert# macro is switched on.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFAccessor {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFAccessor () : mg_dof_handler(0) {
+      Assert (false, ExcInvalidObject());
+    };
+
+                                    /**
+                                     * This should be the default constructor.
+                                     */
+    MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
+                   mg_dof_handler(mg_dof_handler) {};
+
+                                    /**
+                                     * Reset the DoF handler pointer.
+                                     */
+    void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
+      Assert (dh != 0, ExcInvalidObject());
+      mg_dof_handler = dh;
+    };
+
+                                    /**
+                                     * Copy operator.
+                                     */
+    MGDoFAccessor<dim> & operator = (const MGDoFAccessor<dim> &da) {
+      set_dof_handler (da.mg_dof_handler);
+      return *this;
+    };
+
+  protected:
+                                    /**
+                                     * Store the address of the #MGDoFHandler# object
+                                     * to be accessed.
+                                     */
+    MGDoFHandler<dim> *mg_dof_handler;  
+};
+
+
+
+
+
+/**
+ * Grant access to the degrees of freedom located on lines.
+ * This class follows mainly the route laid out by the accessor library
+ * declared in the triangulation library (\Ref{TriaAccessor}). It enables
+ * the user to access the degrees of freedom on the lines (there are similar
+ * versions for the DoFs on quads, etc), where the dimension of the underlying
+ * triangulation does not really matter (i.e. this accessor works with the
+ * lines in 1D-, 2D-, etc dimensions).
+ *
+ *
+ * \subsection{Usage}
+ *
+ * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes
+ * declare typedefs to iterators using the accessors declared in this class
+ * hierarchy tree. Usage is best to happens through these typedefs, since they
+ * are more secure to changes in the class naming and template interface as well
+ * as they provide easier typing (much less complicated names!).
+ * 
+ * 
+ * \subsection{Notes about the class hierarchy structure}
+ *
+ * The class hierarchy seems to be a bit confused here. The reason for this is
+ * that we would really like to derive a #DoFLineAccessor# from a #LineAccessor#.
+ * Unfortunately, we would run into problems, if we wanted a #DoFLineAccessor#
+ * in one spatial dimension, in which case a line is also a cell. The traditional
+ * solution would be to declare a #DoFCellAccessor<1># which is derived from
+ * #DoFLineAccessor<1># and #CellAccessor<1># (the #DoFLineAccessor<dim># cannot
+ * itself be derived from #CellAccessor<dim># since a line is not a cell
+ * unless in one space dimension), but since a #DoFLineAccessor# and a
+ * #CellAccessor# are both derived from #TriaAccessor#, we would have to make
+ * the last derivation virtual.
+ *
+ * Since we want to avoid virtual inheritance since this involves another
+ * indirection in every member variable access, we chose another way: we
+ * pass a second template parameter to a #DoFLineAccessor# which tells it
+ * which class to be derived from: if we are in one spatial dimension, the
+ * base class is to be #CellAccessor<1>#, in two or more dimensions it
+ * is a #LineAccessor<dim>#, i.e. am accessor to lines without the missing
+ * functionality needed for cells (neighbors, etc.).
+ *
+ * This way we can declare a #DoFCellAccessor# in one dimension by deriving
+ * from #DoFLineAccessor<1,CellAccessor<1> >#, thus getting the cell
+ * functionality through the #DoFLineAccessor# instead of through a virtual
+ * multiple inheritance of #DoFLineAccessor# and #CellAccessor<1>#.
+ *
+ * The same concept is used with #DoFQuadAccessor# classes etc.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim, typename BaseClass>
+class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
+                          public DoFLineAccessor<dim, BaseClass> {
+  public:
+                                    /**
+                                     * Default constructor, unused thus
+                                     * not implemented.
+                                     */
+    MGDoFLineAccessor ();
+    
+                                    /**
+                                     * Constructor. The #local_data#
+                                     * argument is assumed to be a pointer
+                                     * to a #MGDoFHandler<dim># object.
+                                     */
+    MGDoFLineAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data);
+    
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * of freedom of this line on the level
+                                     * this line lives on.
+                                     */
+    int mg_dof_index (const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * of freedom of this line on the
+                                     * level this line lives on to #index#.
+                                     */
+    void set_mg_dof_index (const unsigned int i, const int index) const;
+
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * on the #vertex#th vertex for the
+                                     * level this line lives on.
+                                     */
+    int mg_vertex_dof_index (const unsigned int vertex,
+                            const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * on the #vertex#th vertex to #index#
+                                     * for the level this line lives on.
+                                     */
+    void set_mg_vertex_dof_index (const unsigned int vertex,
+                                 const unsigned int i,
+                                 const int          index) const;
+
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * line in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, 
+                                     * dofs on line.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the local numbering
+                                     * for the level this line lives on.
+                                     */
+    void get_mg_dof_indices (vector<int> &dof_indices) const;
+
+                                    /**
+                                     * Return the #i#th child as a MGDoF line
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a line accessor without
+                                     * access to the MG data.
+                                     */
+    TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > child (const unsigned int) const;
+    
+                                    /**
+                                     * Implement the copy operator needed
+                                     * for the iterator classes.
+                                     */
+    void copy_from (const MGDoFLineAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+
+/**
+ * Grant access to the multilevel degrees of freedom located on quads.
+ *
+ * @see DoFLineAccessor
+ */
+template <int dim, typename BaseClass>
+class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
+                          public DoFQuadAccessor<dim, BaseClass> {
+  public:
+                                    /**
+                                     * Default constructor, unused thus
+                                     * not implemented.
+                                     */
+    MGDoFQuadAccessor ();
+
+                                    /**
+                                     * Constructor. The #local_data#
+                                     * argument is assumed to be a pointer
+                                     * to a #DoFHandler<dim># object.
+                                     */
+    MGDoFQuadAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data);
+    
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * of freedom of this quad on the level
+                                     * this quad lives on.
+                                     */
+    int mg_dof_index (const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * of freedom of this quad on the
+                                     * level this quad lives on to #index#.
+                                     */
+    void set_mg_dof_index (const unsigned int i, const int index) const;
+
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * on the #vertex#th vertex for the level
+                                     * this quad lives on.
+                                     */
+    int mg_vertex_dof_index (const unsigned int vertex,
+                            const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * on the #vertex#th vertex for the
+                                     * level this quad lives on to #index#.
+                                     */
+    void set_mg_vertex_dof_index (const unsigned int vertex,
+                                 const unsigned int i,
+                                 const int          index) const;
+
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * quad in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the local numbering
+                                     * for the level this quad lives on.
+                                     */
+    void get_mg_dof_indices (vector<int> &dof_indices) const;
+
+                                    /**
+                                     * Return a pointer to the #i#th line
+                                     * bounding this #Quad#.
+                                     */
+    TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+    line (const unsigned int i) const;
+
+                                    /**
+                                     * Return the #i#th child as a DoF quad
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a quad accessor without
+                                     * access to the DoF data.
+                                     */
+    TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > child (const unsigned int) const;
+    
+                                    /**
+                                     * Implement the copy operator needed
+                                     * for the iterator classes.
+                                     */
+    void copy_from (const DoFQuadAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * Rationale for the declaration of members for this class: gcc 2.8 has a bug
+ * when deriving from explicitely specialized classes which materializes in
+ * the calculation of wrong addresses of member variables. By declaring the
+ * general template of #DoFSubstructAccessor# to have the same object layout as
+ * the specialized versions (using the same base classes), we fool the compiler,
+ * which still looks in the wrong place for the addresses but finds the
+ * right information. This way, at least ot works.
+ *
+ * Insert a guard, however, in the constructor to avoid that anyone (including
+ * the compiler) happens to use this class.
+ */
+template <int dim>
+class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
+                              public TriaAccessor<dim> {
+  public:
+    DoFSubstructAccessor () {
+      Assert (false, ExcInternalError());
+    };
+
+    DeclException0 (ExcInternalError);
+};
+
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * \subsection{Rationale}
+ *
+ * This class is only a wrapper class used to do kind of a typedef
+ * 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> >;
+ * \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
+ * the \Ref{DoFCellAccessor} dimension independent as an inheritance from
+ * #DoFSubstructAccessor<dim>#. If we had not declared these
+ * types, we would have to write two class declarations, one for
+ * #DoFCellAccessor<1>#, derived from #DoFLineAccessor<1,CellAccessor<1> >#
+ * and one for #DoFCellAccessor<2>#, derived from
+ * #DoFQuadAccessor<2,CellAccessor<2> >#.
+ */
+template <>
+class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> > {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFSubstructAccessor (Triangulation<1> *tria,
+                           const int         level,
+                           const int         index,
+                           const void       *local_data) :
+                   MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
+                                    // do this here, since this way the
+                                    // CellAccessor has the possibility to know
+                                    // what a face_iterator is. Otherwise
+                                    // it would have to ask the MGDoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef void * face_iterator;
+};
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * @see MGDoFSubstructAccessor<1>
+ */
+template <>
+class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFSubstructAccessor (Triangulation<2> *tria,
+                           const int         level,
+                           const int         index,
+                           const void       *local_data) :
+                   MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
+                                    // do this here, since this way the
+                                    // CellAccessor has the possibility to know
+                                    // what a face_iterator is. Otherwise
+                                    // it would have to ask the DoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > face_iterator;
+};
+
+
+
+
+
+
+/**
+ * Grant access to the degrees of freedom on a cell. In fact, since all
+ * access to the degrees of freedom has been enabled by the classes
+ * #DoFLineAccessor<1># and #DoFQuadAccessor<2># for the space dimension
+ * one and two, respectively, this class only collects the pieces
+ * together by deriving from the appropriate #DoF*Accessor# and the
+ * right #CellAccessor<dim># and finally adding two functions which give
+ * access to the neighbors and children as #DoFCellAccessor# objects
+ * rather than #CellAccessor# objects (the latter function was inherited
+ * from the #CellAccessor<dim># class).
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFCellAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data) :
+                   MGDoFSubstructAccessor<dim> (tria,level,index,local_data) {};
+
+                                    /**
+                                     * Return the #i#th neighbor as a MGDoF cell
+                                     * iterator. This function is needed since
+                                     * the neighbor function of the base
+                                     * class returns a cell accessor without
+                                     * access to the MGDoF data.
+                                     */
+    TriaIterator<dim,MGDoFCellAccessor<dim> > neighbor (const unsigned int) const;
+
+                                    /**
+                                     * Return the #i#th child as a MGDoF cell
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a cell accessor without
+                                     * access to the DoF data.
+                                     */
+    TriaIterator<dim,MGDoFCellAccessor<dim> > child (const unsigned int) const;
+
+                                    /**
+                                     * Return an iterator to the #i#th face
+                                     * of this cell.
+                                     *
+                                     * This function is not implemented in 1D,
+                                     * and maps to MGDoFQuadAccessor::line in 2D.
+                                     */
+    typename MGDoFSubstructAccessor<dim>::face_iterator
+    face (const unsigned int i) const;
+
+                                    /**
+                                     * Return the value of the given vector
+                                     * restricted to the dofs of this
+                                     * cell in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the multilevel
+                                     * numbering local to the present
+                                     * level of this cell. The vector shall
+                                     * therefore have the same number of
+                                     * entries as there are degrees of
+                                     * freedom on this level.
+                                     */
+    void get_mg_dof_values (const dVector &values,
+                           dVector       &dof_values) const;
+
+
+                                    /**
+                                     *  Exception
+                                     */
+    DeclException0 (ExcNotUsefulForThisDimension);
+};
+
+
+
+
+
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+/* end of #ifndef __mg_dof_accessor_H */
+#endif
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
diff --git a/deal.II/deal.II/include/dofs/mg_dof_handler.h b/deal.II/deal.II/include/dofs/mg_dof_handler.h
new file mode 100644 (file)
index 0000000..21ad59d
--- /dev/null
@@ -0,0 +1,253 @@
+/*----------------------------   mg_dof.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __mg_dof_H
+#define __mg_dof_H
+/*----------------------------   mg_dof.h     ---------------------------*/
+
+
+#include <grid/dof.h>
+
+
+
+// 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 MGDoFHandler : public DoFHandler<dim> {
+  public:
+                                    /**
+                                     * Constructor. Take #tria# as the
+                                     * triangulation to work on.
+                                     */
+    MGDoFHandler (Triangulation<dim> *tria);
+
+                                    /**
+                                     * Destructor
+                                     */
+    virtual ~MGDoFHandler ();
+    
+                                    /**
+                                     * Go through the triangulation and
+                                     * distribute the degrees of freedoms
+                                     * needed for the given finite element
+                                     * according to the given distribution
+                                     * method. We first call the #DoFHandler#'s
+                                     * function and then distribute the
+                                     * levelwise numbers.
+                                     *
+                                     * A copy of the transferred finite
+                                     * element is stored.
+                                     *
+                                     * This function uses the user flags of the
+                                     * triangulation object, so make sure you
+                                     * don't need them after calling this
+                                     * function, or if so store them.
+                                     */
+    virtual void distribute_dofs (const FiniteElementBase<dim> &);
+
+                                    /**
+                                     * Renumber the degrees of freedom according
+                                     * to the given scheme, eventually using
+                                     * constraint information and the given
+                                     * starting points. The starting points
+                                     * default to an empty list, the use of
+                                     * constraint information defaults to
+                                     * false.
+                                     *
+                                     * See the general documentation of the
+                                     * #DoFHandler# class for more details.
+                                     */
+    virtual void renumber_dofs (const RenumberingMethod method,
+                               const bool use_constraints         = false,
+                               const vector<int> &starting_points = vector<int>());
+
+
+  private:
+
+                                    /** We need for each vertex a list of the
+                                     * degree of freedom indices on each of the
+                                     * levels this vertex lives on. Since most
+                                     * vertices live only on a few levels, it
+                                     * is not economical to reserve indices for
+                                     * all the levels there are; rather, we
+                                     * create an object which holds the indices
+                                     * on those levels only where the vertex
+                                     * lives. To construct such an array, it
+                                     * is necessary to know beforehand which
+                                     * is the coarsest level the vertex lives
+                                     * on, how many levels it lives on and
+                                     * how many dofs there are on each vertex.
+                                     * If we have this information, we can
+                                     * allocate exactly the amount of memory
+                                     * which is needed and need not handle
+                                     * growing arrays and the like.
+                                     */
+    class MGVertexDoFs {
+      public:
+       
+                                        /**
+                                         * Constructor. Allocates memory and
+                                         * sets all indices to #-1#.
+                                         */
+       MGVertexDoFs (const unsigned int coarsest_level,
+                     const unsigned int n_levels,
+                     const unsigned int dofs_per_vertex);
+
+                                        /**
+                                         * Destructor
+                                         */
+       ~MGVertexDoFs ();
+
+                                        /**
+                                         * Set the index with number
+                                         * #dof_number# of this vertex on
+                                         * #level# to the given index. To
+                                         * compute the position in the array,
+                                         * one has to specify how many
+                                         * dofs per vertex there are. It is
+                                         * not checked that the level number
+                                         * is below the number of the finest
+                                         * level this vertex lives on.
+                                         *
+                                         * The function is inline, so should
+                                         * be reasonably fast.
+                                         */
+       void set_index (const unsigned int level,
+                       const unsigned int dof_number,
+                       const unsigned int dofs_per_vertex,
+                       const unsigned int index);
+       
+                                        /**
+                                         * Return the index with number
+                                         * #dof_number# of this vertex on
+                                         * #level#. To
+                                         * compute the position in the array,
+                                         * one has to specify how many
+                                         * dofs per vertex there are. It is
+                                         * not checked that the level number
+                                         * is below the number of the finest
+                                         * level this vertex lives on.
+                                         *
+                                         * The function is inline, so should
+                                         * be reasonably fast.
+                                         */
+       int get_index (const unsigned int level,
+                      const unsigned int dof_number,
+                      const unsigned int dofs_per_vertex) const;
+       
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcNoMemory);
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcInvalidIndex);
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException1 (ExcInvalidLevel,
+                       int,
+                       << "The given level number " << arg1 << " is below the "
+                       << "coarsest level this vertex lives on.");
+       
+      private:
+                                        /**
+                                         * Store the coarsest level this
+                                         * vertex lives on. This is used
+                                         * as an offset when accessing the
+                                         * dofs of a specific level.
+                                         */
+       const unsigned int coarsest_level;
+
+                                        /**
+                                         * Array holding the indices.
+                                         */
+       int *const         indices;
+    };
+
+    
+                                    /**
+                                     * Reserve enough space for the MG dof
+                                     * indices for a given triangulation.
+                                     */
+    void reserve_space ();
+
+                                    /**
+                                     * Space to store the DoF numbers for the
+                                     * different levels. Unlike the #levels#
+                                     * object in the #DoFHandler#, these are
+                                     * not global numbers but rather are
+                                     * numbers which start from zero on each
+                                     * level.
+                                     */
+    vector<DoFLevel<dim>*>    mg_levels;
+
+                                    /**
+                                     * For each vertex there is a list of
+                                     * indices of the degrees of freedom indices
+                                     * on the different levels it lives on and
+                                     * which are these levels.
+                                     */
+    vector<MGVertexDoFs>      mg_vertex_dofs;
+    
+                                    /**
+                                     * Vectors storing the number of degrees of
+                                     * freedom on each level.
+                                     */
+    vector<unsigned int>      mg_used_dofs;
+
+    friend class MGDoFLineAccessor<dim, LineAccessor<dim> >;
+    friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
+    friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
+    friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+};
+
+    
+
+
+
+
+
+/* ----------------------- Inline functions of MGVertexDoFs -------------------*/
+
+
+template <int dim>
+inline
+void MGDoFHandler<dim>::MGVertexDoFs::set_index  (const unsigned int level,
+                                                 const unsigned int dof_number,
+                                                 const unsigned int dofs_per_vertex,
+                                                 const unsigned int index) {
+  Assert (level >= coarsest_level, ExcInvalidLevel(level));
+  Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
+  
+  indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index;
+};
+
+
+
+
+template <int dim>
+inline
+int MGDoFHandler<dim>::MGVertexDoFs::get_index  (const unsigned int level,
+                                                const unsigned int dof_number,
+                                                const unsigned int dofs_per_vertex) const {
+  Assert (level >= coarsest_level, ExcInvalidLevel(level));
+  Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
+  
+  return indices[(level-coarsest_level)*dofs_per_vertex + dof_number];
+};
+
+
+
+/*----------------------------   mg_dof.h     ---------------------------*/
+/* end of #ifndef __mg_dof_H */
+#endif
+/*----------------------------   mg_dof.h     ---------------------------*/
diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h
new file mode 100644 (file)
index 0000000..0327475
--- /dev/null
@@ -0,0 +1,505 @@
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __mg_dof_accessor_H
+#define __mg_dof_accessor_H
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+
+
+#include <grid/dof_accessor.h>
+
+
+// forward declaration
+template <int dim> class MGDoFHandler;
+
+
+
+
+/**
+ * Define the basis for accessors to the degrees of freedom for
+ * a multigrid DoF handler object.
+ *
+ * Note that it is allowed to construct an object of which the
+ * #mg_dof_handler# pointer is a Null pointer. Such an object would
+ * result in a strange kind of behaviour, though every reasonable
+ * operating system should disallow access through that pointer.
+ * The reason we do not check for the null pointer in the
+ * constructor which gets passed the #DoFHandler# pointer is that
+ * if we did we could not make dof iterators member of other classes
+ * (like in the #FEValues# class) if we did not know about the
+ * #DoFHandler# object to be used upon construction of that object.
+ * Through the way this class is implemented here, we allow the
+ * creation of a kind of virgin object which only gets useful if
+ * assigned to from another object before first usage.
+ *
+ * Opposite to construction, it is not possible to copy an object
+ * which has an invalid dof handler pointer. This is to guarantee
+ * that every iterator which is once assigned to is a valid
+ * object. However, this assertion only holds in debug mode, when
+ * the #Assert# macro is switched on.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFAccessor {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFAccessor () : mg_dof_handler(0) {
+      Assert (false, ExcInvalidObject());
+    };
+
+                                    /**
+                                     * This should be the default constructor.
+                                     */
+    MGDoFAccessor (MGDoFHandler<dim> *mg_dof_handler) :
+                   mg_dof_handler(mg_dof_handler) {};
+
+                                    /**
+                                     * Reset the DoF handler pointer.
+                                     */
+    void set_mg_dof_handler (MGDoFHandler<dim> *dh) {
+      Assert (dh != 0, ExcInvalidObject());
+      mg_dof_handler = dh;
+    };
+
+                                    /**
+                                     * Copy operator.
+                                     */
+    MGDoFAccessor<dim> & operator = (const MGDoFAccessor<dim> &da) {
+      set_dof_handler (da.mg_dof_handler);
+      return *this;
+    };
+
+  protected:
+                                    /**
+                                     * Store the address of the #MGDoFHandler# object
+                                     * to be accessed.
+                                     */
+    MGDoFHandler<dim> *mg_dof_handler;  
+};
+
+
+
+
+
+/**
+ * Grant access to the degrees of freedom located on lines.
+ * This class follows mainly the route laid out by the accessor library
+ * declared in the triangulation library (\Ref{TriaAccessor}). It enables
+ * the user to access the degrees of freedom on the lines (there are similar
+ * versions for the DoFs on quads, etc), where the dimension of the underlying
+ * triangulation does not really matter (i.e. this accessor works with the
+ * lines in 1D-, 2D-, etc dimensions).
+ *
+ *
+ * \subsection{Usage}
+ *
+ * The \Ref{DoFDimensionInfo} classes inherited by the \Ref{DoFHandler} classes
+ * declare typedefs to iterators using the accessors declared in this class
+ * hierarchy tree. Usage is best to happens through these typedefs, since they
+ * are more secure to changes in the class naming and template interface as well
+ * as they provide easier typing (much less complicated names!).
+ * 
+ * 
+ * \subsection{Notes about the class hierarchy structure}
+ *
+ * The class hierarchy seems to be a bit confused here. The reason for this is
+ * that we would really like to derive a #DoFLineAccessor# from a #LineAccessor#.
+ * Unfortunately, we would run into problems, if we wanted a #DoFLineAccessor#
+ * in one spatial dimension, in which case a line is also a cell. The traditional
+ * solution would be to declare a #DoFCellAccessor<1># which is derived from
+ * #DoFLineAccessor<1># and #CellAccessor<1># (the #DoFLineAccessor<dim># cannot
+ * itself be derived from #CellAccessor<dim># since a line is not a cell
+ * unless in one space dimension), but since a #DoFLineAccessor# and a
+ * #CellAccessor# are both derived from #TriaAccessor#, we would have to make
+ * the last derivation virtual.
+ *
+ * Since we want to avoid virtual inheritance since this involves another
+ * indirection in every member variable access, we chose another way: we
+ * pass a second template parameter to a #DoFLineAccessor# which tells it
+ * which class to be derived from: if we are in one spatial dimension, the
+ * base class is to be #CellAccessor<1>#, in two or more dimensions it
+ * is a #LineAccessor<dim>#, i.e. am accessor to lines without the missing
+ * functionality needed for cells (neighbors, etc.).
+ *
+ * This way we can declare a #DoFCellAccessor# in one dimension by deriving
+ * from #DoFLineAccessor<1,CellAccessor<1> >#, thus getting the cell
+ * functionality through the #DoFLineAccessor# instead of through a virtual
+ * multiple inheritance of #DoFLineAccessor# and #CellAccessor<1>#.
+ *
+ * The same concept is used with #DoFQuadAccessor# classes etc.
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim, typename BaseClass>
+class MGDoFLineAccessor :  public MGDoFAccessor<dim>,
+                          public DoFLineAccessor<dim, BaseClass> {
+  public:
+                                    /**
+                                     * Default constructor, unused thus
+                                     * not implemented.
+                                     */
+    MGDoFLineAccessor ();
+    
+                                    /**
+                                     * Constructor. The #local_data#
+                                     * argument is assumed to be a pointer
+                                     * to a #MGDoFHandler<dim># object.
+                                     */
+    MGDoFLineAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data);
+    
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * of freedom of this line on the level
+                                     * this line lives on.
+                                     */
+    int mg_dof_index (const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * of freedom of this line on the
+                                     * level this line lives on to #index#.
+                                     */
+    void set_mg_dof_index (const unsigned int i, const int index) const;
+
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * on the #vertex#th vertex for the
+                                     * level this line lives on.
+                                     */
+    int mg_vertex_dof_index (const unsigned int vertex,
+                            const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * on the #vertex#th vertex to #index#
+                                     * for the level this line lives on.
+                                     */
+    void set_mg_vertex_dof_index (const unsigned int vertex,
+                                 const unsigned int i,
+                                 const int          index) const;
+
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * line in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, 
+                                     * dofs on line.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the local numbering
+                                     * for the level this line lives on.
+                                     */
+    void get_mg_dof_indices (vector<int> &dof_indices) const;
+
+                                    /**
+                                     * Return the #i#th child as a MGDoF line
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a line accessor without
+                                     * access to the MG data.
+                                     */
+    TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > child (const unsigned int) const;
+    
+                                    /**
+                                     * Implement the copy operator needed
+                                     * for the iterator classes.
+                                     */
+    void copy_from (const MGDoFLineAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+
+/**
+ * Grant access to the multilevel degrees of freedom located on quads.
+ *
+ * @see DoFLineAccessor
+ */
+template <int dim, typename BaseClass>
+class MGDoFQuadAccessor :  public MGDoFAccessor<dim>,
+                          public DoFQuadAccessor<dim, BaseClass> {
+  public:
+                                    /**
+                                     * Default constructor, unused thus
+                                     * not implemented.
+                                     */
+    MGDoFQuadAccessor ();
+
+                                    /**
+                                     * Constructor. The #local_data#
+                                     * argument is assumed to be a pointer
+                                     * to a #DoFHandler<dim># object.
+                                     */
+    MGDoFQuadAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data);
+    
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * of freedom of this quad on the level
+                                     * this quad lives on.
+                                     */
+    int mg_dof_index (const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * of freedom of this quad on the
+                                     * level this quad lives on to #index#.
+                                     */
+    void set_mg_dof_index (const unsigned int i, const int index) const;
+
+                                    /**
+                                     * Return the index of the #i#th degree
+                                     * on the #vertex#th vertex for the level
+                                     * this quad lives on.
+                                     */
+    int mg_vertex_dof_index (const unsigned int vertex,
+                            const unsigned int i) const;
+
+                                    /**
+                                     * Set the index of the #i#th degree
+                                     * on the #vertex#th vertex for the
+                                     * level this quad lives on to #index#.
+                                     */
+    void set_mg_vertex_dof_index (const unsigned int vertex,
+                                 const unsigned int i,
+                                 const int          index) const;
+
+                                    /**
+                                     * Return the indices of the dofs of this
+                                     * quad in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the local numbering
+                                     * for the level this quad lives on.
+                                     */
+    void get_mg_dof_indices (vector<int> &dof_indices) const;
+
+                                    /**
+                                     * Return a pointer to the #i#th line
+                                     * bounding this #Quad#.
+                                     */
+    TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+    line (const unsigned int i) const;
+
+                                    /**
+                                     * Return the #i#th child as a DoF quad
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a quad accessor without
+                                     * access to the DoF data.
+                                     */
+    TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > child (const unsigned int) const;
+    
+                                    /**
+                                     * Implement the copy operator needed
+                                     * for the iterator classes.
+                                     */
+    void copy_from (const DoFQuadAccessor<dim,BaseClass> &a);
+};
+
+
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * Rationale for the declaration of members for this class: gcc 2.8 has a bug
+ * when deriving from explicitely specialized classes which materializes in
+ * the calculation of wrong addresses of member variables. By declaring the
+ * general template of #DoFSubstructAccessor# to have the same object layout as
+ * the specialized versions (using the same base classes), we fool the compiler,
+ * which still looks in the wrong place for the addresses but finds the
+ * right information. This way, at least ot works.
+ *
+ * Insert a guard, however, in the constructor to avoid that anyone (including
+ * the compiler) happens to use this class.
+ */
+template <int dim>
+class MGDoFSubstructAccessor : public MGDoFAccessor<dim>,
+                              public TriaAccessor<dim> {
+  public:
+    DoFSubstructAccessor () {
+      Assert (false, ExcInternalError());
+    };
+
+    DeclException0 (ExcInternalError);
+};
+
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * \subsection{Rationale}
+ *
+ * This class is only a wrapper class used to do kind of a typedef
+ * 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> >;
+ * \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
+ * the \Ref{DoFCellAccessor} dimension independent as an inheritance from
+ * #DoFSubstructAccessor<dim>#. If we had not declared these
+ * types, we would have to write two class declarations, one for
+ * #DoFCellAccessor<1>#, derived from #DoFLineAccessor<1,CellAccessor<1> >#
+ * and one for #DoFCellAccessor<2>#, derived from
+ * #DoFQuadAccessor<2,CellAccessor<2> >#.
+ */
+template <>
+class MGDoFSubstructAccessor<1> :  public MGDoFLineAccessor<1,CellAccessor<1> > {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFSubstructAccessor (Triangulation<1> *tria,
+                           const int         level,
+                           const int         index,
+                           const void       *local_data) :
+                   MGDoFLineAccessor<1,CellAccessor<1> > (tria,level,index,local_data) {};
+                                    // do this here, since this way the
+                                    // CellAccessor has the possibility to know
+                                    // what a face_iterator is. Otherwise
+                                    // it would have to ask the MGDoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef void * face_iterator;
+};
+
+
+
+/**
+ * Intermediate, "typedef"-class, not for public use.
+ *
+ * @see MGDoFSubstructAccessor<1>
+ */
+template <>
+class MGDoFSubstructAccessor<2> : public MGDoFQuadAccessor<2,CellAccessor<2> > {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFSubstructAccessor (Triangulation<2> *tria,
+                           const int         level,
+                           const int         index,
+                           const void       *local_data) :
+                   MGDoFQuadAccessor<2,CellAccessor<2> > (tria,level,index,local_data) {};
+                                    // do this here, since this way the
+                                    // CellAccessor has the possibility to know
+                                    // what a face_iterator is. Otherwise
+                                    // it would have to ask the DoFHandler<dim>
+                                    // which would need another big include
+                                    // file and maybe cyclic interdependence
+    typedef TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > > face_iterator;
+};
+
+
+
+
+
+
+/**
+ * Grant access to the degrees of freedom on a cell. In fact, since all
+ * access to the degrees of freedom has been enabled by the classes
+ * #DoFLineAccessor<1># and #DoFQuadAccessor<2># for the space dimension
+ * one and two, respectively, this class only collects the pieces
+ * together by deriving from the appropriate #DoF*Accessor# and the
+ * right #CellAccessor<dim># and finally adding two functions which give
+ * access to the neighbors and children as #DoFCellAccessor# objects
+ * rather than #CellAccessor# objects (the latter function was inherited
+ * from the #CellAccessor<dim># class).
+ *
+ * @author Wolfgang Bangerth, 1998
+ */
+template <int dim>
+class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
+  public:
+                                    /**
+                                     * Constructor
+                                     */
+    MGDoFCellAccessor (Triangulation<dim> *tria,
+                      const int           level,
+                      const int           index,
+                      const void         *local_data) :
+                   MGDoFSubstructAccessor<dim> (tria,level,index,local_data) {};
+
+                                    /**
+                                     * Return the #i#th neighbor as a MGDoF cell
+                                     * iterator. This function is needed since
+                                     * the neighbor function of the base
+                                     * class returns a cell accessor without
+                                     * access to the MGDoF data.
+                                     */
+    TriaIterator<dim,MGDoFCellAccessor<dim> > neighbor (const unsigned int) const;
+
+                                    /**
+                                     * Return the #i#th child as a MGDoF cell
+                                     * iterator. This function is needed since
+                                     * the child function of the base
+                                     * class returns a cell accessor without
+                                     * access to the DoF data.
+                                     */
+    TriaIterator<dim,MGDoFCellAccessor<dim> > child (const unsigned int) const;
+
+                                    /**
+                                     * Return an iterator to the #i#th face
+                                     * of this cell.
+                                     *
+                                     * This function is not implemented in 1D,
+                                     * and maps to MGDoFQuadAccessor::line in 2D.
+                                     */
+    typename MGDoFSubstructAccessor<dim>::face_iterator
+    face (const unsigned int i) const;
+
+                                    /**
+                                     * Return the value of the given vector
+                                     * restricted to the dofs of this
+                                     * cell in the standard ordering: dofs
+                                     * on vertex 0, dofs on vertex 1, etc,
+                                     * dofs on line 0, dofs on line 1, etc,
+                                     * dofs on quad 0, etc.
+                                     *
+                                     * It is assumed that the vector already
+                                     * has the right size beforehand. The
+                                     * indices refer to the multilevel
+                                     * numbering local to the present
+                                     * level of this cell. The vector shall
+                                     * therefore have the same number of
+                                     * entries as there are degrees of
+                                     * freedom on this level.
+                                     */
+    void get_mg_dof_values (const dVector &values,
+                           dVector       &dof_values) const;
+
+
+                                    /**
+                                     *  Exception
+                                     */
+    DeclException0 (ExcNotUsefulForThisDimension);
+};
+
+
+
+
+
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
+/* end of #ifndef __mg_dof_accessor_H */
+#endif
+/*----------------------------   mg_dof_accessor.h     ---------------------------*/
diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h
new file mode 100644 (file)
index 0000000..21ad59d
--- /dev/null
@@ -0,0 +1,253 @@
+/*----------------------------   mg_dof.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __mg_dof_H
+#define __mg_dof_H
+/*----------------------------   mg_dof.h     ---------------------------*/
+
+
+#include <grid/dof.h>
+
+
+
+// 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 MGDoFHandler : public DoFHandler<dim> {
+  public:
+                                    /**
+                                     * Constructor. Take #tria# as the
+                                     * triangulation to work on.
+                                     */
+    MGDoFHandler (Triangulation<dim> *tria);
+
+                                    /**
+                                     * Destructor
+                                     */
+    virtual ~MGDoFHandler ();
+    
+                                    /**
+                                     * Go through the triangulation and
+                                     * distribute the degrees of freedoms
+                                     * needed for the given finite element
+                                     * according to the given distribution
+                                     * method. We first call the #DoFHandler#'s
+                                     * function and then distribute the
+                                     * levelwise numbers.
+                                     *
+                                     * A copy of the transferred finite
+                                     * element is stored.
+                                     *
+                                     * This function uses the user flags of the
+                                     * triangulation object, so make sure you
+                                     * don't need them after calling this
+                                     * function, or if so store them.
+                                     */
+    virtual void distribute_dofs (const FiniteElementBase<dim> &);
+
+                                    /**
+                                     * Renumber the degrees of freedom according
+                                     * to the given scheme, eventually using
+                                     * constraint information and the given
+                                     * starting points. The starting points
+                                     * default to an empty list, the use of
+                                     * constraint information defaults to
+                                     * false.
+                                     *
+                                     * See the general documentation of the
+                                     * #DoFHandler# class for more details.
+                                     */
+    virtual void renumber_dofs (const RenumberingMethod method,
+                               const bool use_constraints         = false,
+                               const vector<int> &starting_points = vector<int>());
+
+
+  private:
+
+                                    /** We need for each vertex a list of the
+                                     * degree of freedom indices on each of the
+                                     * levels this vertex lives on. Since most
+                                     * vertices live only on a few levels, it
+                                     * is not economical to reserve indices for
+                                     * all the levels there are; rather, we
+                                     * create an object which holds the indices
+                                     * on those levels only where the vertex
+                                     * lives. To construct such an array, it
+                                     * is necessary to know beforehand which
+                                     * is the coarsest level the vertex lives
+                                     * on, how many levels it lives on and
+                                     * how many dofs there are on each vertex.
+                                     * If we have this information, we can
+                                     * allocate exactly the amount of memory
+                                     * which is needed and need not handle
+                                     * growing arrays and the like.
+                                     */
+    class MGVertexDoFs {
+      public:
+       
+                                        /**
+                                         * Constructor. Allocates memory and
+                                         * sets all indices to #-1#.
+                                         */
+       MGVertexDoFs (const unsigned int coarsest_level,
+                     const unsigned int n_levels,
+                     const unsigned int dofs_per_vertex);
+
+                                        /**
+                                         * Destructor
+                                         */
+       ~MGVertexDoFs ();
+
+                                        /**
+                                         * Set the index with number
+                                         * #dof_number# of this vertex on
+                                         * #level# to the given index. To
+                                         * compute the position in the array,
+                                         * one has to specify how many
+                                         * dofs per vertex there are. It is
+                                         * not checked that the level number
+                                         * is below the number of the finest
+                                         * level this vertex lives on.
+                                         *
+                                         * The function is inline, so should
+                                         * be reasonably fast.
+                                         */
+       void set_index (const unsigned int level,
+                       const unsigned int dof_number,
+                       const unsigned int dofs_per_vertex,
+                       const unsigned int index);
+       
+                                        /**
+                                         * Return the index with number
+                                         * #dof_number# of this vertex on
+                                         * #level#. To
+                                         * compute the position in the array,
+                                         * one has to specify how many
+                                         * dofs per vertex there are. It is
+                                         * not checked that the level number
+                                         * is below the number of the finest
+                                         * level this vertex lives on.
+                                         *
+                                         * The function is inline, so should
+                                         * be reasonably fast.
+                                         */
+       int get_index (const unsigned int level,
+                      const unsigned int dof_number,
+                      const unsigned int dofs_per_vertex) const;
+       
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcNoMemory);
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException0 (ExcInvalidIndex);
+                                        /**
+                                         * Exception.
+                                         */
+       DeclException1 (ExcInvalidLevel,
+                       int,
+                       << "The given level number " << arg1 << " is below the "
+                       << "coarsest level this vertex lives on.");
+       
+      private:
+                                        /**
+                                         * Store the coarsest level this
+                                         * vertex lives on. This is used
+                                         * as an offset when accessing the
+                                         * dofs of a specific level.
+                                         */
+       const unsigned int coarsest_level;
+
+                                        /**
+                                         * Array holding the indices.
+                                         */
+       int *const         indices;
+    };
+
+    
+                                    /**
+                                     * Reserve enough space for the MG dof
+                                     * indices for a given triangulation.
+                                     */
+    void reserve_space ();
+
+                                    /**
+                                     * Space to store the DoF numbers for the
+                                     * different levels. Unlike the #levels#
+                                     * object in the #DoFHandler#, these are
+                                     * not global numbers but rather are
+                                     * numbers which start from zero on each
+                                     * level.
+                                     */
+    vector<DoFLevel<dim>*>    mg_levels;
+
+                                    /**
+                                     * For each vertex there is a list of
+                                     * indices of the degrees of freedom indices
+                                     * on the different levels it lives on and
+                                     * which are these levels.
+                                     */
+    vector<MGVertexDoFs>      mg_vertex_dofs;
+    
+                                    /**
+                                     * Vectors storing the number of degrees of
+                                     * freedom on each level.
+                                     */
+    vector<unsigned int>      mg_used_dofs;
+
+    friend class MGDoFLineAccessor<dim, LineAccessor<dim> >;
+    friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
+    friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
+    friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+};
+
+    
+
+
+
+
+
+/* ----------------------- Inline functions of MGVertexDoFs -------------------*/
+
+
+template <int dim>
+inline
+void MGDoFHandler<dim>::MGVertexDoFs::set_index  (const unsigned int level,
+                                                 const unsigned int dof_number,
+                                                 const unsigned int dofs_per_vertex,
+                                                 const unsigned int index) {
+  Assert (level >= coarsest_level, ExcInvalidLevel(level));
+  Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
+  
+  indices[(level-coarsest_level)*dofs_per_vertex + dof_number] = index;
+};
+
+
+
+
+template <int dim>
+inline
+int MGDoFHandler<dim>::MGVertexDoFs::get_index  (const unsigned int level,
+                                                const unsigned int dof_number,
+                                                const unsigned int dofs_per_vertex) const {
+  Assert (level >= coarsest_level, ExcInvalidLevel(level));
+  Assert (dof_number < dofs_per_vertex, ExcInvalidIndex ());
+  
+  return indices[(level-coarsest_level)*dofs_per_vertex + dof_number];
+};
+
+
+
+/*----------------------------   mg_dof.h     ---------------------------*/
+/* end of #ifndef __mg_dof_H */
+#endif
+/*----------------------------   mg_dof.h     ---------------------------*/
diff --git a/deal.II/deal.II/source/dofs/mg_dof_accessor.cc b/deal.II/deal.II/source/dofs/mg_dof_accessor.cc
new file mode 100644 (file)
index 0000000..322b395
--- /dev/null
@@ -0,0 +1,426 @@
+/* $Id$ */
+
+
+#include <grid/mg_dof_accessor.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
+#include <fe/fe.h>
+
+#include <lac/dvector.h>
+#include <lac/dfmatrix.h>
+#include <lac/dsmatrix.h>
+
+
+
+
+
+
+/* ------------------------ MGDoFLineAccessor --------------------------- */
+
+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))) {};
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::mg_dof_index (const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_line,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i];
+};
+
+
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::set_mg_dof_index (const unsigned int i,
+                                                        const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_line,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i] = index;
+};
+
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::mg_vertex_dof_index (const unsigned int vertex,
+                                                          const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<2, ExcInvalidIndex (i,0,2));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::set_mg_vertex_dof_index (const unsigned int vertex,
+                                                               const unsigned int i,
+                                                               const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<2, ExcInvalidIndex (i,0,2));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (2*dof_handler->get_selected_fe().dofs_per_vertex +
+                                dof_handler->get_selected_fe().dofs_per_line),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int d=0; d<dofs_per_line; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> >
+MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > q (tria,
+                                                        present_level+1,
+                                                        child_index (i),
+                                                        mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+
+/* ------------------------ MGDoFQuadAccessor --------------------------- */
+
+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))) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::mg_dof_index (const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_quad,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::set_mg_dof_index (const unsigned int i,
+                                                        const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_quad,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::mg_vertex_dof_index (const unsigned int vertex,
+                                                          const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+  
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::set_mg_vertex_dof_index (const unsigned int vertex,
+                                                               const unsigned int i,
+                                                               const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (4*dof_handler->get_selected_fe().dofs_per_vertex +
+                                4*dof_handler->get_selected_fe().dofs_per_line +
+                                dof_handler->get_selected_fe().dofs_per_quad),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_selected_fe().dofs_per_quad;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->mg_dof_index(d);
+  for (unsigned int d=0; d<dofs_per_quad; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
+  Assert (i<4, ExcInvalidIndex (i, 0, 4));
+
+  return TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+    (
+      tria,
+      present_level,
+      line_index (i),
+      mg_dof_handler
+    );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> >
+MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > q (tria,
+                                                        present_level+1,
+                                                        child_index (i),
+                                                        mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+
+
+/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::neighbor (const unsigned int i) const {
+  TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+                                              neighbor_level (i),
+                                              neighbor_index (i),
+                                              mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
+#endif
+  return q;
+};
+
+
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+                                              present_level+1,
+                                              child_index (i),
+                                              mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+MGDoFSubstructAccessor<1>::face_iterator
+MGDoFCellAccessor<1>::face (const unsigned int) const {
+  Assert (false, ExcNotUsefulForThisDimension());
+  return 0;
+};
+
+
+
+template <>
+void
+MGDoFCellAccessor<1>::get_mg_dof_values (const dVector &values,
+                                        dVector       &dof_values) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs,
+         ExcVectorDoesNotMatch());
+  Assert (values.size() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line;
+  vector<double>::iterator next_dof_value=dof_values.begin();
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+  for (unsigned int d=0; d<dofs_per_line; ++d)
+    *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+MGDoFSubstructAccessor<2>::face_iterator
+MGDoFCellAccessor<2>::face (const unsigned int i) const {
+  return line(i);
+};
+
+
+
+template <>
+void
+MGDoFCellAccessor<2>::get_mg_dof_values (const dVector &values,
+                                        dVector       &dof_values) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs,
+         ExcVectorDoesNotMatch());
+  Assert (values.size() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_selected_fe().dofs_per_quad;
+  vector<double>::iterator next_dof_value=dof_values.begin();
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
+  for (unsigned int d=0; d<dofs_per_quad; ++d)
+    *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+
+
+// explicit instantiations
+#if deal_II_dimension == 1
+template class MGDoFLineAccessor<1,CellAccessor<1> >;
+template class MGDoFCellAccessor<1>;
+
+template class TriaRawIterator<1,MGDoFCellAccessor<1> >;
+template class TriaIterator<1,MGDoFCellAccessor<1> >;
+template class TriaActiveIterator<1,MGDoFCellAccessor<1> >;
+#endif
+
+#if deal_II_dimension == 2
+template class MGDoFLineAccessor<2,LineAccessor<2> >;
+template class MGDoFQuadAccessor<2,QuadAccessor<2> >;
+template class MGDoFQuadAccessor<2,CellAccessor<2> >;
+template class MGDoFCellAccessor<2>;
+
+template class TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaRawIterator<2,MGDoFCellAccessor<2> >;
+template class TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaIterator<2,MGDoFCellAccessor<2> >;
+template class TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaActiveIterator<2,MGDoFCellAccessor<2> >;
+#endif
+
+
diff --git a/deal.II/deal.II/source/dofs/mg_dof_handler.cc b/deal.II/deal.II/source/dofs/mg_dof_handler.cc
new file mode 100644 (file)
index 0000000..14b9a55
--- /dev/null
@@ -0,0 +1,92 @@
+/* $Id$ */
+
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/dof_accessor.h>
+#include <grid/dof_constraints.h>
+#include <grid/tria_levels.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.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])
+{
+  Assert (indices != 0, ExcNoMemory ());
+
+  for (unsigned int i=0; i<n_levels*dofs_per_vertex; ++i)
+    indices[i] = -1;
+};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+  delete[] indices;
+};
+
+
+
+
+
+
+
+/* ------------------------ MGDoFHandler ------------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::MGDoFHandler (Triangulation<dim> *tria) :
+               DoFHandler<dim> (tria)
+{};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::~MGDoFHandler () {};
+
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+                                  // first distribute global dofs
+  DoFHandler<dim>::distribute_dofs (fe);
+
+
+                                  // reserve space for the MG dof numbers
+  reserve_space ();
+  
+};
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void MGDoFHandler<1>::reserve_space () {
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void MGDoFHandler<2>::reserve_space () {
+};
+
+#endif
+
+
+// explicite instantiations
+template class MGDoFHandler<deal_II_dimension>;
diff --git a/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc b/deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
new file mode 100644 (file)
index 0000000..322b395
--- /dev/null
@@ -0,0 +1,426 @@
+/* $Id$ */
+
+
+#include <grid/mg_dof_accessor.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_iterator.templates.h>
+#include <fe/fe.h>
+
+#include <lac/dvector.h>
+#include <lac/dfmatrix.h>
+#include <lac/dsmatrix.h>
+
+
+
+
+
+
+/* ------------------------ MGDoFLineAccessor --------------------------- */
+
+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))) {};
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::mg_dof_index (const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_line,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i];
+};
+
+
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::set_mg_dof_index (const unsigned int i,
+                                                        const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_line,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_line));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->line_dofs[present_index*dof_handler->get_selected_fe().dofs_per_line+i] = index;
+};
+
+
+
+
+template <int dim, typename BaseClass>
+int MGDoFLineAccessor<dim,BaseClass>::mg_vertex_dof_index (const unsigned int vertex,
+                                                          const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<2, ExcInvalidIndex (i,0,2));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFLineAccessor<dim,BaseClass>::set_mg_vertex_dof_index (const unsigned int vertex,
+                                                               const unsigned int i,
+                                                               const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<2, ExcInvalidIndex (i,0,2));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFLineAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (2*dof_handler->get_selected_fe().dofs_per_vertex +
+                                dof_handler->get_selected_fe().dofs_per_line),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int d=0; d<dofs_per_line; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> >
+MGDoFLineAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFLineAccessor<dim,BaseClass> > q (tria,
+                                                        present_level+1,
+                                                        child_index (i),
+                                                        mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+
+/* ------------------------ MGDoFQuadAccessor --------------------------- */
+
+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))) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::mg_dof_index (const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_quad,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::set_mg_dof_index (const unsigned int i,
+                                                        const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+                                  // make sure a FE has been selected
+                                  // and enough room was reserved
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_selected_fe().dofs_per_quad,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_quad));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->quad_dofs[present_index*dof_handler->get_selected_fe().dofs_per_quad+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFQuadAccessor<dim,BaseClass>::mg_vertex_dof_index (const unsigned int vertex,
+                                                          const unsigned int i) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+  
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFQuadAccessor<dim,BaseClass>::set_mg_vertex_dof_index (const unsigned int vertex,
+                                                               const unsigned int i,
+                                                               const int index) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_selected_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_selected_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_selected_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFQuadAccessor<dim,BaseClass>::get_mg_dof_indices (vector<int> &dof_indices) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (4*dof_handler->get_selected_fe().dofs_per_vertex +
+                                4*dof_handler->get_selected_fe().dofs_per_line +
+                                dof_handler->get_selected_fe().dofs_per_quad),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_selected_fe().dofs_per_quad;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->mg_dof_index(d);
+  for (unsigned int d=0; d<dofs_per_quad; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFQuadAccessor<dim,BaseClass>::line (const unsigned int i) const {
+  Assert (i<4, ExcInvalidIndex (i, 0, 4));
+
+  return TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+    (
+      tria,
+      present_level,
+      line_index (i),
+      mg_dof_handler
+    );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> >
+MGDoFQuadAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFQuadAccessor<dim,BaseClass> > q (tria,
+                                                        present_level+1,
+                                                        child_index (i),
+                                                        mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+
+
+/*------------------------- Functions: MGDoFCellAccessor -----------------------*/
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::neighbor (const unsigned int i) const {
+  TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+                                              neighbor_level (i),
+                                              neighbor_index (i),
+                                              mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsNeighbor());
+#endif
+  return q;
+};
+
+
+
+template <int dim>
+TriaIterator<dim,MGDoFCellAccessor<dim> >
+MGDoFCellAccessor<dim>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFCellAccessor<dim> > q (tria,
+                                              present_level+1,
+                                              child_index (i),
+                                              mg_dof_handler);
+  
+#ifdef DEBUG
+  if (q.state() != past_the_end)
+    Assert (q->used(), typename TriaAccessor<dim>::ExcUnusedCellAsChild());
+#endif
+  return q;
+};
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+MGDoFSubstructAccessor<1>::face_iterator
+MGDoFCellAccessor<1>::face (const unsigned int) const {
+  Assert (false, ExcNotUsefulForThisDimension());
+  return 0;
+};
+
+
+
+template <>
+void
+MGDoFCellAccessor<1>::get_mg_dof_values (const dVector &values,
+                                        dVector       &dof_values) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs,
+         ExcVectorDoesNotMatch());
+  Assert (values.size() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line;
+  vector<double>::iterator next_dof_value=dof_values.begin();
+  for (unsigned int vertex=0; vertex<2; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+  for (unsigned int d=0; d<dofs_per_line; ++d)
+    *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 2
+
+template <>
+MGDoFSubstructAccessor<2>::face_iterator
+MGDoFCellAccessor<2>::face (const unsigned int i) const {
+  return line(i);
+};
+
+
+
+template <>
+void
+MGDoFCellAccessor<2>::get_mg_dof_values (const dVector &values,
+                                        dVector       &dof_values) const {
+  Assert (dof_handler != 0, ExcInvalidObject());
+  Assert (mg_dof_handler != 0, ExcInvalidObject());
+  Assert (&dof_handler->get_selected_fe() != 0, ExcInvalidObject());
+  Assert (dof_values.size() == dof_handler->get_selected_fe().total_dofs,
+         ExcVectorDoesNotMatch());
+  Assert (values.size() == dof_handler->n_dofs(),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_selected_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_selected_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_selected_fe().dofs_per_quad;
+  vector<double>::iterator next_dof_value=dof_values.begin();
+  for (unsigned int vertex=0; vertex<4; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next_dof_value++ = values(mg_vertex_dof_index(vertex,d));
+  for (unsigned int line=0; line<4; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next_dof_value++ = values(this->line(line)->mg_dof_index(d));
+  for (unsigned int d=0; d<dofs_per_quad; ++d)
+    *next_dof_value++ = values(mg_dof_index(d));
+};
+
+#endif
+
+
+
+
+
+// explicit instantiations
+#if deal_II_dimension == 1
+template class MGDoFLineAccessor<1,CellAccessor<1> >;
+template class MGDoFCellAccessor<1>;
+
+template class TriaRawIterator<1,MGDoFCellAccessor<1> >;
+template class TriaIterator<1,MGDoFCellAccessor<1> >;
+template class TriaActiveIterator<1,MGDoFCellAccessor<1> >;
+#endif
+
+#if deal_II_dimension == 2
+template class MGDoFLineAccessor<2,LineAccessor<2> >;
+template class MGDoFQuadAccessor<2,QuadAccessor<2> >;
+template class MGDoFQuadAccessor<2,CellAccessor<2> >;
+template class MGDoFCellAccessor<2>;
+
+template class TriaRawIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaRawIterator<2,MGDoFCellAccessor<2> >;
+template class TriaIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaIterator<2,MGDoFCellAccessor<2> >;
+template class TriaActiveIterator<2,MGDoFLineAccessor<2,LineAccessor<2> > >;
+template class TriaActiveIterator<2,MGDoFCellAccessor<2> >;
+#endif
+
+
diff --git a/deal.II/deal.II/source/multigrid/mg_dof_handler.cc b/deal.II/deal.II/source/multigrid/mg_dof_handler.cc
new file mode 100644 (file)
index 0000000..14b9a55
--- /dev/null
@@ -0,0 +1,92 @@
+/* $Id$ */
+
+#include <grid/mg_dof.h>
+#include <grid/dof_levels.h>
+#include <grid/dof_accessor.h>
+#include <grid/dof_constraints.h>
+#include <grid/tria_levels.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria.h>
+#include <grid/geometry_info.h>
+#include <fe/fe.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])
+{
+  Assert (indices != 0, ExcNoMemory ());
+
+  for (unsigned int i=0; i<n_levels*dofs_per_vertex; ++i)
+    indices[i] = -1;
+};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::MGVertexDoFs::~MGVertexDoFs () {
+  delete[] indices;
+};
+
+
+
+
+
+
+
+/* ------------------------ MGDoFHandler ------------------------------------- */
+
+template <int dim>
+MGDoFHandler<dim>::MGDoFHandler (Triangulation<dim> *tria) :
+               DoFHandler<dim> (tria)
+{};
+
+
+
+template <int dim>
+MGDoFHandler<dim>::~MGDoFHandler () {};
+
+
+
+
+template <int dim>
+void MGDoFHandler<dim>::distribute_dofs (const FiniteElementBase<dim> &fe) {
+                                  // first distribute global dofs
+  DoFHandler<dim>::distribute_dofs (fe);
+
+
+                                  // reserve space for the MG dof numbers
+  reserve_space ();
+  
+};
+
+
+
+
+#if deal_II_dimension == 1
+
+template <>
+void MGDoFHandler<1>::reserve_space () {
+};
+
+#endif
+
+
+#if deal_II_dimension == 2
+
+template <>
+void MGDoFHandler<2>::reserve_space () {
+};
+
+#endif
+
+
+// explicite instantiations
+template class MGDoFHandler<deal_II_dimension>;

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.