]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Separate the dof renumbering code from the dof handlers into a new class. Add several...
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Apr 1999 11:10:33 +0000 (11:10 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 20 Apr 1999 11:10:33 +0000 (11:10 +0000)
git-svn-id: https://svn.dealii.org/trunk@1190 0785d39b-7218-0410-832d-ea1e28bc413d

16 files changed:
deal.II/deal.II/Attic/examples/dof/dof_test.cc
deal.II/deal.II/include/dofs/dof_handler.h
deal.II/deal.II/include/dofs/mg_dof_accessor.h
deal.II/deal.II/include/dofs/mg_dof_handler.h
deal.II/deal.II/include/grid/tria_accessor.templates.h
deal.II/deal.II/include/multigrid/mg_dof_accessor.h
deal.II/deal.II/include/multigrid/mg_dof_handler.h
deal.II/deal.II/include/numerics/dof_renumbering.h [new file with mode: 0644]
deal.II/deal.II/source/dofs/dof_handler.cc
deal.II/deal.II/source/dofs/dof_renumbering.cc [new file with mode: 0644]
deal.II/deal.II/source/dofs/mg_dof_accessor.cc
deal.II/deal.II/source/dofs/mg_dof_handler.cc
deal.II/deal.II/source/multigrid/mg_dof_accessor.cc
deal.II/deal.II/source/multigrid/mg_dof_handler.cc
tests/big-tests/dof/dof_test.cc
tests/deal.II/dof_test.cc

index cc8a9f83f524700cd13f3d92c0ce50030fccc51b..8ca5178abfabad62e4f98de2020bfe16b2a5e373 100644 (file)
@@ -12,7 +12,7 @@
 #include <lac/sparsematrix.h>
 #include <base/parameter_handler.h>
 #include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
 
 #include <fstream>
 #include <cmath>
@@ -354,7 +354,7 @@ void TestCases<dim>::run (ParameterHandler &prm) {
   dof->distribute_dofs (fe);
 
   cout << "    Renumbering degrees of freedom..." << endl;
-  dof->renumber_dofs (Cuthill_McKee, false);
+  DoFRenumbering::renumber_Cuthill_McKee (*dof);
     
   SparseMatrixStruct sparsity (dof->n_dofs(),
                               dof->max_couplings_between_dofs());
index b32a35c04ac8f8a617f9a2a08e99bd57a6643d95..6252b03017947c94cab1d671cd11590ec35827c0 100644 (file)
@@ -129,29 +129,6 @@ class DoFDimensionInfo<3> {
 
 
 
-/**
- * Give names to the different possibilities of renumbering the degrees
- * of freedom.
- *
- * \begin{itemize}
- * \item #Cuthill_McKee# and #reverse_Cuthill_McKee# traverse the triangulation
- *   in a diagonal, advancing front like method and produce matrices with an
- *   almost minimal bandwidth.
- * \item #reverse_Cuthill_McKey# does the same thing, but numbers the dofs in
- *   the reverse order.
- * \end{itemize}
- *
- * For a description of the algorithms see the book of Schwarz (H.R.Scharz:
- * Methode der finiten Elemente).
- */
-enum RenumberingMethod {
-      Cuthill_McKee,
-      reverse_Cuthill_McKee
-};
-
-
-
-
 
 /**
  * Manage the distribution and numbering of the degrees of freedom for
@@ -203,117 +180,14 @@ enum RenumberingMethod {
  *
  * This numbering implies very large bandwiths of the resulting matrices and
  * is thus vastly suboptimal for some solution algorithms. For this reason,
- * the #DoFHandler# class offers the function #renumber_dofs# which reorders
- * the dof numbering according to some scheme. Presently available are the
- * Cuthill-McKey (CM) and the Reverse Cuthill-McKey algorithm. These algorithms
- * have one major drawback: they require a good starting point, i.e. the degree
- * of freedom index afterwards to be numbered zero. This can thus be given by
- * the user, e.g. by exploiting knowledge of the actual topology of the
- * domain. It is also possible to give several starting indices, which may
- * be used to simulate a simple upstream numbering (by giving the inflow
- * dofs as starting values) or to make preconditioning faster (by letting
- * the dirichlet boundary indices be starting points).
- *
- * If no starting index is given, one is chosen by the program, namely one
- * with the smallest coordination number (the coordination number is the
- * number of other dofs this dof couples with). This dof is usually located
- * on the boundary of the domain. There is, however, large ambiguity in this
- * when using the hierarchical meshes used in this library, since in most
- * cases the computational domain is not approximated by tilting and deforming
- * elements and by plugging together variable numbers of elements at vertices,
- * but rather by hierarchical refinement. There is therefore a large number
- * of dofs with equal coordination numbers. The renumbering algorithms will
- * therefore not give optimal results.
- *
- * In the book of Schwarz (H.R.Schwarz: Methode der finiten Elemente), it is
- * advised to test many starting points, if possible all with the smallest
- * coordination number and also those with slightly higher numbers. However,
- * this seems only possible for meshes with at most several dozen or a few
- * hundred elements found in small engineering problems of the early 1980s
- * (the second edition was published in 1984), but certainly not with those
- * used in this library, featuring several 10,000 to a few 100,000 elements.
- *
- * On the other hand, the need to reduce the bandwidth has decreased since
- * with the mentioned number of cells, only iterative solution methods are
- * able to solve the resulting matrix systems. These, however, are not so
- * demanding with respect to the bandwidth as direct solvers used for
- * smaller problems. Things like upstream numbering become much more important
- * in recent times, so the suboptimality of the renumbering algorithms is
- * not that important any more.
- *
- * 
- * \subsection{Implementation of renumbering schemes}
- *
- * The renumbering algorithms need quite a lot of memory, since they have
- * to store for each dof with which other dofs it couples. This is done
- * using a #SparseMatrixStruct# object used to store the sparsity pattern of
- * matrices. It
- * is not useful for the user to do anything between distributing the dofs
- * and renumbering, i.e. the calls to #DoFHandler::distribute_dofs# and
- * #DoFHandler::renumber_dofs# should follow each other immediately. If
- * you try to create a sparsity pattern or anything else in between, these
- * will be invalid afterwards.
- *
- * The renumbering may take care of dof-to-dof couplings only induced by
- * eliminating constraints. In addition to the memory consumption mentioned
- * above, this also takes quite some computational time, but it may be
- * switched off upon calling the #renumber_dofs# function. This will then
- * give inferior results, since knots in the graph (representing dofs)
- * are not found to be neighbors even if they would be after condensation.
- * 
- * The renumbering algorithms work on a purely algebraic basis, due to the
- * isomorphism between the graph theoretical groundwork underlying the
- * algorithms and binary matrices (matrices of which the entries are binary
- * values) represented by the sparsity patterns. In special, the algorithms
- * do not try to exploit topological knowledge (e.g. corner detection) to
- * find appropriate starting points. This way, however, they work in
- * arbitrary space dimension.
- *
- * If you want to give starting points, you may give a list of dof indices
- * which will form the first step of the renumbering. The dofs of the list
- * will be consecutively numbered starting with zero, i.e. this list is not
- * renumbered according to the coordination number of the nodes. Indices not
- * in the allowed range are deleted. If no index is allowed, the algorithm
- * will search for its own starting point.
- *
- * 
- * \subsection{Results of renumbering}
- *
- * The renumbering schemes mentioned above do not lead to optimal results.
- * However, after all there is no algorithm that accomplishes this within
- * reasonable time. There are situations where the lack of optimality even
- * leads to worse results than with the original, crude, levelwise numering
- * scheme; one of these examples is a mesh of four cells of which always
- * those cells are refined which are neighbors to the center (you may call
- * this mesh a `zoom in' mesh). In one such example the bandwidth was
- * increased by about 50 per cent.
- *
- * In most other cases, the bandwith is reduced significantly. The reduction
- * is the better the less structured the grid is. With one grid where the
- * cells were refined according to a random driven algorithm, the bandwidth
- * was reduced by a factor of six.
- *
- * Using the constraint information usually leads to reductions in bandwidth
- * of 10 or 20 per cent, but may for some very unstructured grids also lead
- * to an increase. You have to weigh the decrease in your case with the time
- * spent to use the constraint information, which usually is several times
- * longer than the `pure' renumbering algorithm.
- *
- * In almost all cases, the renumbering scheme finds a corner to start with.
- * Since there is more than one corner in most grids and since even an
- * interior degree of freedom may be a better starting point, giving the
- * starting point by the user may be a viable way if you have a simple
- * scheme to derive a suitable point (e.g. by successively taking the
- * third child of the cell top left of the coarsest level, taking its
- * third vertex and the dof index thereof, if you want the top left corner
- * vertex). If you do not know beforehand what your grid will look like
- * (e.g. when using adaptive algorithms), searching a best starting point
- * may be difficult, however, and in many cases will not justify the effort.
+ * the #DoFRenumbering# class offers the function #renumber_dofs# which reorders
+ * the dof numbering according to some scheme. See there for a discussion of
+ * the implemented algorithms.
  *
  *
  * \subsection{User defined renumbering schemes}
  *
- * The #renumber_dofs# function offers a fixed number of renumbering
+ * The #DoFRenumbering# class offers a fixed number of renumbering
  * schemes like the Cuthill-McKey scheme. Basically, the function sets
  * up an array in which for each degree of freedom the index is stored
  * which is to be assigned by the renumbering. Using this array, the
@@ -550,22 +424,6 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      */
     virtual void clear ();
     
-                                    /**
-                                     * 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 this
-                                     * class for more details.
-                                     */
-    void renumber_dofs (const RenumberingMethod method,
-                       const bool use_constraints         = false,
-                       const vector<int> &starting_points = vector<int>());
-    
                                     /**
                                      * Actually do the renumbering based on
                                      * a list of new dof numbers for all the
@@ -1513,7 +1371,7 @@ class DoFHandler : public DoFDimensionInfo<dim> {
                                      * Exception
                                      */
     DeclException0 (ExcNoFESelected);
-                                    /**
+                                    /**
                                      * Exception
                                      */
     DeclException0 (ExcRenumberingIncomplete);
index 31e9920f5ea45ec2273d45582ad84aedbef245ab..36109d757c6f3433551c96d3315f70eb6c231391 100644 (file)
@@ -416,7 +416,7 @@ class MGDoFHexAccessor :  public MGDoFAccessor<dim>, public DoFHexAccessor<dim,
                                      * Return a pointer to the #i#th quad
                                      * bounding this #Hex#.
                                      */
-    TriaIterator<dim,MGDoFLineAccessor<dim,QuadAccessor<dim> > >
+    TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
     quad (const unsigned int i) const;
 
                                     /**
@@ -678,6 +678,8 @@ class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
 
 
 
+
+
 /*----------------------------   mg_dof_accessor.h     ---------------------------*/
 /* end of #ifndef __mg_dof_accessor_H */
 #endif
index edf50654fff2f29efc64beb422673582fc047613..1d6a7d7fb130d57095e2c5982253b7ade5f0a653 100644 (file)
@@ -131,6 +131,10 @@ class MGDoFHandler : public DoFHandler<dim> {
     typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
     typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
 
+    typedef typename MGDoFDimensionInfo<dim>::raw_hex_iterator raw_hex_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::hex_iterator hex_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_hex_iterator active_hex_iterator;
+
     typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
     typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
     typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
@@ -170,20 +174,18 @@ class MGDoFHandler : public DoFHandler<dim> {
     virtual void distribute_dofs (const FiniteElement<dim> &);
 
                                     /**
-                                     * Renumber the degrees of freedom according
-                                     * to the given scheme, eventually
-                                     * using the given
-                                     * starting points. The starting points
-                                     * default to an empty list, the use of
-                                     * constraint information defaults to
-                                     * false.
+                                     * Actually do the renumbering based on
+                                     * a list of new dof numbers for all the
+                                     * dofs. 
                                      *
-                                     * See the general documentation of the
-                                     * #DoFHandler# class for more details.
+                                     * #new_numbers# is an array of integers
+                                     * with size equal to the number of dofs
+                                     * on the present level. It stores the new
+                                     * indices after renumbering in the
+                                     * order of the old indices.
                                      */
-    void renumber_dofs (const unsigned int       level,
-                       const RenumberingMethod  method,
-                       const vector<int>       &starting_points = vector<int>());
+    void renumber_dofs (const unsigned int level,
+                       const vector<int> &new_numbers);
 
                                     /**
                                      * Write the sparsity structure of the
@@ -653,6 +655,103 @@ class MGDoFHandler : public DoFHandler<dim> {
 
                                     /*---------------------------------------*/
 
+                                    /**
+                                     *  @name Hex iterator functions*/
+                                    /*@{
+                                     */
+                                    /**
+                                     *  Return iterator to the first hex, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no hexs, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_hex_iterator    begin_raw_hex   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used hex
+                                     *  on level #level#.
+                                     */
+    hex_iterator        begin_hex       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  hex on level #level#.
+                                     */
+    active_hex_iterator begin_active_hex(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_hex_iterator    end_hex () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    hex_iterator        end_hex (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_hex_iterator    end_raw_hex (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_hex_iterator end_active_hex (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last hex, used or not.
+                                     */
+    raw_hex_iterator    last_raw_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  hex of the level #level#, used or not.
+
+                                     */
+    raw_hex_iterator    last_raw_hex (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used hex.
+                                     */
+    hex_iterator        last_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used hex on level #level#.
+                                     */
+    hex_iterator        last_hex (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active hex.
+                                     */
+    active_hex_iterator last_active_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active hex on level #level#.
+                                     */
+    active_hex_iterator last_active_hex (const unsigned int level) const;
+                                    /*@}*/
+    
+                                    /*---------------------------------------*/
+
                                     /**
                                      * Return the number of degrees of freedom
                                      * on the specified level.
@@ -831,22 +930,6 @@ class MGDoFHandler : public DoFHandler<dim> {
     unsigned int distribute_dofs_on_cell (cell_iterator &cell,
                                          unsigned int   next_free_dof);
     
-                                    /**
-                                     * Actually do the renumbering prepared
-                                     * by the #renumber_dofs# function on
-                                     * the given #level#. Since
-                                     * this is dimension specific, we
-                                     * need to have another function.
-                                     *
-                                     * #new_numbers# is an array of integers
-                                     * with size equal to the number of dofs
-                                     * on the present level. It stores the new
-                                     * indices after renumbering in the
-                                     * order of the old indices.
-                                     */
-    void do_renumbering (const unsigned int level,
-                        const vector<int> &new_numbers);
-
                                     /**
                                      * Reserve enough space for the MG dof
                                      * indices for a given triangulation.
@@ -881,6 +964,8 @@ class MGDoFHandler : public DoFHandler<dim> {
     friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
     friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
     friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+    friend class MGDoFHexAccessor<dim, HexAccessor<dim> >;
+    friend class MGDoFHexAccessor<dim, CellAccessor<dim> >;
 };
 
     
index 6035ebf6bcd1b5abaa2876d3b92b57fc1d877bb4..c1addddbd278eb682fd9284f9e750ec26c0266af 100644 (file)
@@ -473,6 +473,35 @@ HexAccessor<dim>::quad (const unsigned int i) const {
 
 
 
+template <int dim>
+inline
+unsigned int
+HexAccessor<dim>::line_index (unsigned int i) const {
+  Assert (i<12, ExcInvalidIndex(i,0,11));
+
+  if (i<4)
+    return quad(0)->line_index(i);
+  else
+    if (i<8)
+      return quad(1)->line_index(i-4);
+    else
+      switch (i) 
+       {
+         case 8:
+               return quad(2)->line_index(3);
+         case 9:
+               return quad(2)->line_index(1);
+         case 10:
+               return quad(4)->line_index(1);
+         case 11:
+               return quad(4)->line_index(3);
+       };
+  Assert (false, ExcInvalidIndex(i,0,11));
+  return 0;
+};
+
+
+
 template <int dim>
 inline
 unsigned int
index 31e9920f5ea45ec2273d45582ad84aedbef245ab..36109d757c6f3433551c96d3315f70eb6c231391 100644 (file)
@@ -416,7 +416,7 @@ class MGDoFHexAccessor :  public MGDoFAccessor<dim>, public DoFHexAccessor<dim,
                                      * Return a pointer to the #i#th quad
                                      * bounding this #Hex#.
                                      */
-    TriaIterator<dim,MGDoFLineAccessor<dim,QuadAccessor<dim> > >
+    TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
     quad (const unsigned int i) const;
 
                                     /**
@@ -678,6 +678,8 @@ class MGDoFCellAccessor :  public MGDoFSubstructAccessor<dim> {
 
 
 
+
+
 /*----------------------------   mg_dof_accessor.h     ---------------------------*/
 /* end of #ifndef __mg_dof_accessor_H */
 #endif
index edf50654fff2f29efc64beb422673582fc047613..1d6a7d7fb130d57095e2c5982253b7ade5f0a653 100644 (file)
@@ -131,6 +131,10 @@ class MGDoFHandler : public DoFHandler<dim> {
     typedef typename MGDoFDimensionInfo<dim>::quad_iterator quad_iterator;
     typedef typename MGDoFDimensionInfo<dim>::active_quad_iterator active_quad_iterator;
 
+    typedef typename MGDoFDimensionInfo<dim>::raw_hex_iterator raw_hex_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::hex_iterator hex_iterator;
+    typedef typename MGDoFDimensionInfo<dim>::active_hex_iterator active_hex_iterator;
+
     typedef typename MGDoFDimensionInfo<dim>::raw_cell_iterator raw_cell_iterator;
     typedef typename MGDoFDimensionInfo<dim>::cell_iterator cell_iterator;
     typedef typename MGDoFDimensionInfo<dim>::active_cell_iterator active_cell_iterator;
@@ -170,20 +174,18 @@ class MGDoFHandler : public DoFHandler<dim> {
     virtual void distribute_dofs (const FiniteElement<dim> &);
 
                                     /**
-                                     * Renumber the degrees of freedom according
-                                     * to the given scheme, eventually
-                                     * using the given
-                                     * starting points. The starting points
-                                     * default to an empty list, the use of
-                                     * constraint information defaults to
-                                     * false.
+                                     * Actually do the renumbering based on
+                                     * a list of new dof numbers for all the
+                                     * dofs. 
                                      *
-                                     * See the general documentation of the
-                                     * #DoFHandler# class for more details.
+                                     * #new_numbers# is an array of integers
+                                     * with size equal to the number of dofs
+                                     * on the present level. It stores the new
+                                     * indices after renumbering in the
+                                     * order of the old indices.
                                      */
-    void renumber_dofs (const unsigned int       level,
-                       const RenumberingMethod  method,
-                       const vector<int>       &starting_points = vector<int>());
+    void renumber_dofs (const unsigned int level,
+                       const vector<int> &new_numbers);
 
                                     /**
                                      * Write the sparsity structure of the
@@ -653,6 +655,103 @@ class MGDoFHandler : public DoFHandler<dim> {
 
                                     /*---------------------------------------*/
 
+                                    /**
+                                     *  @name Hex iterator functions*/
+                                    /*@{
+                                     */
+                                    /**
+                                     *  Return iterator to the first hex, used
+                                     *  or not, on level #level#. If a level
+                                     *  has no hexs, a past-the-end iterator
+                                     *  is returned.
+                                     */
+    raw_hex_iterator    begin_raw_hex   (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first used hex
+                                     *  on level #level#.
+                                     */
+    hex_iterator        begin_hex       (const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator to the first active
+                                     *  hex on level #level#.
+                                     */
+    active_hex_iterator begin_active_hex(const unsigned int level = 0) const;
+
+                                    /**
+                                     *  Return iterator past the end; this
+                                     *  iterator serves for comparisons of
+                                     *  iterators with past-the-end or
+                                     *  before-the-beginning states.
+                                     */
+    raw_hex_iterator    end_hex () const;
+
+                                    /**
+                                     * Return an iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    hex_iterator        end_hex (const unsigned int level) const;
+    
+                                    /**
+                                     * Return a raw iterator which is the first
+                                     * iterator not on level. If #level# is
+                                     * the last level, then this returns
+                                     * #end()#.
+                                     */
+    raw_hex_iterator    end_raw_hex (const unsigned int level) const;
+
+                                    /**
+                                     * Return an active iterator which is the
+                                     * first iterator not on level. If #level#
+                                     * is the last level, then this returns
+                                     * #end()#.
+                                     */
+    active_hex_iterator end_active_hex (const unsigned int level) const;
+
+
+                                    /**
+                                     *  Return an iterator pointing to the
+                                     *  last hex, used or not.
+                                     */
+    raw_hex_iterator    last_raw_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  hex of the level #level#, used or not.
+
+                                     */
+    raw_hex_iterator    last_raw_hex (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used hex.
+                                     */
+    hex_iterator        last_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  used hex on level #level#.
+                                     */
+    hex_iterator        last_hex (const unsigned int level) const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active hex.
+                                     */
+    active_hex_iterator last_active_hex () const;
+
+                                    /**
+                                     *  Return an iterator pointing to the last
+                                     *  active hex on level #level#.
+                                     */
+    active_hex_iterator last_active_hex (const unsigned int level) const;
+                                    /*@}*/
+    
+                                    /*---------------------------------------*/
+
                                     /**
                                      * Return the number of degrees of freedom
                                      * on the specified level.
@@ -831,22 +930,6 @@ class MGDoFHandler : public DoFHandler<dim> {
     unsigned int distribute_dofs_on_cell (cell_iterator &cell,
                                          unsigned int   next_free_dof);
     
-                                    /**
-                                     * Actually do the renumbering prepared
-                                     * by the #renumber_dofs# function on
-                                     * the given #level#. Since
-                                     * this is dimension specific, we
-                                     * need to have another function.
-                                     *
-                                     * #new_numbers# is an array of integers
-                                     * with size equal to the number of dofs
-                                     * on the present level. It stores the new
-                                     * indices after renumbering in the
-                                     * order of the old indices.
-                                     */
-    void do_renumbering (const unsigned int level,
-                        const vector<int> &new_numbers);
-
                                     /**
                                      * Reserve enough space for the MG dof
                                      * indices for a given triangulation.
@@ -881,6 +964,8 @@ class MGDoFHandler : public DoFHandler<dim> {
     friend class MGDoFLineAccessor<dim, CellAccessor<dim> >;
     friend class MGDoFQuadAccessor<dim, QuadAccessor<dim> >;
     friend class MGDoFQuadAccessor<dim, CellAccessor<dim> >;
+    friend class MGDoFHexAccessor<dim, HexAccessor<dim> >;
+    friend class MGDoFHexAccessor<dim, CellAccessor<dim> >;
 };
 
     
diff --git a/deal.II/deal.II/include/numerics/dof_renumbering.h b/deal.II/deal.II/include/numerics/dof_renumbering.h
new file mode 100644 (file)
index 0000000..62c748c
--- /dev/null
@@ -0,0 +1,205 @@
+/*----------------------------   dof_renumbering.h     ---------------------------*/
+/*      $Id$                 */
+#ifndef __dof_renumbering_H
+#define __dof_renumbering_H
+/*----------------------------   dof_renumbering.h     ---------------------------*/
+
+
+#include <base/exceptions.h>
+#include <basic/forward-declarations.h>
+
+
+
+
+/**
+ * Implementation of a number of renumbering algorithms for the degrees of
+ * freedom on a triangulation.
+ *
+ * \section{Cuthill-McKee like algorithms}
+ *
+ * Within this class, the Cuthill-McKee algorithm is implemented. It starts
+ * at a degree of freedom, searches the other DoFs for those which are couple
+ * with the one we started with and numbers these in a certain way. It then
+ * finds the second level of DoFs, namely those that couple with those of
+ * the previous level (which were those that coupled with the initial DoF)
+ * and numbers these. And so on. For the details of the algorithm, especially
+ * the numbering within each level, we refer the reader to the book of
+ * Schwarz (H.R.Schwarz: Methode der finiten Elemente). The reverse Cuthill-McKee
+ * algorithm does the same job, but numbers all elements in the reverse order.
+ *
+ * These algorithms
+ * have one major drawback: they require a good starting point, i.e. the degree
+ * of freedom index afterwards to be numbered zero. This can thus be given by
+ * the user, e.g. by exploiting knowledge of the actual topology of the
+ * domain. It is also possible to give several starting indices, which may
+ * be used to simulate a simple upstream numbering (by giving the inflow
+ * dofs as starting values) or to make preconditioning faster (by letting
+ * the dirichlet boundary indices be starting points).
+ *
+ * If no starting index is given, one is chosen by the program, namely one
+ * with the smallest coordination number (the coordination number is the
+ * number of other dofs this dof couples with). This dof is usually located
+ * on the boundary of the domain. There is, however, large ambiguity in this
+ * when using the hierarchical meshes used in this library, since in most
+ * cases the computational domain is not approximated by tilting and deforming
+ * elements and by plugging together variable numbers of elements at vertices,
+ * but rather by hierarchical refinement. There is therefore a large number
+ * of dofs with equal coordination numbers. The renumbering algorithms will
+ * therefore not give optimal results.
+ *
+ * In the book of Schwarz (H.R.Schwarz: Methode der finiten Elemente), it is
+ * advised to test many starting points, if possible all with the smallest
+ * coordination number and also those with slightly higher numbers. However,
+ * this seems only possible for meshes with at most several dozen or a few
+ * hundred elements found in small engineering problems of the early 1980s
+ * (the second edition was published in 1984), but certainly not with those
+ * used in this library, featuring several 10,000 to a few 100,000 elements.
+ *
+ * On the other hand, the need to reduce the bandwidth has decreased since
+ * with the mentioned number of cells, only iterative solution methods are
+ * able to solve the resulting matrix systems. These, however, are not so
+ * demanding with respect to the bandwidth as direct solvers used for
+ * smaller problems. Things like upstream numbering become much more important
+ * in recent times, so the suboptimality of the renumbering algorithms is
+ * not that important any more.
+ *
+ * 
+ * \subsection{Implementation of renumbering schemes}
+ *
+ * The renumbering algorithms need quite a lot of memory, since they have
+ * to store for each dof with which other dofs it couples. This is done
+ * using a #SparseMatrixStruct# object used to store the sparsity pattern of
+ * matrices. It
+ * is not useful for the user to do anything between distributing the dofs
+ * and renumbering, i.e. the calls to #DoFHandler::distribute_dofs# and
+ * #DoFHandler::renumber_dofs# should follow each other immediately. If
+ * you try to create a sparsity pattern or anything else in between, these
+ * will be invalid afterwards.
+ *
+ * The renumbering may take care of dof-to-dof couplings only induced by
+ * eliminating constraints. In addition to the memory consumption mentioned
+ * above, this also takes quite some computational time, but it may be
+ * switched off upon calling the #renumber_dofs# function. This will then
+ * give inferior results, since knots in the graph (representing dofs)
+ * are not found to be neighbors even if they would be after condensation.
+ * 
+ * The renumbering algorithms work on a purely algebraic basis, due to the
+ * isomorphism between the graph theoretical groundwork underlying the
+ * algorithms and binary matrices (matrices of which the entries are binary
+ * values) represented by the sparsity patterns. In special, the algorithms
+ * do not try to exploit topological knowledge (e.g. corner detection) to
+ * find appropriate starting points. This way, however, they work in
+ * arbitrary space dimension.
+ *
+ * If you want to give starting points, you may give a list of dof indices
+ * which will form the first step of the renumbering. The dofs of the list
+ * will be consecutively numbered starting with zero, i.e. this list is not
+ * renumbered according to the coordination number of the nodes. Indices not
+ * in the allowed range are deleted. If no index is allowed, the algorithm
+ * will search for its own starting point.
+ *
+ * 
+ * \subsection{Results of renumbering}
+ *
+ * The renumbering schemes mentioned above do not lead to optimal results.
+ * However, after all there is no algorithm that accomplishes this within
+ * reasonable time. There are situations where the lack of optimality even
+ * leads to worse results than with the original, crude, levelwise numering
+ * scheme; one of these examples is a mesh of four cells of which always
+ * those cells are refined which are neighbors to the center (you may call
+ * this mesh a `zoom in' mesh). In one such example the bandwidth was
+ * increased by about 50 per cent.
+ *
+ * In most other cases, the bandwith is reduced significantly. The reduction
+ * is the better the less structured the grid is. With one grid where the
+ * cells were refined according to a random driven algorithm, the bandwidth
+ * was reduced by a factor of six.
+ *
+ * Using the constraint information usually leads to reductions in bandwidth
+ * of 10 or 20 per cent, but may for some very unstructured grids also lead
+ * to an increase. You have to weigh the decrease in your case with the time
+ * spent to use the constraint information, which usually is several times
+ * longer than the `pure' renumbering algorithm.
+ *
+ * In almost all cases, the renumbering scheme finds a corner to start with.
+ * Since there is more than one corner in most grids and since even an
+ * interior degree of freedom may be a better starting point, giving the
+ * starting point by the user may be a viable way if you have a simple
+ * scheme to derive a suitable point (e.g. by successively taking the
+ * third child of the cell top left of the coarsest level, taking its
+ * third vertex and the dof index thereof, if you want the top left corner
+ * vertex). If you do not know beforehand what your grid will look like
+ * (e.g. when using adaptive algorithms), searching a best starting point
+ * may be difficult, however, and in many cases will not justify the effort.
+ *
+ *
+ * \subsection{Multigrid DoF numbering}
+ *
+ * Most algorithms also work on multigrid degree of freedom numberings. Refer
+ * to the actual function declarations to get more information on this.
+ *
+ *
+ * @author Wolfgang Bangerth, 1998, 1999
+ */
+class DoFRenumbering 
+{
+  public:
+                                    /**
+                                     * Renumber the degrees of freedom
+                                     * according to the Cuthill-McKee method,
+                                     * eventually using the reverse numbering
+                                     * scheme.
+                                     *
+                                     * See the general documentation of
+                                     * this class for details on the
+                                     * different methods.
+                                     */
+    template <int dim>
+    static void renumber_Cuthill_McKee (DoFHandler<dim>   &dof_handler,
+                                       const bool         reversed_numbering = false,
+                                       const bool         use_constraints    = false,
+                                       const vector<int> &starting_indices   = vector<int>());
+
+                                    /**
+                                     * Renumber the degrees of freedom
+                                     * according to the Cuthill-McKee method,
+                                     * eventually using the reverse numbering
+                                     * scheme, in this case for a multigrid
+                                     * numbering of degrees of freedom.
+                                     *
+                                     * You can give a triangulation level to
+                                     * which this function is to be applied.
+                                     * Since with a level-wise numbering there
+                                     * are no hanging nodes, no constraints
+                                     * can be used, so the respective
+                                     * parameter of the previous function is
+                                     * ommitted.
+                                     *
+                                     * See the general documentation of
+                                     * this class for details on the
+                                     * different methods.
+                                     */
+    template <int dim>
+    static void renumber_Cuthill_McKee (MGDoFHandler<dim> &dof_handler,
+                                       const unsigned int level,
+                                       const bool         reversed_numbering = false,
+                                       const vector<int> &starting_indices   = vector<int> ());
+
+    
+                                    /**
+                                     * Exception
+                                     */
+    DeclException0 (ExcRenumberingIncomplete);
+};
+
+
+
+
+
+
+
+
+/*----------------------------   dof_renumbering.h     ---------------------------*/
+/* end of #ifndef __dof_renumbering_H */
+#endif
+/*----------------------------   dof_renumbering.h     ---------------------------*/
index 982dae0a86f322613cd1375c29bc80d8482f1d26..fdb33a5dff89c6572615f652e5c0d406dc8117b4 100644 (file)
@@ -16,7 +16,6 @@
 #include <lac/vector.h>
 #include <lac/vector.h>
 
-#include <map>
 #include <set>
 #include <algorithm>
 
@@ -1565,201 +1564,6 @@ unsigned int DoFHandler<3>::distribute_dofs_on_cell (active_cell_iterator &cell,
 
 
 
-template <int dim>
-void DoFHandler<dim>::renumber_dofs (const RenumberingMethod method,
-                                    const bool use_constraints,
-                                    const vector<int> &starting_points) {
-                                  // make the connection graph
-  SparseMatrixStruct sparsity (n_dofs(), max_couplings_between_dofs());
-  make_sparsity_pattern (sparsity);
-
-  if (use_constraints) 
-    {
-      ConstraintMatrix constraints;
-      make_hanging_node_constraints (constraints);
-      constraints.close ();
-      constraints.condense (sparsity);
-    };
-    
-  int n_dofs = sparsity.n_rows();
-                                  // store the new dof numbers; -1 means
-                                  // that no new number was chosen yet
-  vector<int> new_number(sparsity.n_rows(), -1);
-  
-                                  // store the indices of the dofs renumbered
-                                  // in the last round. Default to starting
-                                  // points
-  vector<int> last_round_dofs (starting_points);
-  
-                                  // delete disallowed elements
-  for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-    if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
-      last_round_dofs[i] = -1;
-  
-  remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-            bind2nd(equal_to<int>(), -1));
-  
-                                  // now if no valid points remain:
-                                  // find dof with lowest coordination
-                                  // number
-  
-  if (last_round_dofs.size() == 0)
-    {
-      int          starting_point   = -1;
-      unsigned int min_coordination = n_dofs;
-      for (int row=0; row<n_dofs; ++row) 
-       {
-         unsigned int j;
-
-                                          // loop until we hit the end
-                                          // of this row's entries
-         for (j=sparsity.get_rowstart_indices()[row];
-              j<sparsity.get_rowstart_indices()[row+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-                                          // post-condition after loop:
-                                          // coordination, i.e. the number
-                                          // of entries in this row is now
-                                          // j-rowstart[row]
-         if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
-           {
-             min_coordination = j-sparsity.get_rowstart_indices()[row];
-             starting_point   = row;
-           };
-       };
-      
-                                      // now we still have to care for the
-                                      // case that no dof has a coordination
-                                      // number less than n_dofs. this rather
-                                      // exotic case only happens if we only
-                                      // have one cell, as far as I can see,
-                                      // but there may be others as well.
-                                      //
-                                      // if that should be the case, we can
-                                      // chose an arbitrary dof as starting
-                                      // point, e.g. the one with number zero
-      if (starting_point == -1)
-       starting_point = 0;
-      
-                                      // initialize the first dof
-      last_round_dofs.push_back (starting_point);
-    };
-  
-
-                                  // store next free dof index
-  int         next_free_number = 0;
-
-                                  // enumerate the first round dofs
-  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-    new_number[last_round_dofs[i]] = next_free_number++;
-
-  bool all_dofs_renumbered = false;
-
-                                  // now do as many steps as needed to
-                                  // renumber all dofs
-  while (!all_dofs_renumbered) 
-    {
-                                      // store the indices of the dofs to be
-                                      // renumbered in the next round
-      vector<int> next_round_dofs;
-
-                                      // find all neighbors of the
-                                      // dofs numbered in the last
-                                      // round
-      for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-       for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
-            j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-         if (sparsity.get_column_numbers()[j] == -1)
-           break;
-         else
-           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                      // sort dof numbers
-      sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                      // delete multiple entries
-      vector<int>::iterator end_sorted;
-      end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
-      next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
-                                      // eliminate dofs which are
-                                      // already numbered
-      for (int s=next_round_dofs.size()-1; s>=0; --s)
-       if (new_number[next_round_dofs[s]] != -1)
-         next_round_dofs.erase (&next_round_dofs[s]);
-
-                                      // check whether there are any new
-                                      // dofs
-      all_dofs_renumbered = (next_round_dofs.size() == 0);
-      if (all_dofs_renumbered)
-                                        // end loop if possible
-       continue;
-      
-
-                                      // store for each coordination
-                                      // number the dofs with these
-                                      // coordination number
-      multimap<unsigned int, int> dofs_by_coordination;
-      
-                                      // find coordination number for
-                                      // each of these dofs
-      for (vector<int>::iterator s=next_round_dofs.begin();
-          s!=next_round_dofs.end(); ++s) 
-       {
-         unsigned int coordination = 0;
-         for (unsigned int j=sparsity.get_rowstart_indices()[*s];
-              j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-           else
-             ++coordination;
-
-                                          // insert this dof at its
-                                          // coordination number
-         const pair<const unsigned int, int> new_entry (coordination, *s);
-         dofs_by_coordination.insert (new_entry);
-       };
-      
-                                      ////
-      multimap<unsigned int, int>::iterator i;
-      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-       new_number[i->second] = next_free_number++;
-
-                                      // after that: copy this round's
-                                      // dofs for the next round
-      last_round_dofs = next_round_dofs;
-    };
-
-#ifdef DEBUG
-                                  //  test for all indices numbered
-  if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
-    Assert (false, ExcRenumberingIncomplete());
-  Assert (next_free_number == n_dofs,
-         ExcRenumberingIncomplete());
-#endif
-
-  switch (method) 
-    {
-      case Cuthill_McKee:
-           break;
-      case reverse_Cuthill_McKee:
-      {
-       for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
-         *i = n_dofs-*i;
-       break;
-      };
-      default:
-           Assert (false, ExcNotImplemented());
-    };
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  renumber_dofs (new_number);
-};
-
-
-
 #if deal_II_dimension == 1
 
 template <>
diff --git a/deal.II/deal.II/source/dofs/dof_renumbering.cc b/deal.II/deal.II/source/dofs/dof_renumbering.cc
new file mode 100644 (file)
index 0000000..ca68898
--- /dev/null
@@ -0,0 +1,379 @@
+/* $Id$ */
+
+#include <grid/dof.h>
+#include <grid/mg_dof.h>
+#include <grid/dof_constraints.h>
+#include <numerics/dof_renumbering.h>
+#include <lac/sparsematrix.h>
+
+#include <vector>
+#include <map>
+#include <algorithm>
+
+
+
+template <int dim>
+void DoFRenumbering::renumber_Cuthill_McKee (DoFHandler<dim>   &dof_handler,
+                                            const bool         reversed_numbering,
+                                            const bool         use_constraints,
+                                            const vector<int> &starting_indices) {
+                                  // make the connection graph
+  SparseMatrixStruct sparsity (dof_handler.n_dofs(),
+                              dof_handler.max_couplings_between_dofs());
+  dof_handler.make_sparsity_pattern (sparsity);
+
+  if (use_constraints) 
+    {
+      ConstraintMatrix constraints;
+      dof_handler.make_hanging_node_constraints (constraints);
+      constraints.close ();
+      constraints.condense (sparsity);
+    };
+    
+  const int n_dofs = sparsity.n_rows();
+                                  // store the new dof numbers; -1 means
+                                  // that no new number was chosen yet
+  vector<int> new_number(sparsity.n_rows(), -1);
+  
+                                  // store the indices of the dofs renumbered
+                                  // in the last round. Default to starting
+                                  // points
+  vector<int> last_round_dofs (starting_indices);
+  
+                                  // delete disallowed elements
+  for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+    if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+      last_round_dofs[i] = -1;
+  
+  remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+            bind2nd(equal_to<int>(), -1));
+  
+                                  // now if no valid points remain:
+                                  // find dof with lowest coordination
+                                  // number
+  
+  if (last_round_dofs.size() == 0)
+    {
+      int          starting_point   = -1;
+      unsigned int min_coordination = n_dofs;
+      for (int row=0; row<n_dofs; ++row) 
+       {
+         unsigned int j;
+
+                                          // loop until we hit the end
+                                          // of this row's entries
+         for (j=sparsity.get_rowstart_indices()[row];
+              j<sparsity.get_rowstart_indices()[row+1]; ++j)
+           if (sparsity.get_column_numbers()[j] == -1)
+             break;
+                                          // post-condition after loop:
+                                          // coordination, i.e. the number
+                                          // of entries in this row is now
+                                          // j-rowstart[row]
+         if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
+           {
+             min_coordination = j-sparsity.get_rowstart_indices()[row];
+             starting_point   = row;
+           };
+       };
+      
+                                      // now we still have to care for the
+                                      // case that no dof has a coordination
+                                      // number less than n_dofs. this rather
+                                      // exotic case only happens if we only
+                                      // have one cell, as far as I can see,
+                                      // but there may be others as well.
+                                      //
+                                      // if that should be the case, we can
+                                      // chose an arbitrary dof as starting
+                                      // point, e.g. the one with number zero
+      if (starting_point == -1)
+       starting_point = 0;
+      
+                                      // initialize the first dof
+      last_round_dofs.push_back (starting_point);
+    };
+  
+
+                                  // store next free dof index
+  int         next_free_number = 0;
+
+                                  // enumerate the first round dofs
+  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+    new_number[last_round_dofs[i]] = next_free_number++;
+
+  bool all_dofs_renumbered = false;
+
+                                  // now do as many steps as needed to
+                                  // renumber all dofs
+  while (!all_dofs_renumbered) 
+    {
+                                      // store the indices of the dofs to be
+                                      // renumbered in the next round
+      vector<int> next_round_dofs;
+
+                                      // find all neighbors of the
+                                      // dofs numbered in the last
+                                      // round
+      for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+       for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
+            j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
+         if (sparsity.get_column_numbers()[j] == -1)
+           break;
+         else
+           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+      
+                                      // sort dof numbers
+      sort (next_round_dofs.begin(), next_round_dofs.end());
+
+                                      // delete multiple entries
+      vector<int>::iterator end_sorted;
+      end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
+      next_round_dofs.erase (end_sorted, next_round_dofs.end());
+
+                                      // eliminate dofs which are
+                                      // already numbered
+      for (int s=next_round_dofs.size()-1; s>=0; --s)
+       if (new_number[next_round_dofs[s]] != -1)
+         next_round_dofs.erase (&next_round_dofs[s]);
+
+                                      // check whether there are any new
+                                      // dofs
+      all_dofs_renumbered = (next_round_dofs.size() == 0);
+      if (all_dofs_renumbered)
+                                        // end loop if possible
+       continue;
+      
+
+                                      // store for each coordination
+                                      // number the dofs with these
+                                      // coordination number
+      multimap<unsigned int, int> dofs_by_coordination;
+      
+                                      // find coordination number for
+                                      // each of these dofs
+      for (vector<int>::iterator s=next_round_dofs.begin();
+          s!=next_round_dofs.end(); ++s) 
+       {
+         unsigned int coordination = 0;
+         for (unsigned int j=sparsity.get_rowstart_indices()[*s];
+              j<sparsity.get_rowstart_indices()[*s+1]; ++j)
+           if (sparsity.get_column_numbers()[j] == -1)
+             break;
+           else
+             ++coordination;
+
+                                          // insert this dof at its
+                                          // coordination number
+         const pair<const unsigned int, int> new_entry (coordination, *s);
+         dofs_by_coordination.insert (new_entry);
+       };
+      
+                                      ////
+      multimap<unsigned int, int>::iterator i;
+      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
+       new_number[i->second] = next_free_number++;
+
+                                      // after that: copy this round's
+                                      // dofs for the next round
+      last_round_dofs = next_round_dofs;
+    };
+
+#ifdef DEBUG
+                                  //  test for all indices numbered
+  if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+    Assert (false, ExcRenumberingIncomplete());
+  Assert (next_free_number == n_dofs,
+         ExcRenumberingIncomplete());
+#endif
+
+  if (reversed_numbering)
+    for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+      *i = n_dofs-*i;
+
+                                  // actually perform renumbering;
+                                  // this is dimension specific and
+                                  // thus needs an own function
+  dof_handler.renumber_dofs (new_number);
+};
+
+
+
+
+template <int dim>
+void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler<dim>      &dof_handler,
+                                            const unsigned int      level,
+                                            const bool              reversed_numbering,
+                                            const vector<int>      &starting_indices) {
+                                  // make the connection graph
+  SparseMatrixStruct sparsity (dof_handler.n_dofs(level),
+                              dof_handler.max_couplings_between_dofs());
+  dof_handler.make_sparsity_pattern (level, sparsity);
+    
+  const int n_dofs = sparsity.n_rows();
+                                  // store the new dof numbers; -1 means
+                                  // that no new number was chosen yet
+  vector<int> new_number(n_dofs, -1);
+  
+                                  // store the indices of the dofs renumbered
+                                  // in the last round. Default to starting
+                                  // points
+  vector<int> last_round_dofs (starting_indices);
+  
+                                  // delete disallowed elements
+  for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+    if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
+      last_round_dofs[i] = -1;
+  
+  remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+            bind2nd(equal_to<int>(), -1));
+  
+                                  // now if no valid points remain:
+                                  // find dof with lowest coordination
+                                  // number
+  
+  if (last_round_dofs.size() == 0)
+    {
+      int          starting_point   = -1;
+      unsigned int min_coordination = n_dofs;
+      for (int row=0; row<n_dofs; ++row) 
+       {
+         unsigned int j;
+         for (j=sparsity.get_rowstart_indices()[row];
+              j<sparsity.get_rowstart_indices()[row+1]; ++j)
+           if (sparsity.get_column_numbers()[j] == -1)
+             break;
+                                          // post-condition after loop:
+                                          // coordination is now
+                                          // j-rowstart[row]
+         if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
+           {
+             min_coordination = j-sparsity.get_rowstart_indices()[row];
+             starting_point   = row;
+           };
+       };
+                                      // initialize the first dof
+      last_round_dofs.push_back (starting_point);
+    };
+  
+
+                                  // store next free dof index
+  int         next_free_number = 0;
+
+                                  // enumerate the first round dofs
+  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+    new_number[last_round_dofs[i]] = next_free_number++;
+
+  bool all_dofs_renumbered = false;
+
+                                  // now do as many steps as needed to
+                                  // renumber all dofs
+  while (!all_dofs_renumbered) 
+    {
+                                      // store the indices of the dofs to be
+                                      // renumbered in the next round
+      vector<int> next_round_dofs;
+
+                                      // find all neighbors of the
+                                      // dofs numbered in the last
+                                      // round
+      for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+       for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
+            j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
+         if (sparsity.get_column_numbers()[j] == -1)
+           break;
+         else
+           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
+      
+                                      // sort dof numbers
+      sort (next_round_dofs.begin(), next_round_dofs.end());
+
+                                      // delete multiple entries
+      vector<int>::iterator end_sorted;
+      end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
+      next_round_dofs.erase (end_sorted, next_round_dofs.end());
+
+                                      // eliminate dofs which are
+                                      // already numbered
+      for (int s=next_round_dofs.size()-1; s>=0; --s)
+       if (new_number[next_round_dofs[s]] != -1)
+         next_round_dofs.erase (&next_round_dofs[s]);
+
+                                      // check whether there are any new
+                                      // dofs
+      all_dofs_renumbered = (next_round_dofs.size() == 0);
+      if (all_dofs_renumbered)
+                                        // end loop if possible
+       continue;
+      
+
+                                      // store for each coordination
+                                      // number the dofs with these
+                                      // coordination number
+      multimap<unsigned int, int> dofs_by_coordination;
+      
+                                      // find coordination number for
+                                      // each of these dofs
+      for (vector<int>::iterator s=next_round_dofs.begin();
+          s!=next_round_dofs.end(); ++s) 
+       {
+         unsigned int coordination = 0;
+         for (unsigned int j=sparsity.get_rowstart_indices()[*s];
+              j<sparsity.get_rowstart_indices()[*s+1]; ++j)
+           if (sparsity.get_column_numbers()[j] == -1)
+             break;
+           else
+             ++coordination;
+
+                                          // insert this dof at its
+                                          // coordination number
+         const pair<const unsigned int, int> new_entry (coordination, *s);
+         dofs_by_coordination.insert (new_entry);
+       };
+      
+                                      ////
+      multimap<unsigned int, int>::iterator i;
+      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
+       new_number[i->second] = next_free_number++;
+
+                                      // after that: copy this round's
+                                      // dofs for the next round
+      last_round_dofs = next_round_dofs;
+    };
+
+#ifdef DEBUG
+                                  //  test for all indices numbered
+  if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
+    Assert (false, ExcRenumberingIncomplete());
+  Assert (next_free_number == n_dofs,
+         ExcRenumberingIncomplete());
+#endif
+
+  if (reversed_numbering)
+    for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
+      *i = n_dofs-*i;
+
+                                  // actually perform renumbering;
+                                  // this is dimension specific and
+                                  // thus needs an own function
+  dof_handler.renumber_dofs (level, new_number);
+};
+
+
+
+
+
+
+// explicit instantiations
+template
+void DoFRenumbering::renumber_Cuthill_McKee (DoFHandler<deal_II_dimension> &dof_handler,
+                                            const bool                     reversed_numbering,
+                                            const bool                     use_constraints,
+                                            const vector<int>             &starting_indices);
+
+template
+void DoFRenumbering::renumber_Cuthill_McKee (MGDoFHandler<deal_II_dimension> &dof_handler,
+                                            const unsigned int               level,
+                                            const bool                       reversed_numbering,
+                                            const vector<int>               &starting_indices);
+
index 101f9344fea71eaa9415595b2296e732d5d8d9e3..d6fd3b0b2445aa06160a9d3d78252609f555d2c9 100644 (file)
@@ -299,6 +299,177 @@ MGDoFQuadAccessor<dim,BaseClass>::copy_from (const MGDoFQuadAccessor<dim,BaseCla
 
 
 
+/* ------------------------ MGDoFHexAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFHexAccessor<dim,BaseClass>::MGDoFHexAccessor (Triangulation<dim> *tria,
+                                                  const int           level,
+                                                  const int           index,
+                                                  const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFHexAccessor(tria,level,index,local_data) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_fe().dofs_per_hex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_fe().dofs_per_hex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (vertex<8, ExcInvalidIndex (i,0,8));
+  Assert (i<dof_handler->get_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+  
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex +
+                                12*dof_handler->get_fe().dofs_per_line +
+                                6*dof_handler->get_fe().dofs_per_quad +
+                                dof_handler->get_fe().dofs_per_hex),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_fe().dofs_per_quad,
+                    dofs_per_hex    = dof_handler->get_fe().dofs_per_hex;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<8; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<12; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->mg_dof_index(d);
+  for (unsigned int quad=0; quad<12; ++quad)
+    for (unsigned int d=0; d<dofs_per_quad; ++d)
+      *next++ = this->quad(quad)->mg_dof_index(d);
+  for (unsigned int d=0; d<dofs_per_hex; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::line (const unsigned int i) const {
+  Assert (i<12, ExcInvalidIndex (i, 0, 12));
+
+  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,QuadAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::quad (const unsigned int i) const {
+  Assert (i<12, ExcInvalidIndex (i, 0, 6));
+
+  return TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
+    (
+      tria,
+      present_level,
+      quad_index (i),
+      mg_dof_handler
+    );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFHexAccessor<dim,BaseClass> >
+MGDoFHexAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFHexAccessor<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;
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<dim,BaseClass>::copy_from (const MGDoFHexAccessor<dim,BaseClass> &a) {
+  DoFHexAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
+
 
 /*------------------------- Functions: MGDoFCellAccessor -----------------------*/
 
@@ -440,3 +611,23 @@ template class TriaActiveIterator<2,MGDoFCellAccessor<2> >;
 #endif
 
 
+#if deal_II_dimension == 3
+template class MGDoFLineAccessor<3,LineAccessor<3> >;
+template class MGDoFQuadAccessor<3,QuadAccessor<3> >;
+template class MGDoFHexAccessor<3,HexAccessor<3> >;
+template class MGDoFHexAccessor<3,CellAccessor<3> >;
+template class MGDoFCellAccessor<3>;
+
+template class TriaRawIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFHexAccessor<3,HexAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFCellAccessor<3> >;
+template class TriaIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaIterator<3,MGDoFCellAccessor<3> >;
+template class TriaActiveIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFCellAccessor<3> >;
+#endif
+
+
index 2d6f4ffaf4a2c8c8f584460ea04e8bd1631ccf27..26c8d83933ce6bdaa4b088e59f75679f53b93b2f 100644 (file)
@@ -353,6 +353,94 @@ MGDoFHandler<1>::last_active_quad () const {
   return 0;
 };
 
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::begin_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::begin_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::begin_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::end_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
 #endif
 
 
@@ -516,6 +604,261 @@ MGDoFHandler<2>::last_active_face (const unsigned int level) const {
   return last_active_line (level);
 };
 
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::begin_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::begin_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::begin_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::end_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+#endif
+
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::begin_raw (const unsigned int level) const {
+  return begin_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::begin (const unsigned int level) const {
+  return begin_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::begin_active (const unsigned int level) const {
+  return begin_active_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::end () const {
+  return end_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw () const {
+  return last_raw_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw (const unsigned int level) const {
+  return last_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last () const {
+  return last_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last (const unsigned int level) const {
+  return last_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active () const {
+  return last_active_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active (const unsigned int level) const {
+  return last_active_hex (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::begin_face (const unsigned int level) const {
+  return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::begin_active_face (const unsigned int level) const {
+  return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::end_face () const {
+  return end_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face () const {
+  return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face (const unsigned int level) const {
+  return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face () const {
+  return last_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face (const unsigned int level) const {
+  return last_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face () const {
+  return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face (const unsigned int level) const {
+  return last_active_quad (level);
+};
+
 #endif
 
 
@@ -588,6 +931,39 @@ MGDoFHandler<dim>::begin_active_quad (const unsigned int level) const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::begin_raw_hex (const unsigned int level) const {
+  return raw_hex_iterator (tria,
+                          tria->begin_raw_hex(level)->level(),
+                          tria->begin_raw_hex(level)->index(),
+                          this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::begin_hex (const unsigned int level) const {
+  return hex_iterator (tria,
+                      tria->begin_hex(level)->level(),
+                      tria->begin_hex(level)->index(),
+                      this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::begin_active_hex (const unsigned int level) const {
+  return active_hex_iterator (tria,
+                             tria->begin_active_hex(level)->level(),
+                             tria->begin_active_hex(level)->index(),
+                             this);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::raw_line_iterator
 MGDoFHandler<dim>::end_line () const {
@@ -604,6 +980,14 @@ MGDoFHandler<dim>::end_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_hex () const {
+  return raw_hex_iterator (tria, -1, -1, this);
+};
+
+
+
 template <int dim>
 typename MGDoFDimensionInfo<dim>::raw_cell_iterator
 MGDoFHandler<dim>::end_raw (const unsigned int level) const {
@@ -694,34 +1078,65 @@ MGDoFHandler<dim>::end_active_line (const unsigned int level) const {
 
 
 
-
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_quad() :
+         begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         quad_iterator(end_quad()) :
+         begin_quad (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_quad_iterator(end_quad()) :
+         begin_active_quad (level+1));
+};
+
+
+
 template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_raw_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         end_quad() :
-         begin_raw_quad (level+1));
+         end_hex() :
+         begin_raw_hex (level+1));
 };
 
 
 
 template <int dim>
-typename MGDoFDimensionInfo<dim>::quad_iterator
-MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::hex_iterator
+MGDoFHandler<dim>::end_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         quad_iterator(end_quad()) :
-         begin_quad (level+1));
+         hex_iterator(end_hex()) :
+         begin_hex (level+1));
 };
 
 
 
 
 template <int dim>
-typename MGDoFDimensionInfo<dim>::active_quad_iterator
-MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::active_hex_iterator
+MGDoFHandler<dim>::end_active_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         active_quad_iterator(end_quad()) :
-         begin_active_quad (level+1));
+         active_hex_iterator(end_hex()) :
+         begin_active_hex (level+1));
 };
 
 
@@ -795,6 +1210,42 @@ MGDoFHandler<dim>::last_active_quad (const unsigned int level) const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex (const unsigned int level) const {
+  return raw_hex_iterator (tria,
+                          tria->last_raw_hex(level)->level(),
+                          tria->last_raw_hex(level)->index(),
+                          this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex (const unsigned int level) const {
+  return hex_iterator (tria,
+                      tria->last_hex(level)->level(),
+                      tria->last_hex(level)->index(),
+                      this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex (const unsigned int level) const {
+  return active_hex_iterator (tria,
+                             tria->last_active_hex(level)->level(),
+                             tria->last_active_hex(level)->index(),
+                             this);
+};
+
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::raw_line_iterator
 MGDoFHandler<dim>::last_raw_line () const {
@@ -811,6 +1262,14 @@ MGDoFHandler<dim>::last_raw_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex () const {
+  return last_raw_hex (mg_levels.size()-1);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::line_iterator
 MGDoFHandler<dim>::last_line () const {
@@ -827,6 +1286,14 @@ MGDoFHandler<dim>::last_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex () const {
+  return last_hex (mg_levels.size()-1);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::active_line_iterator
 MGDoFHandler<dim>::last_active_line () const {
@@ -844,6 +1311,14 @@ MGDoFHandler<dim>::last_active_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex () const {
+  return last_active_hex (mg_levels.size()-1);
+};
+
+
+
 
 
 
@@ -989,6 +1464,69 @@ MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell,
 
 
 
+#if deal_II_dimension == 3
+
+template <>
+unsigned int
+MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell,
+                                         unsigned int   next_free_dof) {
+  if (selected_fe->dofs_per_vertex > 0)
+                                    // number dofs on vertices
+    for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+                                      // check whether dofs for this
+                                      // vertex have been distributed
+                                      // (only check the first dof)
+      if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+    
+                                  // for the lines
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+      {
+       line_iterator line = cell->line(l);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (line->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           line->set_mg_dof_index (d, next_free_dof++);            
+      };
+
+                                  // for the quads
+  if (selected_fe->dofs_per_quad > 0)
+    for (unsigned int q=0; q<GeometryInfo<3>::lines_per_cell; ++q)
+      {
+       quad_iterator quad = cell->quad(q);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (quad->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           quad->set_mg_dof_index (d, next_free_dof++);            
+      };
+  
+
+                                      // dofs of cell
+  if (selected_fe->dofs_per_hex > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_hex; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+  
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
+  
+  return next_free_dof;
+};
+
+#endif
+
+
+
 
 template <int dim>
 unsigned int MGDoFHandler<dim>::n_dofs (const unsigned int level) const {
@@ -1025,182 +1563,11 @@ void MGDoFHandler<dim>::make_sparsity_pattern (const unsigned int  level,
 
 
 
-
-template <int dim>
-void MGDoFHandler<dim>::renumber_dofs (const unsigned int      level,
-                                      const RenumberingMethod method,
-                                      const vector<int>      &starting_points) {
-                                  // make the connection graph
-  SparseMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
-  make_sparsity_pattern (level, sparsity);
-    
-  int n_dofs = sparsity.n_rows();
-                                  // store the new dof numbers; -1 means
-                                  // that no new number was chosen yet
-  vector<int> new_number(n_dofs, -1);
-  
-                                  // store the indices of the dofs renumbered
-                                  // in the last round. Default to starting
-                                  // points
-  vector<int> last_round_dofs (starting_points);
-  
-                                  // delete disallowed elements
-  for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-    if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
-      last_round_dofs[i] = -1;
-  
-  remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-            bind2nd(equal_to<int>(), -1));
-  
-                                  // now if no valid points remain:
-                                  // find dof with lowest coordination
-                                  // number
-  
-  if (last_round_dofs.size() == 0)
-    {
-      int          starting_point   = -1;
-      unsigned int min_coordination = n_dofs;
-      for (int row=0; row<n_dofs; ++row) 
-       {
-         unsigned int j;
-         for (j=sparsity.get_rowstart_indices()[row];
-              j<sparsity.get_rowstart_indices()[row+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-                                          // post-condition after loop:
-                                          // coordination is now
-                                          // j-rowstart[row]
-         if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
-           {
-             min_coordination = j-sparsity.get_rowstart_indices()[row];
-             starting_point   = row;
-           };
-       };
-                                      // initialize the first dof
-      last_round_dofs.push_back (starting_point);
-    };
-  
-
-                                  // store next free dof index
-  int         next_free_number = 0;
-
-                                  // enumerate the first round dofs
-  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-    new_number[last_round_dofs[i]] = next_free_number++;
-
-  bool all_dofs_renumbered = false;
-
-                                  // now do as many steps as needed to
-                                  // renumber all dofs
-  while (!all_dofs_renumbered) 
-    {
-                                      // store the indices of the dofs to be
-                                      // renumbered in the next round
-      vector<int> next_round_dofs;
-
-                                      // find all neighbors of the
-                                      // dofs numbered in the last
-                                      // round
-      for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-       for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
-            j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-         if (sparsity.get_column_numbers()[j] == -1)
-           break;
-         else
-           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                      // sort dof numbers
-      sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                      // delete multiple entries
-      vector<int>::iterator end_sorted;
-      end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
-      next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
-                                      // eliminate dofs which are
-                                      // already numbered
-      for (int s=next_round_dofs.size()-1; s>=0; --s)
-       if (new_number[next_round_dofs[s]] != -1)
-         next_round_dofs.erase (&next_round_dofs[s]);
-
-                                      // check whether there are any new
-                                      // dofs
-      all_dofs_renumbered = (next_round_dofs.size() == 0);
-      if (all_dofs_renumbered)
-                                        // end loop if possible
-       continue;
-      
-
-                                      // store for each coordination
-                                      // number the dofs with these
-                                      // coordination number
-      multimap<unsigned int, int> dofs_by_coordination;
-      
-                                      // find coordination number for
-                                      // each of these dofs
-      for (vector<int>::iterator s=next_round_dofs.begin();
-          s!=next_round_dofs.end(); ++s) 
-       {
-         unsigned int coordination = 0;
-         for (unsigned int j=sparsity.get_rowstart_indices()[*s];
-              j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-           else
-             ++coordination;
-
-                                          // insert this dof at its
-                                          // coordination number
-         const pair<const unsigned int, int> new_entry (coordination, *s);
-         dofs_by_coordination.insert (new_entry);
-       };
-      
-                                      ////
-      multimap<unsigned int, int>::iterator i;
-      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-       new_number[i->second] = next_free_number++;
-
-                                      // after that: copy this round's
-                                      // dofs for the next round
-      last_round_dofs = next_round_dofs;
-    };
-
-#ifdef DEBUG
-                                  //  test for all indices numbered
-  if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
-    Assert (false, ExcRenumberingIncomplete());
-  Assert (next_free_number == n_dofs,
-         ExcRenumberingIncomplete());
-#endif
-
-  switch (method) 
-    {
-      case Cuthill_McKee:
-           break;
-      case reverse_Cuthill_McKee:
-      {
-       for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
-         *i = n_dofs-*i;
-       break;
-      };
-      default:
-           Assert (false, ExcNotImplemented());
-    };
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  do_renumbering (level, new_number);
-};
-
-
-
-
 #if deal_II_dimension == 1
 
 template <>
-void MGDoFHandler<1>::do_renumbering (const unsigned int level,
-                                     const vector<int> &new_numbers) {
+void MGDoFHandler<1>::renumber_dofs (const unsigned int level,
+                                    const vector<int> &new_numbers) {
   Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
   
                                   // note that we can not use cell iterators
@@ -1236,8 +1603,45 @@ void MGDoFHandler<1>::do_renumbering (const unsigned int level,
 #if deal_II_dimension == 2
 
 template <>
-void MGDoFHandler<2>::do_renumbering (const unsigned int  level,
-                                     const vector<int>  &new_numbers) {
+void MGDoFHandler<2>::renumber_dofs (const unsigned int  level,
+                                    const vector<int>  &new_numbers) {
+  Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+  
+  for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+       i!=mg_vertex_dofs.end(); ++i)
+                                    // if the present vertex lives on
+                                    // the present level
+    if ((i->get_coarsest_level() <= level) &&
+       (i->get_finest_level() >= level))
+      for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+       i->set_index (level, d, selected_fe->dofs_per_vertex,
+                     new_numbers[i->get_index (level, d,
+                                               selected_fe->dofs_per_vertex)]);
+  
+  for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+       i!=mg_levels[level]->line_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
+
+  for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+       i!=mg_levels[level]->quad_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::renumber_dofs (const unsigned int  level,
+                                    const vector<int>  &new_numbers) {
   Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
   
   for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
@@ -1257,12 +1661,20 @@ void MGDoFHandler<2>::do_renumbering (const unsigned int  level,
       Assert (*i != -1, ExcInternalError());
       *i = new_numbers[*i];
     };
+
   for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
        i!=mg_levels[level]->quad_dofs.end(); ++i)
     {
       Assert (*i != -1, ExcInternalError());
       *i = new_numbers[*i];
     };
+
+  for (vector<int>::iterator i=mg_levels[level]->hex_dofs.begin();
+       i!=mg_levels[level]->hex_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
 };
 
 #endif
@@ -1445,5 +1857,100 @@ void MGDoFHandler<2>::reserve_space () {
 #endif
 
 
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::reserve_space () {
+  const unsigned int dim = 3;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+  
+                                  ////////////////////////////
+                                  // DESTRUCTION
+
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.clear ();
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.clear ();
+  
+  
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines and quads on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<3>);
+
+      mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+                                                selected_fe->dofs_per_line,
+                                                -1);
+      mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+                                                selected_fe->dofs_per_quad,
+                                                -1);
+      mg_levels.back()->hex_dofs = vector<int> (tria->levels[i]->hexes.hexes.size() *
+                                               selected_fe->dofs_per_hex,
+                                               -1);
+    };
+
+  
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+  
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex],
+                                  selected_fe->dofs_per_vertex);
+    };
+};
+
+#endif
+
+
 // explicite instantiations
 template class MGDoFHandler<deal_II_dimension>;
index 101f9344fea71eaa9415595b2296e732d5d8d9e3..d6fd3b0b2445aa06160a9d3d78252609f555d2c9 100644 (file)
@@ -299,6 +299,177 @@ MGDoFQuadAccessor<dim,BaseClass>::copy_from (const MGDoFQuadAccessor<dim,BaseCla
 
 
 
+/* ------------------------ MGDoFHexAccessor --------------------------- */
+
+template <int dim, typename BaseClass>
+MGDoFHexAccessor<dim,BaseClass>::MGDoFHexAccessor (Triangulation<dim> *tria,
+                                                  const int           level,
+                                                  const int           index,
+                                                  const AccessorData *local_data) :
+               MGDoFAccessor<dim> (local_data),
+               DoFHexAccessor(tria,level,index,local_data) {};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_fe().dofs_per_hex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+  return mg_dof_handler->mg_levels[present_level]
+    ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i];
+};
+
+
+
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (i<dof_handler->get_fe().dofs_per_hex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_hex));
+
+  mg_dof_handler->mg_levels[present_level]
+    ->hex_dofs[present_index*dof_handler->get_fe().dofs_per_hex+i] = index;
+};
+
+
+
+template <int dim, typename BaseClass>
+inline
+int MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (vertex<8, ExcInvalidIndex (i,0,8));
+  Assert (i<dof_handler->get_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+  
+  return (mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+         .get_index (present_level, i, dof_handler->get_fe().dofs_per_vertex));
+};
+
+
+  
+template <int dim, typename BaseClass>
+void MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (vertex<4, ExcInvalidIndex (i,0,4));
+  Assert (i<dof_handler->get_fe().dofs_per_vertex,
+         ExcInvalidIndex (i, 0, dof_handler->get_fe().dofs_per_vertex));
+
+  mg_dof_handler->mg_vertex_dofs[vertex_index(vertex)]
+    .set_index (present_level, i, dof_handler->get_fe().dofs_per_vertex, index);
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<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_fe() != 0, ExcInvalidObject());
+  Assert (dof_indices.size() == (8*dof_handler->get_fe().dofs_per_vertex +
+                                12*dof_handler->get_fe().dofs_per_line +
+                                6*dof_handler->get_fe().dofs_per_quad +
+                                dof_handler->get_fe().dofs_per_hex),
+         ExcVectorDoesNotMatch());
+
+  const unsigned int dofs_per_vertex = dof_handler->get_fe().dofs_per_vertex,
+                    dofs_per_line   = dof_handler->get_fe().dofs_per_line,
+                    dofs_per_quad   = dof_handler->get_fe().dofs_per_quad,
+                    dofs_per_hex    = dof_handler->get_fe().dofs_per_hex;
+  vector<int>::iterator next = dof_indices.begin();
+  for (unsigned int vertex=0; vertex<8; ++vertex)
+    for (unsigned int d=0; d<dofs_per_vertex; ++d)
+      *next++ = mg_vertex_dof_index(vertex,d);
+  for (unsigned int line=0; line<12; ++line)
+    for (unsigned int d=0; d<dofs_per_line; ++d)
+      *next++ = this->line(line)->mg_dof_index(d);
+  for (unsigned int quad=0; quad<12; ++quad)
+    for (unsigned int d=0; d<dofs_per_quad; ++d)
+      *next++ = this->quad(quad)->mg_dof_index(d);
+  for (unsigned int d=0; d<dofs_per_hex; ++d)
+    *next++ = mg_dof_index(d);
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFLineAccessor<dim,LineAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::line (const unsigned int i) const {
+  Assert (i<12, ExcInvalidIndex (i, 0, 12));
+
+  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,QuadAccessor<dim> > >
+MGDoFHexAccessor<dim,BaseClass>::quad (const unsigned int i) const {
+  Assert (i<12, ExcInvalidIndex (i, 0, 6));
+
+  return TriaIterator<dim,MGDoFQuadAccessor<dim,QuadAccessor<dim> > >
+    (
+      tria,
+      present_level,
+      quad_index (i),
+      mg_dof_handler
+    );
+};
+
+
+
+template <int dim, typename BaseClass>
+TriaIterator<dim,MGDoFHexAccessor<dim,BaseClass> >
+MGDoFHexAccessor<dim,BaseClass>::child (const unsigned int i) const {
+  TriaIterator<dim,MGDoFHexAccessor<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;
+};
+
+
+
+template <int dim, typename BaseClass>
+void
+MGDoFHexAccessor<dim,BaseClass>::copy_from (const MGDoFHexAccessor<dim,BaseClass> &a) {
+  DoFHexAccessor::copy_from (a);
+  set_mg_dof_handler (a.mg_dof_handler);
+};
+
+
+
+
 
 /*------------------------- Functions: MGDoFCellAccessor -----------------------*/
 
@@ -440,3 +611,23 @@ template class TriaActiveIterator<2,MGDoFCellAccessor<2> >;
 #endif
 
 
+#if deal_II_dimension == 3
+template class MGDoFLineAccessor<3,LineAccessor<3> >;
+template class MGDoFQuadAccessor<3,QuadAccessor<3> >;
+template class MGDoFHexAccessor<3,HexAccessor<3> >;
+template class MGDoFHexAccessor<3,CellAccessor<3> >;
+template class MGDoFCellAccessor<3>;
+
+template class TriaRawIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFHexAccessor<3,HexAccessor<3> > >;
+template class TriaRawIterator<3,MGDoFCellAccessor<3> >;
+template class TriaIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaIterator<3,MGDoFCellAccessor<3> >;
+template class TriaActiveIterator<3,MGDoFLineAccessor<3,LineAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFQuadAccessor<3,QuadAccessor<3> > >;
+template class TriaActiveIterator<3,MGDoFCellAccessor<3> >;
+#endif
+
+
index 2d6f4ffaf4a2c8c8f584460ea04e8bd1631ccf27..26c8d83933ce6bdaa4b088e59f75679f53b93b2f 100644 (file)
@@ -353,6 +353,94 @@ MGDoFHandler<1>::last_active_quad () const {
   return 0;
 };
 
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::begin_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::begin_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::begin_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::end_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::raw_hex_iterator
+MGDoFHandler<1>::last_raw_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::hex_iterator
+MGDoFHandler<1>::last_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<1>::active_hex_iterator
+MGDoFHandler<1>::last_active_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
 #endif
 
 
@@ -516,6 +604,261 @@ MGDoFHandler<2>::last_active_face (const unsigned int level) const {
   return last_active_line (level);
 };
 
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::begin_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::begin_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::begin_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::end_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex (const unsigned int) const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::raw_hex_iterator
+MGDoFHandler<2>::last_raw_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::hex_iterator
+MGDoFHandler<2>::last_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+template <>
+MGDoFHandler<2>::active_hex_iterator
+MGDoFHandler<2>::last_active_hex () const {
+  Assert (false, ExcNotImplemented());
+  return 0;
+};
+
+
+
+#endif
+
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::begin_raw (const unsigned int level) const {
+  return begin_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::begin (const unsigned int level) const {
+  return begin_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::begin_active (const unsigned int level) const {
+  return begin_active_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::end () const {
+  return end_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw () const {
+  return last_raw_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::raw_cell_iterator
+MGDoFHandler<3>::last_raw (const unsigned int level) const {
+  return last_raw_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last () const {
+  return last_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::cell_iterator
+MGDoFHandler<3>::last (const unsigned int level) const {
+  return last_hex (level);
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active () const {
+  return last_active_hex ();
+};
+
+
+
+template <>
+MGDoFHandler<3>::active_cell_iterator
+MGDoFHandler<3>::last_active (const unsigned int level) const {
+  return last_active_hex (level);
+};
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::begin_raw_face (const unsigned int level) const {
+  return begin_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::begin_face (const unsigned int level) const {
+  return begin_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::begin_active_face (const unsigned int level) const {
+  return begin_active_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::end_face () const {
+  return end_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face () const {
+  return last_raw_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::raw_face_iterator
+MGDoFHandler<3>::last_raw_face (const unsigned int level) const {
+  return last_raw_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face () const {
+  return last_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::face_iterator
+MGDoFHandler<3>::last_face (const unsigned int level) const {
+  return last_quad (level);
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face () const {
+  return last_active_quad ();
+};
+
+
+
+template <>
+MGDoFDimensionInfo<3>::active_face_iterator
+MGDoFHandler<3>::last_active_face (const unsigned int level) const {
+  return last_active_quad (level);
+};
+
 #endif
 
 
@@ -588,6 +931,39 @@ MGDoFHandler<dim>::begin_active_quad (const unsigned int level) const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::begin_raw_hex (const unsigned int level) const {
+  return raw_hex_iterator (tria,
+                          tria->begin_raw_hex(level)->level(),
+                          tria->begin_raw_hex(level)->index(),
+                          this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::begin_hex (const unsigned int level) const {
+  return hex_iterator (tria,
+                      tria->begin_hex(level)->level(),
+                      tria->begin_hex(level)->index(),
+                      this);
+};
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::begin_active_hex (const unsigned int level) const {
+  return active_hex_iterator (tria,
+                             tria->begin_active_hex(level)->level(),
+                             tria->begin_active_hex(level)->index(),
+                             this);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::raw_line_iterator
 MGDoFHandler<dim>::end_line () const {
@@ -604,6 +980,14 @@ MGDoFHandler<dim>::end_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_hex () const {
+  return raw_hex_iterator (tria, -1, -1, this);
+};
+
+
+
 template <int dim>
 typename MGDoFDimensionInfo<dim>::raw_cell_iterator
 MGDoFHandler<dim>::end_raw (const unsigned int level) const {
@@ -694,34 +1078,65 @@ MGDoFHandler<dim>::end_active_line (const unsigned int level) const {
 
 
 
-
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::raw_quad_iterator
+MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         end_quad() :
+         begin_raw_quad (level+1));
+};
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::quad_iterator
+MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         quad_iterator(end_quad()) :
+         begin_quad (level+1));
+};
+
+
+
+
+template <int dim>
+typename MGDoFDimensionInfo<dim>::active_quad_iterator
+MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+  return (level == mg_levels.size()-1 ?
+         active_quad_iterator(end_quad()) :
+         begin_active_quad (level+1));
+};
+
+
+
 template <int dim>
-typename MGDoFDimensionInfo<dim>::raw_quad_iterator
-MGDoFHandler<dim>::end_raw_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::raw_hex_iterator
+MGDoFHandler<dim>::end_raw_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         end_quad() :
-         begin_raw_quad (level+1));
+         end_hex() :
+         begin_raw_hex (level+1));
 };
 
 
 
 template <int dim>
-typename MGDoFDimensionInfo<dim>::quad_iterator
-MGDoFHandler<dim>::end_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::hex_iterator
+MGDoFHandler<dim>::end_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         quad_iterator(end_quad()) :
-         begin_quad (level+1));
+         hex_iterator(end_hex()) :
+         begin_hex (level+1));
 };
 
 
 
 
 template <int dim>
-typename MGDoFDimensionInfo<dim>::active_quad_iterator
-MGDoFHandler<dim>::end_active_quad (const unsigned int level) const {
+typename MGDoFDimensionInfo<dim>::active_hex_iterator
+MGDoFHandler<dim>::end_active_hex (const unsigned int level) const {
   return (level == mg_levels.size()-1 ?
-         active_quad_iterator(end_quad()) :
-         begin_active_quad (level+1));
+         active_hex_iterator(end_hex()) :
+         begin_active_hex (level+1));
 };
 
 
@@ -795,6 +1210,42 @@ MGDoFHandler<dim>::last_active_quad (const unsigned int level) const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex (const unsigned int level) const {
+  return raw_hex_iterator (tria,
+                          tria->last_raw_hex(level)->level(),
+                          tria->last_raw_hex(level)->index(),
+                          this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex (const unsigned int level) const {
+  return hex_iterator (tria,
+                      tria->last_hex(level)->level(),
+                      tria->last_hex(level)->index(),
+                      this);
+};
+
+
+
+
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex (const unsigned int level) const {
+  return active_hex_iterator (tria,
+                             tria->last_active_hex(level)->level(),
+                             tria->last_active_hex(level)->index(),
+                             this);
+};
+
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::raw_line_iterator
 MGDoFHandler<dim>::last_raw_line () const {
@@ -811,6 +1262,14 @@ MGDoFHandler<dim>::last_raw_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::raw_hex_iterator
+MGDoFHandler<dim>::last_raw_hex () const {
+  return last_raw_hex (mg_levels.size()-1);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::line_iterator
 MGDoFHandler<dim>::last_line () const {
@@ -827,6 +1286,14 @@ MGDoFHandler<dim>::last_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::hex_iterator
+MGDoFHandler<dim>::last_hex () const {
+  return last_hex (mg_levels.size()-1);
+};
+
+
+
 template <int dim>
 typename MGDoFHandler<dim>::active_line_iterator
 MGDoFHandler<dim>::last_active_line () const {
@@ -844,6 +1311,14 @@ MGDoFHandler<dim>::last_active_quad () const {
 
 
 
+template <int dim>
+typename MGDoFHandler<dim>::active_hex_iterator
+MGDoFHandler<dim>::last_active_hex () const {
+  return last_active_hex (mg_levels.size()-1);
+};
+
+
+
 
 
 
@@ -989,6 +1464,69 @@ MGDoFHandler<2>::distribute_dofs_on_cell (cell_iterator &cell,
 
 
 
+#if deal_II_dimension == 3
+
+template <>
+unsigned int
+MGDoFHandler<3>::distribute_dofs_on_cell (cell_iterator &cell,
+                                         unsigned int   next_free_dof) {
+  if (selected_fe->dofs_per_vertex > 0)
+                                    // number dofs on vertices
+    for (unsigned int vertex=0; vertex<GeometryInfo<3>::vertices_per_cell; ++vertex)
+                                      // check whether dofs for this
+                                      // vertex have been distributed
+                                      // (only check the first dof)
+      if (cell->mg_vertex_dof_index(vertex, 0) == -1)
+       for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+         cell->set_mg_vertex_dof_index (vertex, d, next_free_dof++);
+    
+                                  // for the lines
+  if (selected_fe->dofs_per_line > 0)
+    for (unsigned int l=0; l<GeometryInfo<3>::lines_per_cell; ++l)
+      {
+       line_iterator line = cell->line(l);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (line->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           line->set_mg_dof_index (d, next_free_dof++);            
+      };
+
+                                  // for the quads
+  if (selected_fe->dofs_per_quad > 0)
+    for (unsigned int q=0; q<GeometryInfo<3>::lines_per_cell; ++q)
+      {
+       quad_iterator quad = cell->quad(q);
+       
+                                        // distribute dofs if necessary:
+                                        // check whether line dof is already
+                                        // numbered (check only first dof)
+       if (quad->mg_dof_index(0) == -1)
+                                          // if not: distribute dofs
+         for (unsigned int d=0; d<selected_fe->dofs_per_line; ++d)
+           quad->set_mg_dof_index (d, next_free_dof++);            
+      };
+  
+
+                                      // dofs of cell
+  if (selected_fe->dofs_per_hex > 0)
+    for (unsigned int d=0; d<selected_fe->dofs_per_hex; ++d)
+      cell->set_mg_dof_index (d, next_free_dof++);
+
+  
+                                  // note that this cell has been processed
+  cell->set_user_flag ();
+  
+  return next_free_dof;
+};
+
+#endif
+
+
+
 
 template <int dim>
 unsigned int MGDoFHandler<dim>::n_dofs (const unsigned int level) const {
@@ -1025,182 +1563,11 @@ void MGDoFHandler<dim>::make_sparsity_pattern (const unsigned int  level,
 
 
 
-
-template <int dim>
-void MGDoFHandler<dim>::renumber_dofs (const unsigned int      level,
-                                      const RenumberingMethod method,
-                                      const vector<int>      &starting_points) {
-                                  // make the connection graph
-  SparseMatrixStruct sparsity (n_dofs(level), max_couplings_between_dofs());
-  make_sparsity_pattern (level, sparsity);
-    
-  int n_dofs = sparsity.n_rows();
-                                  // store the new dof numbers; -1 means
-                                  // that no new number was chosen yet
-  vector<int> new_number(n_dofs, -1);
-  
-                                  // store the indices of the dofs renumbered
-                                  // in the last round. Default to starting
-                                  // points
-  vector<int> last_round_dofs (starting_points);
-  
-                                  // delete disallowed elements
-  for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-    if ((last_round_dofs[i]<0) || (last_round_dofs[i]>=n_dofs))
-      last_round_dofs[i] = -1;
-  
-  remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-            bind2nd(equal_to<int>(), -1));
-  
-                                  // now if no valid points remain:
-                                  // find dof with lowest coordination
-                                  // number
-  
-  if (last_round_dofs.size() == 0)
-    {
-      int          starting_point   = -1;
-      unsigned int min_coordination = n_dofs;
-      for (int row=0; row<n_dofs; ++row) 
-       {
-         unsigned int j;
-         for (j=sparsity.get_rowstart_indices()[row];
-              j<sparsity.get_rowstart_indices()[row+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-                                          // post-condition after loop:
-                                          // coordination is now
-                                          // j-rowstart[row]
-         if (j-sparsity.get_rowstart_indices()[row] <  min_coordination)
-           {
-             min_coordination = j-sparsity.get_rowstart_indices()[row];
-             starting_point   = row;
-           };
-       };
-                                      // initialize the first dof
-      last_round_dofs.push_back (starting_point);
-    };
-  
-
-                                  // store next free dof index
-  int         next_free_number = 0;
-
-                                  // enumerate the first round dofs
-  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-    new_number[last_round_dofs[i]] = next_free_number++;
-
-  bool all_dofs_renumbered = false;
-
-                                  // now do as many steps as needed to
-                                  // renumber all dofs
-  while (!all_dofs_renumbered) 
-    {
-                                      // store the indices of the dofs to be
-                                      // renumbered in the next round
-      vector<int> next_round_dofs;
-
-                                      // find all neighbors of the
-                                      // dofs numbered in the last
-                                      // round
-      for (unsigned int i=0; i<last_round_dofs.size(); ++i)
-       for (unsigned int j=sparsity.get_rowstart_indices()[last_round_dofs[i]];
-            j<sparsity.get_rowstart_indices()[last_round_dofs[i]+1]; ++j)
-         if (sparsity.get_column_numbers()[j] == -1)
-           break;
-         else
-           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                      // sort dof numbers
-      sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                      // delete multiple entries
-      vector<int>::iterator end_sorted;
-      end_sorted = unique (next_round_dofs.begin(), next_round_dofs.end());
-      next_round_dofs.erase (end_sorted, next_round_dofs.end());
-
-                                      // eliminate dofs which are
-                                      // already numbered
-      for (int s=next_round_dofs.size()-1; s>=0; --s)
-       if (new_number[next_round_dofs[s]] != -1)
-         next_round_dofs.erase (&next_round_dofs[s]);
-
-                                      // check whether there are any new
-                                      // dofs
-      all_dofs_renumbered = (next_round_dofs.size() == 0);
-      if (all_dofs_renumbered)
-                                        // end loop if possible
-       continue;
-      
-
-                                      // store for each coordination
-                                      // number the dofs with these
-                                      // coordination number
-      multimap<unsigned int, int> dofs_by_coordination;
-      
-                                      // find coordination number for
-                                      // each of these dofs
-      for (vector<int>::iterator s=next_round_dofs.begin();
-          s!=next_round_dofs.end(); ++s) 
-       {
-         unsigned int coordination = 0;
-         for (unsigned int j=sparsity.get_rowstart_indices()[*s];
-              j<sparsity.get_rowstart_indices()[*s+1]; ++j)
-           if (sparsity.get_column_numbers()[j] == -1)
-             break;
-           else
-             ++coordination;
-
-                                          // insert this dof at its
-                                          // coordination number
-         const pair<const unsigned int, int> new_entry (coordination, *s);
-         dofs_by_coordination.insert (new_entry);
-       };
-      
-                                      ////
-      multimap<unsigned int, int>::iterator i;
-      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-       new_number[i->second] = next_free_number++;
-
-                                      // after that: copy this round's
-                                      // dofs for the next round
-      last_round_dofs = next_round_dofs;
-    };
-
-#ifdef DEBUG
-                                  //  test for all indices numbered
-  if (find (new_number.begin(), new_number.end(), -1) != new_number.end())
-    Assert (false, ExcRenumberingIncomplete());
-  Assert (next_free_number == n_dofs,
-         ExcRenumberingIncomplete());
-#endif
-
-  switch (method) 
-    {
-      case Cuthill_McKee:
-           break;
-      case reverse_Cuthill_McKee:
-      {
-       for (vector<int>::iterator i=new_number.begin(); i!=new_number.end(); ++i)
-         *i = n_dofs-*i;
-       break;
-      };
-      default:
-           Assert (false, ExcNotImplemented());
-    };
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  do_renumbering (level, new_number);
-};
-
-
-
-
 #if deal_II_dimension == 1
 
 template <>
-void MGDoFHandler<1>::do_renumbering (const unsigned int level,
-                                     const vector<int> &new_numbers) {
+void MGDoFHandler<1>::renumber_dofs (const unsigned int level,
+                                    const vector<int> &new_numbers) {
   Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
   
                                   // note that we can not use cell iterators
@@ -1236,8 +1603,45 @@ void MGDoFHandler<1>::do_renumbering (const unsigned int level,
 #if deal_II_dimension == 2
 
 template <>
-void MGDoFHandler<2>::do_renumbering (const unsigned int  level,
-                                     const vector<int>  &new_numbers) {
+void MGDoFHandler<2>::renumber_dofs (const unsigned int  level,
+                                    const vector<int>  &new_numbers) {
+  Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
+  
+  for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
+       i!=mg_vertex_dofs.end(); ++i)
+                                    // if the present vertex lives on
+                                    // the present level
+    if ((i->get_coarsest_level() <= level) &&
+       (i->get_finest_level() >= level))
+      for (unsigned int d=0; d<selected_fe->dofs_per_vertex; ++d)
+       i->set_index (level, d, selected_fe->dofs_per_vertex,
+                     new_numbers[i->get_index (level, d,
+                                               selected_fe->dofs_per_vertex)]);
+  
+  for (vector<int>::iterator i=mg_levels[level]->line_dofs.begin();
+       i!=mg_levels[level]->line_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
+
+  for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
+       i!=mg_levels[level]->quad_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
+};
+
+#endif
+
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::renumber_dofs (const unsigned int  level,
+                                    const vector<int>  &new_numbers) {
   Assert (new_numbers.size() == n_dofs(level), ExcRenumberingIncomplete());
   
   for (vector<MGVertexDoFs>::iterator i=mg_vertex_dofs.begin();
@@ -1257,12 +1661,20 @@ void MGDoFHandler<2>::do_renumbering (const unsigned int  level,
       Assert (*i != -1, ExcInternalError());
       *i = new_numbers[*i];
     };
+
   for (vector<int>::iterator i=mg_levels[level]->quad_dofs.begin();
        i!=mg_levels[level]->quad_dofs.end(); ++i)
     {
       Assert (*i != -1, ExcInternalError());
       *i = new_numbers[*i];
     };
+
+  for (vector<int>::iterator i=mg_levels[level]->hex_dofs.begin();
+       i!=mg_levels[level]->hex_dofs.end(); ++i)
+    {
+      Assert (*i != -1, ExcInternalError());
+      *i = new_numbers[*i];
+    };
 };
 
 #endif
@@ -1445,5 +1857,100 @@ void MGDoFHandler<2>::reserve_space () {
 #endif
 
 
+
+
+#if deal_II_dimension == 3
+
+template <>
+void MGDoFHandler<3>::reserve_space () {
+  const unsigned int dim = 3;
+  
+  Assert (selected_fe != 0, ExcNoFESelected());
+  Assert (tria->n_levels() > 0, ExcInvalidTriangulation());
+  
+                                  ////////////////////////////
+                                  // DESTRUCTION
+
+                                   // delete all levels and set them up
+                                   // newly, since vectors are
+                                   // troublesome if you want to change
+                                   // their size
+  for (unsigned int i=0; i<mg_levels.size(); ++i)
+    delete mg_levels[i];
+  mg_levels.clear ();
+
+                                  // also delete vector of vertex indices
+                                  // this calls the destructor which
+                                  // must free the space
+  mg_vertex_dofs.clear ();
+  
+  
+                                  ////////////////////////////
+                                  // CONSTRUCTION
+  
+                                  // first allocate space for the
+                                  // lines and quads on each level
+  for (unsigned int i=0; i<tria->n_levels(); ++i) 
+    {
+      mg_levels.push_back (new DoFLevel<3>);
+
+      mg_levels.back()->line_dofs = vector<int> (tria->levels[i]->lines.lines.size() *
+                                                selected_fe->dofs_per_line,
+                                                -1);
+      mg_levels.back()->quad_dofs = vector<int> (tria->levels[i]->quads.quads.size() *
+                                                selected_fe->dofs_per_quad,
+                                                -1);
+      mg_levels.back()->hex_dofs = vector<int> (tria->levels[i]->hexes.hexes.size() *
+                                               selected_fe->dofs_per_hex,
+                                               -1);
+    };
+
+  
+                                  // now allocate space for the
+                                  // vertices. To this end, we need
+                                  // to construct as many objects as
+                                  // there are vertices and let them
+                                  // allocate enough space for their
+                                  // vertex indices on the levels they
+                                  // live on. We need therefore to
+                                  // count to how many levels a cell
+                                  // belongs to, which we do by looping
+                                  // over all cells and storing the
+                                  // maximum and minimum level each
+                                  // vertex we pass by  belongs to
+  mg_vertex_dofs.resize (tria->vertices.size());
+
+  vector<unsigned int> min_level (tria->vertices.size(), tria->n_levels());
+  vector<unsigned int> max_level (tria->vertices.size(), 0);
+
+  Triangulation<dim>::cell_iterator cell = tria->begin(),
+                                   endc = tria->end();
+  for (; cell!=endc; ++cell)
+    for (unsigned int vertex=0; vertex<GeometryInfo<dim>::vertices_per_cell;
+        ++vertex)
+      {
+       const unsigned int vertex_index = cell->vertex_index(vertex);
+       if (min_level[vertex_index] > static_cast<unsigned int>(cell->level()))
+         min_level[vertex_index] = cell->level();
+       if (max_level[vertex_index] < static_cast<unsigned int>(cell->level()))
+         max_level[vertex_index] = cell->level();
+      };
+
+  
+                                  // now allocate the needed space
+  for (unsigned int vertex=0; vertex<tria->vertices.size(); ++vertex)
+    {
+      Assert (min_level[vertex] < tria->n_levels(),   ExcInternalError());
+      Assert (max_level[vertex] >= min_level[vertex], ExcInternalError());
+
+      mg_vertex_dofs[vertex].init (min_level[vertex],
+                                  max_level[vertex],
+                                  selected_fe->dofs_per_vertex);
+    };
+};
+
+#endif
+
+
 // explicite instantiations
 template class MGDoFHandler<deal_II_dimension>;
index cc8a9f83f524700cd13f3d92c0ce50030fccc51b..8ca5178abfabad62e4f98de2020bfe16b2a5e373 100644 (file)
@@ -12,7 +12,7 @@
 #include <lac/sparsematrix.h>
 #include <base/parameter_handler.h>
 #include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
 
 #include <fstream>
 #include <cmath>
@@ -354,7 +354,7 @@ void TestCases<dim>::run (ParameterHandler &prm) {
   dof->distribute_dofs (fe);
 
   cout << "    Renumbering degrees of freedom..." << endl;
-  dof->renumber_dofs (Cuthill_McKee, false);
+  DoFRenumbering::renumber_Cuthill_McKee (*dof);
     
   SparseMatrixStruct sparsity (dof->n_dofs(),
                               dof->max_couplings_between_dofs());
index f73b99509bf8051741677e9b48782f56018cfb5c..e5f74da24382dd6a46235b38ce2617ce52cf5f1e 100644 (file)
@@ -14,7 +14,7 @@
 #include <grid/tria_accessor.h>
 #include <lac/sparsematrix.h>
 #include <grid/dof_constraints.h>
-
+#include <numerics/dof_renumbering.h>
 
 #include <cmath>
 #include <cstdlib>
@@ -287,7 +287,7 @@ void TestCases<dim>::run (const unsigned int test_case) {
   dof->distribute_dofs (fe);
 
   cout << "    Renumbering degrees of freedom..." << endl;
-  dof->renumber_dofs (Cuthill_McKee, false);
+  DoFRenumbering::renumber_Cuthill_McKee (*dof);
     
   SparseMatrixStruct sparsity (dof->n_dofs(),
                               dof->max_couplings_between_dofs());

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.