]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Convert DoFRenumbering from a class with all static members to a namespace.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Mar 2008 22:11:09 +0000 (22:11 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Fri, 7 Mar 2008 22:11:09 +0000 (22:11 +0000)
git-svn-id: https://svn.dealii.org/trunk@15878 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_renumbering.h
deal.II/deal.II/source/dofs/dof_renumbering.cc
deal.II/doc/news/changes.h

index 632ca14b95d1dc9575e8962b4995184efee02e05..a85a2bd14453bc97d6c4c345c8fd94255494a961 100644 (file)
@@ -367,169 +367,33 @@ DEAL_II_NAMESPACE_OPEN
  * @ingroup dofs
  * @author Wolfgang Bangerth, Guido Kanschat, 1998, 1999, 2000, 2004, 2007, 2008
  */
-class DoFRenumbering 
+namespace DoFRenumbering 
 {
-  public:
-                                    /**
-                                     * A namespace for the
-                                     * implementation of some
-                                     * renumbering algorithms based
-                                     * on algorithms implemented in
-                                     * the Boost Graph Library (BGL)
-                                     * by Jeremy Siek and others.
-                                     *
-                                     * While often slighty slower to
-                                     * compute, the algorithms using
-                                     * BOOST often lead to matrices
-                                     * with smaller bandwidths and
-                                     * sparse ILUs based on this
-                                     * numbering are therefore more
-                                     * efficient.
-                                     *
-                                     * For a comparison of these
-                                     * algorithms with the ones
-                                     * defined in DoFRenumbering, see
-                                     * the comparison section in the
-                                     * documentation of the
-                                     * DoFRenumbering namespace.
-                                     */
-    class boost
-    {
-      public:
-                                        /**
-                                         * Renumber the degrees of
-                                         * freedom according to the
-                                         * Cuthill-McKee method,
-                                         * eventually using the reverse
-                                         * numbering scheme.
-                                         *
-                                         * See the general
-                                         * documentation of the
-                                         * parent class for details
-                                         * on the different methods.
-                                         *
-                                         * As an example of the
-                                         * results of this algorithm,
-                                         * take a look at the
-                                         * comparison of various
-                                         * algorithms in the
-                                         * documentation of the
-                                         * DoFRenumbering namespace.
-                                         */
-       template <class DH>
-       static void 
-       Cuthill_McKee (DH&                              dof_handler,
-                      const bool                       reversed_numbering = false,
-                      const bool                       use_constraints    = false);
-
-                                        /**
-                                         * Computes the renumbering
-                                         * vector needed by the
-                                         * Cuthill_McKee() function. Does
-                                         * not perform the renumbering on
-                                         * the DoFHandler dofs but
-                                         * returns the renumbering
-                                         * vector.
-                                         */    
-       template <class DH>
-       static void
-       compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
-                              const DH&,
-                              const bool reversed_numbering = false,
-                              const bool use_constraints    = false);
-
-                                        /**
-                                         * Renumber the degrees of
-                                         * freedom based on the BOOST
-                                         * implementation of the King
-                                         * algorithm. This often
-                                         * results in slightly larger
-                                         * (by a few percent)
-                                         * bandwidths than the
-                                         * Cuthill-McKee algorithm,
-                                         * but sparse ILUs are often
-                                         * slightly (also by a few
-                                         * percent) better
-                                         * preconditioners.
-                                         *
-                                         * As an example of the
-                                         * results of this algorithm,
-                                         * take a look at the
-                                         * comparison of various
-                                         * algorithms in the
-                                         * documentation of the
-                                         * DoFRenumbering namespace.
-                                         *
-                                         * This algorithm is used in
-                                         * @ref step_22 "step-22".
-                                         */
-       template <class DH>
-       static void 
-       king_ordering (DH&                              dof_handler,
-                      const bool                       reversed_numbering = false,
-                      const bool                       use_constraints    = false);
-
-                                        /**
-                                         * Compute the renumbering
-                                         * for the King algorithm but
-                                         * do not actually renumber
-                                         * the degrees of freedom in
-                                         * the DoF handler argument.
-                                         */
-       template <class DH>
-       static void
-       compute_king_ordering (std::vector<unsigned int>& new_dof_indices,
-                              const DH&,
-                              const bool reversed_numbering = false,
-                              const bool use_constraints    = false);
-
-                                        /**
-                                         * Renumber the degrees of
-                                         * freedom based on the BOOST
-                                         * implementation of the
-                                         * minimum degree
-                                         * algorithm. Unlike the
-                                         * Cuthill-McKee algorithm,
-                                         * this algorithm does not
-                                         * attempt to minimize the
-                                         * bandwidth of a matrix but
-                                         * to minimize the amount of
-                                         * fill-in when doing an LU
-                                         * decomposition. It may
-                                         * sometimes yield better
-                                         * ILUs because of this
-                                         * property.
-                                         *
-                                         * As an example of the
-                                         * results of this algorithm,
-                                         * take a look at the
-                                         * comparison of various
-                                         * algorithms in the
-                                         * documentation of the
-                                         * DoFRenumbering namespace.
-                                         */
-       template <class DH>
-       static void 
-       minimum_degree (DH&                              dof_handler,
-                       const bool                       reversed_numbering = false,
-                       const bool                       use_constraints    = false);
-
-                                        /**
-                                         * Compute the renumbering
-                                         * for the minimum degree
-                                         * algorithm but do not
-                                         * actually renumber the
-                                         * degrees of freedom in the
-                                         * DoF handler argument.
-                                         */
-       template <class DH>
-       static void
-       compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
-                               const DH&,
-                               const bool reversed_numbering = false,
-                               const bool use_constraints    = false);
-    };
-    
+                                  /**
+                                   * A namespace for the
+                                   * implementation of some
+                                   * renumbering algorithms based
+                                   * on algorithms implemented in
+                                   * the Boost Graph Library (BGL)
+                                   * by Jeremy Siek and others.
+                                   *
+                                   * While often slighty slower to
+                                   * compute, the algorithms using
+                                   * BOOST often lead to matrices
+                                   * with smaller bandwidths and
+                                   * sparse ILUs based on this
+                                   * numbering are therefore more
+                                   * efficient.
+                                   *
+                                   * For a comparison of these
+                                   * algorithms with the ones
+                                   * defined in DoFRenumbering, see
+                                   * the comparison section in the
+                                   * documentation of the
+                                   * DoFRenumbering namespace.
+                                   */
+  namespace boost
+  {
                                     /**
                                      * Renumber the degrees of
                                      * freedom according to the
@@ -537,23 +401,24 @@ class DoFRenumbering
                                      * eventually using the reverse
                                      * numbering scheme.
                                      *
-                                     * See the general documentation
-                                     * of this class for details on
-                                     * the different methods.
+                                     * See the general
+                                     * documentation of the
+                                     * parent class for details
+                                     * on the different methods.
                                      *
-                                     * As an example of the results
-                                     * of this algorithm, take a look
-                                     * at the comparison of various
+                                     * As an example of the
+                                     * results of this algorithm,
+                                     * take a look at the
+                                     * comparison of various
                                      * algorithms in the
                                      * documentation of the
                                      * DoFRenumbering namespace.
                                      */
     template <class DH>
-    static void 
+    void 
     Cuthill_McKee (DH&                              dof_handler,
                   const bool                       reversed_numbering = false,
-                  const bool                       use_constraints    = false,
-                  const std::vector<unsigned int> &starting_indices   = std::vector<unsigned int>());
+                  const bool                       use_constraints    = false);
 
                                     /**
                                      * Computes the renumbering
@@ -565,598 +430,730 @@ class DoFRenumbering
                                      * vector.
                                      */    
     template <class DH>
-    static void
+    void
     compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
                           const DH&,
                           const bool reversed_numbering = false,
-                          const bool use_constraints    = false,
-                          const std::vector<unsigned int> &starting_indices   = std::vector<unsigned int>());
+                          const bool use_constraints    = false);
 
                                     /**
                                      * 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.
+                                     * freedom based on the BOOST
+                                     * implementation of the King
+                                     * algorithm. This often
+                                     * results in slightly larger
+                                     * (by a few percent)
+                                     * bandwidths than the
+                                     * Cuthill-McKee algorithm,
+                                     * but sparse ILUs are often
+                                     * slightly (also by a few
+                                     * percent) better
+                                     * preconditioners.
                                      *
-                                     * See the general documentation
-                                     * of this class for details on
-                                     * the different methods.
-                                     */
-    template <int dim>
-    static void 
-    Cuthill_McKee (MGDoFHandler<dim>          &dof_handler,
-                  const unsigned int          level,
-                  const bool                  reversed_numbering = false,
-                  const std::vector<unsigned int> &starting_indices   = std::vector<unsigned int> ());
-
-                                    /**
-                                     * Sort the degrees of freedom by
-                                     * vector component. The
-                                     * numbering within each
-                                     * component is not touched, so a
-                                     * degree of freedom with index
-                                     * $i$, belonging to some
-                                     * component, and another degree
-                                     * of freedom with index $j$
-                                     * belonging to the same
-                                     * component will be assigned new
-                                     * indices $n(i)$ and $n(j)$ with
-                                     * $n(i)<n(j)$ if $i<j$ and
-                                     * $n(i)>n(j)$ if $i>j$.
-                                     *
-                                     * You can specify that the
-                                     * components are ordered in a
-                                     * different way than suggested
-                                     * by the FESystem object you
-                                     * use. To this end, set up the
-                                     * vector @p target_component
-                                     * such that the entry at index
-                                     * @p i denotes the number of the
-                                     * target component for dofs with
-                                     * component @p i in the
-                                     * FESystem. Naming the same
-                                     * component more than once is
-                                     * possible and results in a
-                                     * blocking of several components
-                                     * into one. This is discussed in
-                                     * @ref step_22 "step-22". If you
-                                     * omit this argument, the same
-                                     * order as given by the finite
-                                     * element is used.
+                                     * As an example of the
+                                     * results of this algorithm,
+                                     * take a look at the
+                                     * comparison of various
+                                     * algorithms in the
+                                     * documentation of the
+                                     * DoFRenumbering namespace.
                                      *
-                                     * If one of the base finite
-                                     * elements from which the global
-                                     * finite element under
-                                     * consideration here, is a
-                                     * non-primitive one, i.e. its
-                                     * shape functions have more than
-                                     * one non-zero component, then
-                                     * it is not possible to
-                                     * associate these degrees of
-                                     * freedom with a single vector
-                                     * component. In this case, they
-                                     * are associated with the first
-                                     * vector component to which they
-                                     * belong.
-                                     * 
-                                     * For finite elements with only
-                                     * one component, or a single
-                                     * non-primitive base element,
-                                     * this function is the identity
-                                     * operation.
-                                     */
-    template <int dim>
-    static void
-    component_wise (DoFHandler<dim>                 &dof_handler,
-                   const std::vector<unsigned int> &target_component
-                   = std::vector<unsigned int>());
-
-                                    /**
-                                     * Sort the degrees of freedom by
-                                     * component. It does the same
-                                     * thing as the above function,
-                                     * only that it does this for one
-                                     * single level of a multi-level
-                                     * discretization. The
-                                     * non-multigrid part of the
-                                     * MGDoFHandler is not touched.
+                                     * This algorithm is used in
+                                     * @ref step_22 "step-22".
                                      */
-    template <int dim>
-    static void
-    component_wise (MGDoFHandler<dim>&               dof_handler,
-                   unsigned int                     level,
-                   const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
-
+    template <class DH>
+    void 
+    king_ordering (DH&                              dof_handler,
+                  const bool                       reversed_numbering = false,
+                  const bool                       use_constraints    = false);
 
                                     /**
-                                     * Sort the degrees of freedom by
-                                     * component. It does the same
-                                     * thing as the previous
-                                     * functions, but more: it
-                                     * renumbers not only every level
-                                     * of the multigrid part, but
-                                     * also the global,
-                                     * i.e. non-multigrid components.
-                                     */
-    template <int dim>
-    static void
-    component_wise (MGDoFHandler<dim>&               dof_handler,
-                   const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
-
-    
-  protected:
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * component_wise()
-                                     * functions. Does not perform
-                                     * the renumbering on the
-                                     * DoFHandler dofs but returns
-                                     * the renumbering vector.
-                                     */    
-    template <int dim, class ITERATOR, class ENDITERATOR>
-    static unsigned int
-    compute_component_wise (std::vector<unsigned int>& new_dof_indices,
-                           ITERATOR& start,
-                           const ENDITERATOR& end,
-                           const std::vector<unsigned int> &target_component);
-  public:
-                                    /**
-                                     * Cell-wise renumbering for DG
-                                     * elements.  This function takes
-                                     * the ordered set of cells in
-                                     * <tt>cell_order</tt>, and makes
-                                     * sure that all degrees of
-                                     * freedom in a cell with higher
-                                     * index are behind all degrees
-                                     * of freedom of a cell with
-                                     * lower index. The order inside
-                                     * a cell bloock will be the same
-                                     * as before this renumbering.
-                                     *
-                                     * This function only works with
-                                     * Discontinuous Galerkin Finite
-                                     * Elements, i.e. all degrees of
-                                     * freedom have to be associated
-                                     * with the interior of the cell.
+                                     * Compute the renumbering
+                                     * for the King algorithm but
+                                     * do not actually renumber
+                                     * the degrees of freedom in
+                                     * the DoF handler argument.
                                      */
     template <class DH>
-    static void
-    cell_wise (DH &dof_handler,
-              const std::vector<typename DH::cell_iterator> &cell_order);
+    void
+    compute_king_ordering (std::vector<unsigned int>& new_dof_indices,
+                          const DH&,
+                          const bool reversed_numbering = false,
+                          const bool use_constraints    = false);
 
                                     /**
-                                     * @deprecated Use cell_wise() instead.
+                                     * Renumber the degrees of
+                                     * freedom based on the BOOST
+                                     * implementation of the
+                                     * minimum degree
+                                     * algorithm. Unlike the
+                                     * Cuthill-McKee algorithm,
+                                     * this algorithm does not
+                                     * attempt to minimize the
+                                     * bandwidth of a matrix but
+                                     * to minimize the amount of
+                                     * fill-in when doing an LU
+                                     * decomposition. It may
+                                     * sometimes yield better
+                                     * ILUs because of this
+                                     * property.
                                      *
-                                     * Cell-wise numbering only for DG.
+                                     * As an example of the
+                                     * results of this algorithm,
+                                     * take a look at the
+                                     * comparison of various
+                                     * algorithms in the
+                                     * documentation of the
+                                     * DoFRenumbering namespace.
                                      */
     template <class DH>
-    static void
-    cell_wise_dg (DH &dof_handler,
-                 const std::vector<typename DH::cell_iterator> &cell_order);
+    void 
+    minimum_degree (DH&                              dof_handler,
+                   const bool                       reversed_numbering = false,
+                   const bool                       use_constraints    = false);
 
                                     /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * cell_wise_dg() function. Does
-                                     * not perform the renumbering on
-                                     * the DoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
+                                     * Compute the renumbering
+                                     * for the minimum degree
+                                     * algorithm but do not
+                                     * actually renumber the
+                                     * degrees of freedom in the
+                                     * DoF handler argument.
                                      */
     template <class DH>
-    static void
-    compute_cell_wise_dg (std::vector<unsigned int>& renumbering,
-                         std::vector<unsigned int>& inverse_renumbering,
-                         const DH &dof_handler,
-                         const std::vector<typename DH::cell_iterator> &cell_order);
+    void
+    compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
+                           const DH&,
+                           const bool reversed_numbering = false,
+                           const bool use_constraints    = false);
+  }
+    
+                                  /**
+                                   * 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.
+                                   *
+                                   * As an example of the results
+                                   * of this algorithm, take a look
+                                   * at the comparison of various
+                                   * algorithms in the
+                                   * documentation of the
+                                   * DoFRenumbering namespace.
+                                   */
+  template <class DH>
+  void 
+  Cuthill_McKee (DH&                              dof_handler,
+                const bool                       reversed_numbering = false,
+                const bool                       use_constraints    = false,
+                const std::vector<unsigned int> &starting_indices   = std::vector<unsigned int>());
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * cell_wise() function. Does
-                                     * not perform the renumbering on
-                                     * the DoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <class DH>
-    static void
-    compute_cell_wise (std::vector<unsigned int>& renumbering,
-                      std::vector<unsigned int>& inverse_renumbering,
-                      const DH &dof_handler,
-                      const std::vector<typename DH::cell_iterator> &cell_order);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * Cuthill_McKee() function. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */    
+  template <class DH>
+  void
+  compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
+                        const DH&,
+                        const bool reversed_numbering = false,
+                        const bool use_constraints    = false,
+                        const std::vector<unsigned int> &starting_indices   = std::vector<unsigned 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>
+  void 
+  Cuthill_McKee (MGDoFHandler<dim>          &dof_handler,
+                const unsigned int          level,
+                const bool                  reversed_numbering = false,
+                const std::vector<unsigned int> &starting_indices   = std::vector<unsigned int> ());
+
+                                  /**
+                                   * Sort the degrees of freedom by
+                                   * vector component. The
+                                   * numbering within each
+                                   * component is not touched, so a
+                                   * degree of freedom with index
+                                   * $i$, belonging to some
+                                   * component, and another degree
+                                   * of freedom with index $j$
+                                   * belonging to the same
+                                   * component will be assigned new
+                                   * indices $n(i)$ and $n(j)$ with
+                                   * $n(i)<n(j)$ if $i<j$ and
+                                   * $n(i)>n(j)$ if $i>j$.
+                                   *
+                                   * You can specify that the
+                                   * components are ordered in a
+                                   * different way than suggested
+                                   * by the FESystem object you
+                                   * use. To this end, set up the
+                                   * vector @p target_component
+                                   * such that the entry at index
+                                   * @p i denotes the number of the
+                                   * target component for dofs with
+                                   * component @p i in the
+                                   * FESystem. Naming the same
+                                   * component more than once is
+                                   * possible and results in a
+                                   * blocking of several components
+                                   * into one. This is discussed in
+                                   * @ref step_22 "step-22". If you
+                                   * omit this argument, the same
+                                   * order as given by the finite
+                                   * element is used.
+                                   *
+                                   * If one of the base finite
+                                   * elements from which the global
+                                   * finite element under
+                                   * consideration here, is a
+                                   * non-primitive one, i.e. its
+                                   * shape functions have more than
+                                   * one non-zero component, then
+                                   * it is not possible to
+                                   * associate these degrees of
+                                   * freedom with a single vector
+                                   * component. In this case, they
+                                   * are associated with the first
+                                   * vector component to which they
+                                   * belong.
+                                   * 
+                                   * For finite elements with only
+                                   * one component, or a single
+                                   * non-primitive base element,
+                                   * this function is the identity
+                                   * operation.
+                                   */
+  template <int dim>
+  void
+  component_wise (DoFHandler<dim>                 &dof_handler,
+                 const std::vector<unsigned int> &target_component
+                 = std::vector<unsigned int>());
+
+                                  /**
+                                   * Sort the degrees of freedom by
+                                   * component. It does the same
+                                   * thing as the above function,
+                                   * only that it does this for one
+                                   * single level of a multi-level
+                                   * discretization. The
+                                   * non-multigrid part of the
+                                   * MGDoFHandler is not touched.
+                                   */
+  template <int dim>
+  void
+  component_wise (MGDoFHandler<dim>&               dof_handler,
+                 unsigned int                     level,
+                 const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
+
+
+                                  /**
+                                   * Sort the degrees of freedom by
+                                   * component. It does the same
+                                   * thing as the previous
+                                   * functions, but more: it
+                                   * renumbers not only every level
+                                   * of the multigrid part, but
+                                   * also the global,
+                                   * i.e. non-multigrid components.
+                                   */
+  template <int dim>
+  void
+  component_wise (MGDoFHandler<dim>&               dof_handler,
+                 const std::vector<unsigned int>& target_component = std::vector<unsigned int>());
 
-                                    /**
-                                     * Cell-wise renumbering on one
-                                     * level. See the other function
-                                     * with the same name.
-                                     */
-    template <int dim>
-    static void
-    cell_wise (MGDoFHandler<dim>   &dof_handler,
-              const unsigned int   level,
-              const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
     
-                                    /**
-                                     * @deprecated Use cell_wise() instead.
-                                     *
-                                     * Cell-wise renumbering on one
-                                     * level for DG elements. See the
-                                     * other function with the same
-                                     * name.
-                                     */
-    template <int dim>
-    static void
-    cell_wise_dg (MGDoFHandler<dim>   &dof_handler,
-                 const unsigned int   level,
-                 const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * component_wise()
+                                   * functions. Does not perform
+                                   * the renumbering on the
+                                   * DoFHandler dofs but returns
+                                   * the renumbering vector.
+                                   */    
+  template <int dim, class ITERATOR, class ENDITERATOR>
+  static unsigned int
+  compute_component_wise (std::vector<unsigned int>& new_dof_indices,
+                         ITERATOR& start,
+                         const ENDITERATOR& end,
+                         const std::vector<unsigned int> &target_component);
+
+                                  /**
+                                   * Cell-wise renumbering for DG
+                                   * elements.  This function takes
+                                   * the ordered set of cells in
+                                   * <tt>cell_order</tt>, and makes
+                                   * sure that all degrees of
+                                   * freedom in a cell with higher
+                                   * index are behind all degrees
+                                   * of freedom of a cell with
+                                   * lower index. The order inside
+                                   * a cell bloock will be the same
+                                   * as before this renumbering.
+                                   *
+                                   * This function only works with
+                                   * Discontinuous Galerkin Finite
+                                   * Elements, i.e. all degrees of
+                                   * freedom have to be associated
+                                   * with the interior of the cell.
+                                   */
+  template <class DH>
+  void
+  cell_wise (DH &dof_handler,
+            const std::vector<typename DH::cell_iterator> &cell_order);
+
+                                  /**
+                                   * @deprecated Use cell_wise() instead.
+                                   *
+                                   * Cell-wise numbering only for DG.
+                                   */
+  template <class DH>
+  void
+  cell_wise_dg (DH &dof_handler,
+               const std::vector<typename DH::cell_iterator> &cell_order);
+
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * cell_wise_dg() function. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <class DH>
+  void
+  compute_cell_wise_dg (std::vector<unsigned int>& renumbering,
+                       std::vector<unsigned int>& inverse_renumbering,
+                       const DH &dof_handler,
+                       const std::vector<typename DH::cell_iterator> &cell_order);
+
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * cell_wise() function. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <class DH>
+  void
+  compute_cell_wise (std::vector<unsigned int>& renumbering,
+                    std::vector<unsigned int>& inverse_renumbering,
+                    const DH &dof_handler,
+                    const std::vector<typename DH::cell_iterator> &cell_order);
+
+                                  /**
+                                   * Cell-wise renumbering on one
+                                   * level. See the other function
+                                   * with the same name.
+                                   */
+  template <int dim>
+  void
+  cell_wise (MGDoFHandler<dim>   &dof_handler,
+            const unsigned int   level,
+            const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
     
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * cell_wise_dg() level renumbering function. Does
-                                     * not perform the renumbering on
-                                     * the MGDoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <int dim>
-    static void
-    compute_cell_wise_dg (std::vector<unsigned int>& renumbering,
-                         std::vector<unsigned int>& inverse_renumbering,
-                         const MGDoFHandler<dim>&   dof_handler,
-                         const unsigned int         level,
-                         const std::vector<typename MGDoFHandler<dim>::cell_iterator>& cell_order);
+                                  /**
+                                   * @deprecated Use cell_wise() instead.
+                                   *
+                                   * Cell-wise renumbering on one
+                                   * level for DG elements. See the
+                                   * other function with the same
+                                   * name.
+                                   */
+  template <int dim>
+  void
+  cell_wise_dg (MGDoFHandler<dim>   &dof_handler,
+               const unsigned int   level,
+               const std::vector<typename MGDoFHandler<dim>::cell_iterator> &cell_order);
+    
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * cell_wise_dg() level renumbering function. Does
+                                   * not perform the renumbering on
+                                   * the MGDoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <int dim>
+  void
+  compute_cell_wise_dg (std::vector<unsigned int>& renumbering,
+                       std::vector<unsigned int>& inverse_renumbering,
+                       const MGDoFHandler<dim>&   dof_handler,
+                       const unsigned int         level,
+                       const std::vector<typename MGDoFHandler<dim>::cell_iterator>& cell_order);
 
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * cell_wise() level renumbering function. Does
-                                     * not perform the renumbering on
-                                     * the MGDoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <int dim>
-    static void
-    compute_cell_wise (std::vector<unsigned int>& renumbering,
-                      std::vector<unsigned int>& inverse_renumbering,
-                      const MGDoFHandler<dim>&   dof_handler,
-                      const unsigned int         level,
-                      const std::vector<typename MGDoFHandler<dim>::cell_iterator>& cell_order);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * cell_wise() level renumbering function. Does
+                                   * not perform the renumbering on
+                                   * the MGDoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <int dim>
+  void
+  compute_cell_wise (std::vector<unsigned int>& renumbering,
+                    std::vector<unsigned int>& inverse_renumbering,
+                    const MGDoFHandler<dim>&   dof_handler,
+                    const unsigned int         level,
+                    const std::vector<typename MGDoFHandler<dim>::cell_iterator>& cell_order);
 
 
-                                    /**
-                                     * Cell-wise downstream numbering
-                                     * with respect to a constant
-                                     * flow direction.
-                                     *
-                                     * The cells are sorted such that
-                                     * the centers of higher numbers
-                                     * are further downstream with
-                                     * respect to the constant vector
-                                     * @p direction than the centers
-                                     * of lower numbers. Even if this
-                                     * yields a downstream numbering
-                                     * with respect to the flux on
-                                     * the edges for fairly general
-                                     * grids, this might not be
-                                     * guaranteed for all meshes.
-                                     *
-                                     * This function produces a
-                                     * downstream ordering of the
-                                     * mesh cells and calls
-                                     * cell_wise_dg().  Therefore, it
-                                     * only works with Discontinuous
-                                     * Galerkin Finite Elements,
-                                     * i.e. all degrees of freedom
-                                     * have to be associated with the
-                                     * interior of the cell.
-                                     */
-    template <class DH, int dim>
-    static void
-    downstream (DH&               dof_handler,
-               const Point<dim>& direction);
+                                  /**
+                                   * Cell-wise downstream numbering
+                                   * with respect to a constant
+                                   * flow direction.
+                                   *
+                                   * The cells are sorted such that
+                                   * the centers of higher numbers
+                                   * are further downstream with
+                                   * respect to the constant vector
+                                   * @p direction than the centers
+                                   * of lower numbers. Even if this
+                                   * yields a downstream numbering
+                                   * with respect to the flux on
+                                   * the edges for fairly general
+                                   * grids, this might not be
+                                   * guaranteed for all meshes.
+                                   *
+                                   * This function produces a
+                                   * downstream ordering of the
+                                   * mesh cells and calls
+                                   * cell_wise_dg().  Therefore, it
+                                   * only works with Discontinuous
+                                   * Galerkin Finite Elements,
+                                   * i.e. all degrees of freedom
+                                   * have to be associated with the
+                                   * interior of the cell.
+                                   */
+  template <class DH, int dim>
+  void
+  downstream (DH&               dof_handler,
+             const Point<dim>& direction);
 
-                                    /**
-                                     * @deprecated Use downstream() instead.
-                                     */
-    template <class DH, int dim>
-    static void
-    downstream_dg (DH&               dof_handler,
-                  const Point<dim>& direction);
+                                  /**
+                                   * @deprecated Use downstream() instead.
+                                   */
+  template <class DH, int dim>
+  void
+  downstream_dg (DH&               dof_handler,
+                const Point<dim>& direction);
 
-                                    /**
-                                     * Cell-wise downstream numbering
-                                     * with respect to a constant
-                                     * flow direction on one
-                                     * level. See the other function
-                                     * with the same name.
-                                     */
-    template <int dim>
-    static void
-    downstream (MGDoFHandler<dim>  &dof_handler,
-               const unsigned int level,
-               const Point<dim> &direction);
-                                    /**
-                                     * @deprecated Use downstream()
-                                     * instead.
-                                     */
-    template <int dim>
-    static void
-    downstream_dg (MGDoFHandler<dim>  &dof_handler,
-                  const unsigned int level,
-                  const Point<dim> &direction);
+                                  /**
+                                   * Cell-wise downstream numbering
+                                   * with respect to a constant
+                                   * flow direction on one
+                                   * level. See the other function
+                                   * with the same name.
+                                   */
+  template <int dim>
+  void
+  downstream (MGDoFHandler<dim>  &dof_handler,
+             const unsigned int level,
+             const Point<dim> &direction);
+                                  /**
+                                   * @deprecated Use downstream()
+                                   * instead.
+                                   */
+  template <int dim>
+  void
+  downstream_dg (MGDoFHandler<dim>  &dof_handler,
+                const unsigned int level,
+                const Point<dim> &direction);
 
-                                    /**
-                                     * @deprecated The new function
-                                     * of this name computes the
-                                     * renumbering and its inverse at
-                                     * the same time. So, at least if
-                                     * you need both, you should use
-                                     * the other one.
-                                     *
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * downstream_dg() function. Does
-                                     * not perform the renumbering on
-                                     * the DoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <class DH, int dim>
-    static void
-    compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
-                          const DH&                  dof_handler,
-                          const Point<dim>&          direction);
+                                  /**
+                                   * @deprecated The new function
+                                   * of this name computes the
+                                   * renumbering and its inverse at
+                                   * the same time. So, at least if
+                                   * you need both, you should use
+                                   * the other one.
+                                   *
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * downstream_dg() function. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <class DH, int dim>
+  void
+  compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
+                        const DH&                  dof_handler,
+                        const Point<dim>&          direction);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * downstream_dg() function. Does
-                                     * not perform the renumbering on
-                                     * the DoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <class DH, int dim>
-    static void
-    compute_downstream (std::vector<unsigned int>& new_dof_indices,
-                       std::vector<unsigned int>& reverse,
-                       const DH&                  dof_handler,
-                       const Point<dim>&          direction);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * downstream_dg() function. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <class DH, int dim>
+  void
+  compute_downstream (std::vector<unsigned int>& new_dof_indices,
+                     std::vector<unsigned int>& reverse,
+                     const DH&                  dof_handler,
+                     const Point<dim>&          direction);
     
-                                    /**
-                                     * @deprecated Use
-                                     * compute_downstream() instead
-                                     */
-    template <class DH, int dim>
-    static void
-    compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
-                          std::vector<unsigned int>& reverse,
-                          const DH&                  dof_handler,
-                          const Point<dim>&          direction);
+                                  /**
+                                   * @deprecated Use
+                                   * compute_downstream() instead
+                                   */
+  template <class DH, int dim>
+  void
+  compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
+                        std::vector<unsigned int>& reverse,
+                        const DH&                  dof_handler,
+                        const Point<dim>&          direction);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * downstream_dg() function. Does
-                                     * not perform the renumbering on
-                                     * the MGDoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <int dim>
-    static void
-    compute_downstream (std::vector<unsigned int>& new_dof_indices,
-                       std::vector<unsigned int>& reverse,
-                       const MGDoFHandler<dim>&   dof_handler,
-                       const unsigned int         level,
-                       const Point<dim>&          direction);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * downstream_dg() function. Does
+                                   * not perform the renumbering on
+                                   * the MGDoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <int dim>
+  void
+  compute_downstream (std::vector<unsigned int>& new_dof_indices,
+                     std::vector<unsigned int>& reverse,
+                     const MGDoFHandler<dim>&   dof_handler,
+                     const unsigned int         level,
+                     const Point<dim>&          direction);
     
-                                    /**
-                                     * @deprecated Use
-                                     * compute_downstream() instead
-                                     */
-    template <int dim>
-    static void
-    compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
-                          std::vector<unsigned int>& reverse,
-                          const MGDoFHandler<dim>&   dof_handler,
-                          const unsigned int         level,
-                          const Point<dim>&          direction);
+                                  /**
+                                   * @deprecated Use
+                                   * compute_downstream() instead
+                                   */
+  template <int dim>
+  void
+  compute_downstream_dg (std::vector<unsigned int>& new_dof_indices,
+                        std::vector<unsigned int>& reverse,
+                        const MGDoFHandler<dim>&   dof_handler,
+                        const unsigned int         level,
+                        const Point<dim>&          direction);
 
-                                    /**
-                                     * Cell-wise clockwise numbering.
-                                     *
-                                     * This function produces a
-                                     * (counter)clockwise ordering of
-                                     * the mesh cells with respect to
-                                     * the hub @p center and calls
-                                     * cell_wise_dg().  Therefore, it
-                                     * only works with Discontinuous
-                                     * Galerkin Finite Elements,
-                                     * i.e. all degrees of freedom
-                                     * have to be associated with the
-                                     * interior of the cell.
-                                     */
-    template <class DH, int dim>
-    static void
-    clockwise_dg (DH&               dof_handler,
-                 const Point<dim>& center,
-                 const bool        counter = false);
+                                  /**
+                                   * Cell-wise clockwise numbering.
+                                   *
+                                   * This function produces a
+                                   * (counter)clockwise ordering of
+                                   * the mesh cells with respect to
+                                   * the hub @p center and calls
+                                   * cell_wise_dg().  Therefore, it
+                                   * only works with Discontinuous
+                                   * Galerkin Finite Elements,
+                                   * i.e. all degrees of freedom
+                                   * have to be associated with the
+                                   * interior of the cell.
+                                   */
+  template <class DH, int dim>
+  void
+  clockwise_dg (DH&               dof_handler,
+               const Point<dim>& center,
+               const bool        counter = false);
 
-                                    /**
-                                     * Cell-wise clockwise numbering
-                                     * on one level. See the other
-                                     * function with the same name.
-                                     */
-    template <int dim>
-    static void
-    clockwise_dg (MGDoFHandler<dim>  &dof_handler,
-                 const unsigned int level,
-                 const Point<dim> &center,
-                 const bool counter = false);
+                                  /**
+                                   * Cell-wise clockwise numbering
+                                   * on one level. See the other
+                                   * function with the same name.
+                                   */
+  template <int dim>
+  void
+  clockwise_dg (MGDoFHandler<dim>  &dof_handler,
+               const unsigned int level,
+               const Point<dim> &center,
+               const bool counter = false);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * clockwise_dg() functions. Does
-                                     * not perform the renumbering on
-                                     * the DoFHandler dofs but
-                                     * returns the renumbering
-                                     * vector.
-                                     */
-    template <class DH, int dim>
-    static void
-    compute_clockwise_dg (std::vector<unsigned int>& new_dof_indices,
-                         const DH&                  dof_handler,
-                         const Point<dim>&          center,
-                         const bool                 counter);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * clockwise_dg() functions. Does
+                                   * not perform the renumbering on
+                                   * the DoFHandler dofs but
+                                   * returns the renumbering
+                                   * vector.
+                                   */
+  template <class DH, int dim>
+  void
+  compute_clockwise_dg (std::vector<unsigned int>& new_dof_indices,
+                       const DH&                  dof_handler,
+                       const Point<dim>&          center,
+                       const bool                 counter);
 
-                                    /**
-                                     * Sort those degrees of freedom
-                                     * which are tagged with @p true
-                                     * in the @p selected_dofs array
-                                     * to the back of the DoF
-                                     * numbers. The sorting is
-                                     * stable, i.e. the relative
-                                     * order within the tagged
-                                     * degrees of freedom is
-                                     * preserved, as is the relative
-                                     * order within the untagged
-                                     * ones.
-                                     */
-    template <class DH>
-    static void
-    sort_selected_dofs_back (DH&                      dof_handler,
-                            const std::vector<bool>& selected_dofs);
+                                  /**
+                                   * Sort those degrees of freedom
+                                   * which are tagged with @p true
+                                   * in the @p selected_dofs array
+                                   * to the back of the DoF
+                                   * numbers. The sorting is
+                                   * stable, i.e. the relative
+                                   * order within the tagged
+                                   * degrees of freedom is
+                                   * preserved, as is the relative
+                                   * order within the untagged
+                                   * ones.
+                                   */
+  template <class DH>
+  void
+  sort_selected_dofs_back (DH&                      dof_handler,
+                          const std::vector<bool>& selected_dofs);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * sort_selected_dofs_back()
-                                     * function. Does not perform the
-                                     * renumbering on the DoFHandler
-                                     * dofs but returns the
-                                     * renumbering vector.
-                                     */
-    template <class DH>
-    static void
-    compute_sort_selected_dofs_back (std::vector<unsigned int>& new_dof_indices,
-                                    const DH&                  dof_handler,
-                                    const std::vector<bool>&   selected_dofs);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * sort_selected_dofs_back()
+                                   * function. Does not perform the
+                                   * renumbering on the DoFHandler
+                                   * dofs but returns the
+                                   * renumbering vector.
+                                   */
+  template <class DH>
+  void
+  compute_sort_selected_dofs_back (std::vector<unsigned int>& new_dof_indices,
+                                  const DH&                  dof_handler,
+                                  const std::vector<bool>&   selected_dofs);
 
-                                    /**
-                                     * Renumber the degrees of
-                                     * freedom in a random way.
-                                     */
-    template <class DH>
-    static void
-    random (DH& dof_handler);
+                                  /**
+                                   * Renumber the degrees of
+                                   * freedom in a random way.
+                                   */
+  template <class DH>
+  void
+  random (DH& dof_handler);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the random()
-                                     * function. Does not perform the
-                                     * renumbering on the DoFHandler
-                                     * dofs but returns the
-                                     * renumbering vector.
-                                     */   
-    template <class DH>
-    static void
-    compute_random (std::vector<unsigned int> &new_dof_indices,
-                   const DH& dof_handler);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the random()
+                                   * function. Does not perform the
+                                   * renumbering on the DoFHandler
+                                   * dofs but returns the
+                                   * renumbering vector.
+                                   */   
+  template <class DH>
+  void
+  compute_random (std::vector<unsigned int> &new_dof_indices,
+                 const DH& dof_handler);
 
-                                    /**
-                                     * Renumber the degrees of
-                                     * freedom such that they are
-                                     * associated with the subdomain
-                                     * id of the cells they are
-                                     * living on, i.e. first all
-                                     * degrees of freedom that belong
-                                     * to cells with subdomain zero,
-                                     * then all with subdomain one,
-                                     * etc. This is useful when doing
-                                     * parallel computations after
-                                     * assigning subdomain ids using
-                                     * a partitioner (see the
-                                     * GridTools::partition_triangulation
-                                     * function for this).
-                                     *
-                                     * Note that degrees of freedom
-                                     * associated with faces, edges,
-                                     * and vertices may be associated
-                                     * with multiple subdomains if
-                                     * they are sitting on partition
-                                     * boundaries. It would therefore
-                                     * be undefined with which
-                                     * subdomain they have to be
-                                     * associated. For this, we use
-                                     * what we get from the
-                                     * DoFTools::get_subdomain_association
-                                     * function.
-                                     *
-                                     * The algorithm is stable, i.e. if
-                                      * two dofs i,j have <tt>i<j</tt> and belong
-                                      * to the same subdomain, then they
-                                      * will be in this order also after
-                                      * reordering.
-                                     */
-    template <int dim>
-    static void
-    subdomain_wise (DoFHandler<dim> &dof_handler);
+                                  /**
+                                   * Renumber the degrees of
+                                   * freedom such that they are
+                                   * associated with the subdomain
+                                   * id of the cells they are
+                                   * living on, i.e. first all
+                                   * degrees of freedom that belong
+                                   * to cells with subdomain zero,
+                                   * then all with subdomain one,
+                                   * etc. This is useful when doing
+                                   * parallel computations after
+                                   * assigning subdomain ids using
+                                   * a partitioner (see the
+                                   * GridTools::partition_triangulation
+                                   * function for this).
+                                   *
+                                   * Note that degrees of freedom
+                                   * associated with faces, edges,
+                                   * and vertices may be associated
+                                   * with multiple subdomains if
+                                   * they are sitting on partition
+                                   * boundaries. It would therefore
+                                   * be undefined with which
+                                   * subdomain they have to be
+                                   * associated. For this, we use
+                                   * what we get from the
+                                   * DoFTools::get_subdomain_association
+                                   * function.
+                                   *
+                                   * The algorithm is stable, i.e. if
+                                   * two dofs i,j have <tt>i<j</tt> and belong
+                                   * to the same subdomain, then they
+                                   * will be in this order also after
+                                   * reordering.
+                                   */
+  template <int dim>
+  void
+  subdomain_wise (DoFHandler<dim> &dof_handler);
 
-                                    /**
-                                     * Computes the renumbering
-                                     * vector needed by the
-                                     * subdomain_wise()
-                                     * function. Does not perform the
-                                     * renumbering on the @p
-                                     * DoFHandler dofs but returns
-                                     * the renumbering vector.
-                                     */   
-    template <int dim>
-    static void
-    compute_subdomain_wise (std::vector<unsigned int> &new_dof_indices,
-                            const DoFHandler<dim>     &dof_handler);
+                                  /**
+                                   * Computes the renumbering
+                                   * vector needed by the
+                                   * subdomain_wise()
+                                   * function. Does not perform the
+                                   * renumbering on the @p
+                                   * DoFHandler dofs but returns
+                                   * the renumbering vector.
+                                   */   
+  template <int dim>
+  void
+  compute_subdomain_wise (std::vector<unsigned int> &new_dof_indices,
+                         const DoFHandler<dim>     &dof_handler);
     
-                                    /**
-                                     * Exception
-                                     *
-                                     * @ingroup Exceptions
-                                     */
-    DeclException0 (ExcRenumberingIncomplete);
-                                    /**
-                                     * Exception
-                                     *
-                                     * @ingroup Exceptions
-                                     */
-    DeclException0 (ExcInvalidComponentOrder);
-                                    /**
-                                     * The function is only
-                                     * implemented for Discontinuous
-                                     * Galerkin Finite elements.
-                                     *
-                                     * @ingroup Exceptions
-                                     */
-    DeclException0 (ExcNotDGFEM);
-};
+                                  /**
+                                   * Exception
+                                   *
+                                   * @ingroup Exceptions
+                                   */
+  DeclException0 (ExcRenumberingIncomplete);
+                                  /**
+                                   * Exception
+                                   *
+                                   * @ingroup Exceptions
+                                   */
+  DeclException0 (ExcInvalidComponentOrder);
+                                  /**
+                                   * The function is only
+                                   * implemented for Discontinuous
+                                   * Galerkin Finite elements.
+                                   *
+                                   * @ingroup Exceptions
+                                   */
+  DeclException0 (ExcNotDGFEM);
+}
 
 DEAL_II_NAMESPACE_CLOSE
 
index b0e1eebd95b9405cc024d52094e76286ad8f8c0f..b38a596eabaac8f0e04159c18fe1451ca8770267 100644 (file)
@@ -63,55 +63,83 @@ extern "C" long int lrand48 (void);
 DEAL_II_NAMESPACE_OPEN
 
 
+namespace internal
+{
+                                  // the following function is sorta
+                                  // stupid: gcc complains that the
+                                  // variables _3, _4, ... declared
+                                  // in
+                                  // contrib/boost/include/boost/bind/placeholders.hpp
+                                  // are unused. _1 and _2 are
+                                  // apparently used, so this isn't a
+                                  // problem. to suppress the warning
+                                  // in this case, write a dummy
+                                  // function that simply references
+                                  // these variables in an
+                                  // unobtrusive way...
+  void avoid_warning_about_unused_variables ()
+  {
+    std::printf ("%p %p %p %p %p %p %p",
+                (void*)&_3, (void*)&_4,
+                (void*)&_5,(void*) &_6,
+                (void*)&_7,(void*) &_8,
+                (void*)&_9);
+  }
+}
+
 
+namespace DoFRenumbering
+{
+  namespace internal
+  {
 // The following two classes are defined to be used in the compute_*
 // functions. Using them allows to use the same function to compute
 // level numbering for multigrid as well as numbering for the global
 // system matrix.
-template <class T>
-class WrapDoFIterator : private T
-{
-  public:
-    WrapDoFIterator (const T& t) : T(t) {}
-
-    void get_dof_indices (std::vector<unsigned int>& v) const
-      {
-       (*this)->get_dof_indices(v);
-      }
-
-    template <class T2>
-    bool operator != (const T2& i) const
-      {
-       return (! (T::operator==(i)));
-      }
-                                    // Allow access to these operators of T
-    T::operator->;
-    T::operator++;
-    T::operator==;
-};
-
-
+    template <class T>
+    class WrapDoFIterator : private T
+    {
+      public:
+       WrapDoFIterator (const T& t) : T(t) {}
+
+       void get_dof_indices (std::vector<unsigned int>& v) const
+         {
+           (*this)->get_dof_indices(v);
+         }
+
+       template <class T2>
+       bool operator != (const T2& i) const
+         {
+           return (! (T::operator==(i)));
+         }
+                                        // Allow access to these operators of T
+       T::operator->;
+       T::operator++;
+       T::operator==;
+    };
 
-template <class T>
-class WrapMGDoFIterator : private T
-{
-  public:
-    WrapMGDoFIterator (const T& t) : T(t) {}
 
-    void get_dof_indices (std::vector<unsigned int>& v) const
-      {
-       (*this)->get_mg_dof_indices(v);
-      }
 
-    bool operator != (const WrapMGDoFIterator<T>& i) const
-      {
-       return (! (T::operator==(i)));
-      }
-                                    // Allow access to these operators of T
-    T::operator->;
-    T::operator++;
-    T::operator==;
-};
+    template <class T>
+    class WrapMGDoFIterator : private T
+    {
+      public:
+       WrapMGDoFIterator (const T& t) : T(t) {}
+
+       void get_dof_indices (std::vector<unsigned int>& v) const
+         {
+           (*this)->get_mg_dof_indices(v);
+         }
+
+       bool operator != (const WrapMGDoFIterator<T>& i) const
+         {
+           return (! (T::operator==(i)));
+         }
+                                        // Allow access to these operators of T
+       T::operator->;
+       T::operator++;
+       T::operator==;
+    };
 
 
 
@@ -120,2016 +148,2032 @@ class WrapMGDoFIterator : private T
  * of the second cell is downstream of the center of the first one with
  * respect to the direction given to the constructor.
  */
-template <int dim>
-struct CompareDownstream
-{
-                                    /**
-                                     * Constructor.
-                                     */
-    CompareDownstream (const Point<dim> &dir)
-                    :
-                   dir(dir) 
-      {}
-                                    /**
-                                     * Return true if c1 < c2.
-                                     */
-    template <class DHCellIterator>
-    bool operator () (const DHCellIterator &c1,
-                     const DHCellIterator &c2) const
-      {
-       const Point<dim> diff = c2->center() - c1->center();
-       return (diff*dir > 0);
-      }
-
-  private:
-                                    /**
-                                     * Flow direction.
-                                     */
-    const Point<dim> dir;
-};
-
-
-namespace types
-{
-  using namespace boost;
-  using namespace std;
+    template <int dim>
+    struct CompareDownstream
+    {
+                                        /**
+                                         * Constructor.
+                                         */
+       CompareDownstream (const Point<dim> &dir)
+                       :
+                       dir(dir) 
+         {}
+                                        /**
+                                         * Return true if c1 < c2.
+                                         */
+       template <class DHCellIterator>
+       bool operator () (const DHCellIterator &c1,
+                         const DHCellIterator &c2) const
+         {
+           const Point<dim> diff = c2->center() - c1->center();
+           return (diff*dir > 0);
+         }
+
+      private:
+                                        /**
+                                         * Flow direction.
+                                         */
+       const Point<dim> dir;
+    };
+  }
+  
+  namespace boost
+  {
+    namespace types
+    {
+      using namespace ::boost;
+      using namespace std;
 
-  typedef vecS x;
+      typedef vecS x;
     
-  typedef adjacency_list<vecS, vecS, undirectedS, 
-                        property<vertex_color_t, default_color_type,
-                                 property<vertex_degree_t,int> > > Graph;
-  typedef graph_traits<Graph>::vertex_descriptor Vertex;
-  typedef graph_traits<Graph>::vertices_size_type size_type;
+      typedef adjacency_list<vecS, vecS, undirectedS, 
+                            property<vertex_color_t, default_color_type,
+                                     property<vertex_degree_t,int> > > Graph;
+      typedef graph_traits<Graph>::vertex_descriptor Vertex;
+      typedef graph_traits<Graph>::vertices_size_type size_type;
 
-  typedef std::pair<size_type, size_type> Pair;
-}
+      typedef std::pair<size_type, size_type> Pair;
+    }
   
 
-namespace
-{
-  template <class DH>
-  void create_graph (const DH                                                       &dof_handler,
-                    const bool                                                      use_constraints,
-                    types::Graph                                                   &graph,
-                    types::property_map<types::Graph,types::vertex_degree_t>::type &graph_degree)
-  {
-    Assert (use_constraints == false, ExcNotImplemented());
-    
-    std::vector<unsigned int> dofs_on_this_cell;
-    dofs_on_this_cell.reserve (DoFTools::max_dofs_per_cell(dof_handler));
-    typename DH::active_cell_iterator cell = dof_handler.begin_active(),
-                                     endc = dof_handler.end();
-    for (; cell!=endc; ++cell) 
+    namespace internal
+    {
+      template <class DH>
+      void create_graph (const DH                                                       &dof_handler,
+                        const bool                                                      use_constraints,
+                        types::Graph                                                   &graph,
+                        types::property_map<types::Graph,types::vertex_degree_t>::type &graph_degree)
       {
-       const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
-       dofs_on_this_cell.resize (dofs_per_cell);
-       cell->get_dof_indices (dofs_on_this_cell);
+       Assert (use_constraints == false, ExcNotImplemented());
+    
+       std::vector<unsigned int> dofs_on_this_cell;
+       dofs_on_this_cell.reserve (DoFTools::max_dofs_per_cell(dof_handler));
+       typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+                                         endc = dof_handler.end();
+       for (; cell!=endc; ++cell) 
+         {
+           const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+           dofs_on_this_cell.resize (dofs_per_cell);
+           cell->get_dof_indices (dofs_on_this_cell);
 
-       for (unsigned int i=0; i<dofs_per_cell; ++i)
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           add_edge(dofs_on_this_cell[i], dofs_on_this_cell[j], graph);
-      }
+           for (unsigned int i=0; i<dofs_per_cell; ++i)
+             for (unsigned int j=0; j<dofs_per_cell; ++j)
+               add_edge(dofs_on_this_cell[i], dofs_on_this_cell[j], graph);
+         }
 
 
-    types::graph_traits<types::Graph>::vertex_iterator ui, ui_end;
+       types::graph_traits<types::Graph>::vertex_iterator ui, ui_end;
 
-    graph_degree = get(boost::vertex_degree, graph);
-    for (boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
-      graph_degree[*ui] = degree(*ui, graph);
+       graph_degree = get(::boost::vertex_degree, graph);
+       for (::boost::tie(ui, ui_end) = vertices(graph); ui != ui_end; ++ui)
+         graph_degree[*ui] = degree(*ui, graph);
 
-  }
-}
+      }
+    }
 
   
-template <class DH>
-void
-DoFRenumbering::boost::
-Cuthill_McKee (DH&              dof_handler,
-              const bool       reversed_numbering,
-              const bool       use_constraints)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
-                       use_constraints);
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  dof_handler.renumber_dofs (renumbering);
-}
+    template <class DH>
+    void
+    Cuthill_McKee (DH&              dof_handler,
+                  const bool       reversed_numbering,
+                  const bool       use_constraints)
+    {
+      std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                           DH::invalid_dof_index);
+      compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
+                           use_constraints);
+
+                                      // actually perform renumbering;
+                                      // this is dimension specific and
+                                      // thus needs an own function
+      dof_handler.renumber_dofs (renumbering);
+    }
 
 
-template <class DH>
-void
-DoFRenumbering::boost::
-compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
-                      const DH        &dof_handler,
-                      const bool       reversed_numbering,
-                      const bool       use_constraints)
-{
-  types::Graph
-    graph(dof_handler.n_dofs());
-  types::property_map<types::Graph,types::vertex_degree_t>::type
-    graph_degree;
+    template <class DH>
+    void
+    compute_Cuthill_McKee (std::vector<unsigned int>& new_dof_indices,
+                          const DH        &dof_handler,
+                          const bool       reversed_numbering,
+                          const bool       use_constraints)
+    {
+      types::Graph
+       graph(dof_handler.n_dofs());
+      types::property_map<types::Graph,types::vertex_degree_t>::type
+       graph_degree;
 
-  create_graph (dof_handler, use_constraints, graph, graph_degree);
+      internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
     
-  types::property_map<types::Graph, types::vertex_index_t>::type
-    index_map = get(::boost::vertex_index, graph);
+      types::property_map<types::Graph, types::vertex_index_t>::type
+       index_map = get(::boost::vertex_index, graph);
 
 
-  std::vector<types::Vertex> inv_perm(num_vertices(graph));
+      std::vector<types::Vertex> inv_perm(num_vertices(graph));
 
-  if (reversed_numbering == false)
-    ::boost::cuthill_mckee_ordering(graph, inv_perm.rbegin(),
-                                   get(::boost::vertex_color, graph),
-                                   make_degree_map(graph));
-  else
-    ::boost::cuthill_mckee_ordering(graph, inv_perm.begin(),
-                                   get(::boost::vertex_color, graph),
-                                   make_degree_map(graph));
+      if (reversed_numbering == false)
+       ::boost::cuthill_mckee_ordering(graph, inv_perm.rbegin(),
+                                       get(::boost::vertex_color, graph),
+                                       make_degree_map(graph));
+      else
+       ::boost::cuthill_mckee_ordering(graph, inv_perm.begin(),
+                                       get(::boost::vertex_color, graph),
+                                       make_degree_map(graph));
     
-  for (types::size_type c = 0; c != inv_perm.size(); ++c)
-    new_dof_indices[index_map[inv_perm[c]]] = c;
+      for (types::size_type c = 0; c != inv_perm.size(); ++c)
+       new_dof_indices[index_map[inv_perm[c]]] = c;
 
-  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
-                    DH::invalid_dof_index) == new_dof_indices.end(),
-         ExcInternalError());
-}
+      Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
+                        DH::invalid_dof_index) == new_dof_indices.end(),
+             ExcInternalError());
+    }
 
 
 
-template <class DH>
-void
-DoFRenumbering::boost::
-king_ordering (DH&              dof_handler,
-              const bool       reversed_numbering,
-              const bool       use_constraints)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_king_ordering(renumbering, dof_handler, reversed_numbering,
-                       use_constraints);
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  dof_handler.renumber_dofs (renumbering);
-}
+    template <class DH>
+    void
+    king_ordering (DH&              dof_handler,
+                  const bool       reversed_numbering,
+                  const bool       use_constraints)
+    {
+      std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                           DH::invalid_dof_index);
+      compute_king_ordering(renumbering, dof_handler, reversed_numbering,
+                           use_constraints);
+
+                                      // actually perform renumbering;
+                                      // this is dimension specific and
+                                      // thus needs an own function
+      dof_handler.renumber_dofs (renumbering);
+    }
 
 
-template <class DH>
-void
-DoFRenumbering::boost::
-compute_king_ordering (std::vector<unsigned int>& new_dof_indices,
-                      const DH        &dof_handler,
-                      const bool       reversed_numbering,
-                      const bool       use_constraints)
-{
-  types::Graph
-    graph(dof_handler.n_dofs());
-  types::property_map<types::Graph,types::vertex_degree_t>::type
-    graph_degree;
+    template <class DH>
+    void
+    compute_king_ordering (std::vector<unsigned int>& new_dof_indices,
+                          const DH        &dof_handler,
+                          const bool       reversed_numbering,
+                          const bool       use_constraints)
+    {
+      types::Graph
+       graph(dof_handler.n_dofs());
+      types::property_map<types::Graph,types::vertex_degree_t>::type
+       graph_degree;
 
-  create_graph (dof_handler, use_constraints, graph, graph_degree);
+      internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
     
-  types::property_map<types::Graph, types::vertex_index_t>::type
-    index_map = get(::boost::vertex_index, graph);
+      types::property_map<types::Graph, types::vertex_index_t>::type
+       index_map = get(::boost::vertex_index, graph);
 
 
-  std::vector<types::Vertex> inv_perm(num_vertices(graph));
+      std::vector<types::Vertex> inv_perm(num_vertices(graph));
 
-  if (reversed_numbering == false)
-    ::boost::king_ordering(graph, inv_perm.rbegin());
-  else
-    ::boost::king_ordering(graph, inv_perm.begin());
+      if (reversed_numbering == false)
+       ::boost::king_ordering(graph, inv_perm.rbegin());
+      else
+       ::boost::king_ordering(graph, inv_perm.begin());
     
-  for (types::size_type c = 0; c != inv_perm.size(); ++c)
-    new_dof_indices[index_map[inv_perm[c]]] = c;
+      for (types::size_type c = 0; c != inv_perm.size(); ++c)
+       new_dof_indices[index_map[inv_perm[c]]] = c;
 
-  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
-                    DH::invalid_dof_index) == new_dof_indices.end(),
-         ExcInternalError());
-}
+      Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
+                        DH::invalid_dof_index) == new_dof_indices.end(),
+             ExcInternalError());
+    }
 
 
 
-template <class DH>
-void
-DoFRenumbering::boost::
-minimum_degree (DH&              dof_handler,
-               const bool       reversed_numbering,
-               const bool       use_constraints)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_minimum_degree(renumbering, dof_handler, reversed_numbering,
-                        use_constraints);
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  dof_handler.renumber_dofs (renumbering);
-}
+    template <class DH>
+    void
+    minimum_degree (DH&              dof_handler,
+                   const bool       reversed_numbering,
+                   const bool       use_constraints)
+    {
+      std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                           DH::invalid_dof_index);
+      compute_minimum_degree(renumbering, dof_handler, reversed_numbering,
+                            use_constraints);
+
+                                      // actually perform renumbering;
+                                      // this is dimension specific and
+                                      // thus needs an own function
+      dof_handler.renumber_dofs (renumbering);
+    }
 
 
-template <class DH>
-void
-DoFRenumbering::boost::
-compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
-                       const DH        &dof_handler,
-                       const bool       reversed_numbering,
-                       const bool       use_constraints)
-{
-  Assert (use_constraints == false, ExcNotImplemented());
+    template <class DH>
+    void
+    compute_minimum_degree (std::vector<unsigned int>& new_dof_indices,
+                           const DH        &dof_handler,
+                           const bool       reversed_numbering,
+                           const bool       use_constraints)
+    {
+      Assert (use_constraints == false, ExcNotImplemented());
 
-                                  // the following code is pretty
-                                  // much a verbatim copy of the
-                                  // sample code for the
-                                  // minimum_degree_ordering manual
-                                  // page from the BOOST Graph
-                                  // Library
-  using namespace boost;
+                                      // the following code is pretty
+                                      // much a verbatim copy of the
+                                      // sample code for the
+                                      // minimum_degree_ordering manual
+                                      // page from the BOOST Graph
+                                      // Library
+      using namespace ::boost;
 
-  int delta = 0;
+      int delta = 0;
 
-  typedef double Type;
+      typedef double Type;
 
-                                  // must be BGL directed graph now
-  typedef adjacency_list<vecS, vecS, directedS>  Graph;
-  typedef graph_traits<Graph>::vertex_descriptor Vertex;
+                                      // must be BGL directed graph now
+      typedef adjacency_list<vecS, vecS, directedS>  Graph;
+      typedef graph_traits<Graph>::vertex_descriptor Vertex;
 
-  int n = dof_handler.n_dofs();
+      int n = dof_handler.n_dofs();
 
-  Graph G(n);
+      Graph G(n);
 
-  std::vector<unsigned int> dofs_on_this_cell;
+      std::vector<unsigned int> dofs_on_this_cell;
 
-  typename DH::active_cell_iterator cell = dof_handler.begin_active(),
-                                   endc = dof_handler.end();
+      typename DH::active_cell_iterator cell = dof_handler.begin_active(),
+                                       endc = dof_handler.end();
 
-  for (; cell!=endc; ++cell)
-    {
+      for (; cell!=endc; ++cell)
+       {
 
-      const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
+         const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
 
-      dofs_on_this_cell.resize (dofs_per_cell);
+         dofs_on_this_cell.resize (dofs_per_cell);
 
-      cell->get_dof_indices (dofs_on_this_cell);
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
-           {
-             add_edge (dofs_on_this_cell[i], dofs_on_this_cell[j], G);
-             add_edge (dofs_on_this_cell[j], dofs_on_this_cell[i], G);
-           }
-    }
+         cell->get_dof_indices (dofs_on_this_cell);
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           for (unsigned int j=0; j<dofs_per_cell; ++j)
+             if (dofs_on_this_cell[i] > dofs_on_this_cell[j])
+               {
+                 add_edge (dofs_on_this_cell[i], dofs_on_this_cell[j], G);
+                 add_edge (dofs_on_this_cell[j], dofs_on_this_cell[i], G);
+               }
+       }
   
   
-  typedef std::vector<int> Vector;
+      typedef std::vector<int> Vector;
 
 
-  Vector inverse_perm(n, 0);
+      Vector inverse_perm(n, 0);
 
-  Vector perm(n, 0);
+      Vector perm(n, 0);
 
 
-  Vector supernode_sizes(n, 1);
-                                  // init has to be 1
+      Vector supernode_sizes(n, 1);
+                                      // init has to be 1
 
-  ::boost::property_map<Graph, vertex_index_t>::type
-    id = get(vertex_index, G);
+      ::boost::property_map<Graph, vertex_index_t>::type
+         id = get(vertex_index, G);
 
 
-  Vector degree(n, 0);
+      Vector degree(n, 0);
 
 
-    minimum_degree_ordering
-      (G,
-       make_iterator_property_map(&degree[0], id, degree[0]),
-       &inverse_perm[0],
-       &perm[0],
-       make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]),
-       delta, id);
+      minimum_degree_ordering
+       (G,
+        make_iterator_property_map(&degree[0], id, degree[0]),
+        &inverse_perm[0],
+        &perm[0],
+        make_iterator_property_map(&supernode_sizes[0], id, supernode_sizes[0]),
+        delta, id);
 
 
-    for (int i=0; i<n; ++i)
-      {
-       Assert (std::find (perm.begin(), perm.end(), i)
-               != inverse_perm.end(),
-               ExcInternalError());
-       Assert (std::find (inverse_perm.begin(), inverse_perm.end(), i)
-               != inverse_perm.end(),
-               ExcInternalError());
-       Assert (inverse_perm[perm[i]] == i, ExcInternalError());
-      }
-
-    if (reversed_numbering == true)
-      std::copy (perm.begin(), perm.end(),
-                new_dof_indices.begin());
-    else
-      std::copy (inverse_perm.begin(), inverse_perm.end(),
-                new_dof_indices.begin());
-}
+      for (int i=0; i<n; ++i)
+       {
+         Assert (std::find (perm.begin(), perm.end(), i)
+                 != inverse_perm.end(),
+                 ExcInternalError());
+         Assert (std::find (inverse_perm.begin(), inverse_perm.end(), i)
+                 != inverse_perm.end(),
+                 ExcInternalError());
+         Assert (inverse_perm[perm[i]] == i, ExcInternalError());
+       }
 
+      if (reversed_numbering == true)
+       std::copy (perm.begin(), perm.end(),
+                  new_dof_indices.begin());
+      else
+       std::copy (inverse_perm.begin(), inverse_perm.end(),
+                  new_dof_indices.begin());
+    }
 
+  }  // namespace boost
   
 
 
-template <class DH>
-void
-DoFRenumbering::
-Cuthill_McKee (DH&              dof_handler,
-               const bool       reversed_numbering,
-               const bool       use_constraints,
-               const std::vector<unsigned int> &starting_indices)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
-                       use_constraints, starting_indices);
-
-                                  // actually perform renumbering;
-                                  // this is dimension specific and
-                                  // thus needs an own function
-  dof_handler.renumber_dofs (renumbering);
-}
+  template <class DH>
+  void
+  Cuthill_McKee (DH&              dof_handler,
+                const bool       reversed_numbering,
+                const bool       use_constraints,
+                const std::vector<unsigned int> &starting_indices)
+  {
+    std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                         DH::invalid_dof_index);
+    compute_Cuthill_McKee(renumbering, dof_handler, reversed_numbering,
+                         use_constraints, starting_indices);
+
+                                    // actually perform renumbering;
+                                    // this is dimension specific and
+                                    // thus needs an own function
+    dof_handler.renumber_dofs (renumbering);
+  }
 
 
 
-template <class DH>
-void
-DoFRenumbering::
-compute_Cuthill_McKee (std::vector<unsigned int>& new_indices,
-                       const DH&                  dof_handler,
-                       const bool                 reversed_numbering,
-                       const bool                 use_constraints,
-                       const std::vector<unsigned int>& starting_indices)
-{
-                                  // make the connection graph. in
-                                  // more than 2d use an intermediate
-                                  // compressed sparsity pattern
-                                  // since the we don't have very
-                                  // good estimates for
-                                  // max_couplings_between_dofs() in
-                                  // 3d and this then leads to
-                                  // excessive memory consumption
-                                  //
-                                  // note that if constraints are not
-                                  // requested, then the
-                                  // 'constraints' object will be
-                                  // empty, and calling condense with
-                                  // it is a no-op
-  ConstraintMatrix constraints;
-  if (use_constraints) 
-    DoFTools::make_hanging_node_constraints (dof_handler, constraints);
-  constraints.close ();
+  template <class DH>
+  void
+  compute_Cuthill_McKee (std::vector<unsigned int>& new_indices,
+                        const DH&                  dof_handler,
+                        const bool                 reversed_numbering,
+                        const bool                 use_constraints,
+                        const std::vector<unsigned int>& starting_indices)
+  {
+                                    // make the connection graph. in
+                                    // more than 2d use an intermediate
+                                    // compressed sparsity pattern
+                                    // since the we don't have very
+                                    // good estimates for
+                                    // max_couplings_between_dofs() in
+                                    // 3d and this then leads to
+                                    // excessive memory consumption
+                                    //
+                                    // note that if constraints are not
+                                    // requested, then the
+                                    // 'constraints' object will be
+                                    // empty, and calling condense with
+                                    // it is a no-op
+    ConstraintMatrix constraints;
+    if (use_constraints) 
+      DoFTools::make_hanging_node_constraints (dof_handler, constraints);
+    constraints.close ();
     
-  SparsityPattern sparsity;
-  if (DH::dimension <= 2)
-    {
-      sparsity.reinit (dof_handler.n_dofs(),
-                      dof_handler.n_dofs(),
-                      dof_handler.max_couplings_between_dofs());
-      DoFTools::make_sparsity_pattern (dof_handler, sparsity);
-      constraints.condense (sparsity);
-    }
-  else
-    {
-      CompressedSparsityPattern csp (dof_handler.n_dofs(),
-                                    dof_handler.n_dofs());
-      DoFTools::make_sparsity_pattern (dof_handler, csp);
-      constraints.condense (csp);
-      sparsity.copy_from (csp);
-    }
+    SparsityPattern sparsity;
+    if (DH::dimension <= 2)
+      {
+       sparsity.reinit (dof_handler.n_dofs(),
+                        dof_handler.n_dofs(),
+                        dof_handler.max_couplings_between_dofs());
+       DoFTools::make_sparsity_pattern (dof_handler, sparsity);
+       constraints.condense (sparsity);
+      }
+    else
+      {
+       CompressedSparsityPattern csp (dof_handler.n_dofs(),
+                                      dof_handler.n_dofs());
+       DoFTools::make_sparsity_pattern (dof_handler, csp);
+       constraints.condense (csp);
+       sparsity.copy_from (csp);
+      }
 
-                                  // constraints are not needed anymore
-  constraints.clear ();
+                                    // constraints are not needed anymore
+    constraints.clear ();
   
-  const unsigned int n_dofs = sparsity.n_rows();
-                                  // store the new dof numbers;
-                                  // invalid_dof_index means that no new
-                                  // number was chosen yet
-  Assert(new_indices.size() == n_dofs,
-        ExcDimensionMismatch(new_indices.size(), n_dofs));
+    const unsigned int n_dofs = sparsity.n_rows();
+                                    // store the new dof numbers;
+                                    // invalid_dof_index means that no new
+                                    // number was chosen yet
+    Assert(new_indices.size() == n_dofs,
+          ExcDimensionMismatch(new_indices.size(), n_dofs));
   
-                                  // store the indices of the dofs renumbered
-                                  // in the last round. Default to starting
-                                  // points
-  std::vector<unsigned int> last_round_dofs (starting_indices);
+                                    // store the indices of the dofs renumbered
+                                    // in the last round. Default to starting
+                                    // points
+    std::vector<unsigned 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]==DH::invalid_dof_index) ||
-       (last_round_dofs[i]>=n_dofs))
-      last_round_dofs[i] = DH::invalid_dof_index;
+                                    // delete disallowed elements
+    for (unsigned int i=0; i<last_round_dofs.size(); ++i)
+      if ((last_round_dofs[i]==DH::invalid_dof_index) ||
+         (last_round_dofs[i]>=n_dofs))
+       last_round_dofs[i] = DH::invalid_dof_index;
   
-  std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-                 std::bind2nd(std::equal_to<unsigned int>(),
-                              DH::invalid_dof_index));
+    std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+                   std::bind2nd(std::equal_to<unsigned int>(),
+                                DH::invalid_dof_index));
   
-                                  // now if no valid points remain:
-                                  // find dof with lowest coordination
-                                  // number
+                                    // now if no valid points remain:
+                                    // find dof with lowest coordination
+                                    // number
   
-  if (last_round_dofs.size() == 0)
-    {
-      unsigned int starting_point   = DH::invalid_dof_index;
-      unsigned int min_coordination = n_dofs;
-      for (unsigned 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] == SparsityPattern::invalid_entry)
-             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;
-           };
-       };
+    if (last_round_dofs.size() == 0)
+      {
+       unsigned int starting_point   = DH::invalid_dof_index;
+       unsigned int min_coordination = n_dofs;
+       for (unsigned 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] == SparsityPattern::invalid_entry)
+               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 == DH::invalid_dof_index)
-       starting_point = 0;
+                                        // 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 == DH::invalid_dof_index)
+         starting_point = 0;
       
-                                      // initialize the first dof
-      last_round_dofs.push_back (starting_point);
-    };
+                                        // initialize the first dof
+       last_round_dofs.push_back (starting_point);
+      };
 
 
-                                  // store next free dof index
-  unsigned int next_free_number = 0;
+                                    // store next free dof index
+    unsigned int next_free_number = 0;
 
-                                  // enumerate the first round dofs
-  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-    new_indices[last_round_dofs[i]] = next_free_number++;
+                                    // enumerate the first round dofs
+    for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+      new_indices[last_round_dofs[i]] = next_free_number++;
 
-  bool all_dofs_renumbered = false;
+    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
-      std::vector<unsigned 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] == SparsityPattern::invalid_entry)
-           break;
-         else
-           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                      // sort dof numbers
-      std::sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                      // delete multiple entries
-      std::vector<unsigned int>::iterator end_sorted;
-      end_sorted = std::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_indices[next_round_dofs[s]] != DH::invalid_dof_index)
-         next_round_dofs.erase (next_round_dofs.begin() + 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
-      std::multimap<unsigned int, int> dofs_by_coordination;
-      
-                                      // find coordination number for
-                                      // each of these dofs
-      for (std::vector<unsigned 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)
+                                    // 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
+       std::vector<unsigned 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] == SparsityPattern::invalid_entry)
              break;
            else
-             ++coordination;
-
-                                          // insert this dof at its
-                                          // coordination number
-         const std::pair<const unsigned int, int> new_entry (coordination, *s);
-         dofs_by_coordination.insert (new_entry);
-       };
+             next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
       
-                                      // assign new DoF numbers to
-                                      // the elements of the present
-                                      // front:
-      std::multimap<unsigned int, int>::iterator i;
-      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-       new_indices[i->second] = next_free_number++;
-
-                                      // after that: copy this round's
-                                      // dofs for the next round
-      last_round_dofs = next_round_dofs;
-    };
+                                        // sort dof numbers
+       std::sort (next_round_dofs.begin(), next_round_dofs.end());
+
+                                        // delete multiple entries
+       std::vector<unsigned int>::iterator end_sorted;
+       end_sorted = std::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_indices[next_round_dofs[s]] != DH::invalid_dof_index)
+           next_round_dofs.erase (next_round_dofs.begin() + 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
+       std::multimap<unsigned int, int> dofs_by_coordination;
+      
+                                        // find coordination number for
+                                        // each of these dofs
+       for (std::vector<unsigned 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] == SparsityPattern::invalid_entry)
+               break;
+             else
+               ++coordination;
+
+                                            // insert this dof at its
+                                            // coordination number
+           const std::pair<const unsigned int, int> new_entry (coordination, *s);
+           dofs_by_coordination.insert (new_entry);
+         };
+      
+                                        // assign new DoF numbers to
+                                        // the elements of the present
+                                        // front:
+       std::multimap<unsigned int, int>::iterator i;
+       for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
+         new_indices[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. this mostly tests
-                                  // whether the
-                                  // front-marching-algorithm (which
-                                  // Cuthill-McKee actually is) has
-                                  // reached all points. it should
-                                  // usually do so, but might not for
-                                  // two reasons:
-                                  //
-                                  // - The algorithm above has a bug, or
-                                  // - The domain is not connected and
-                                  // consists of separate parts.
-                                  //
-                                  // In any case, if not all DoFs
-                                  // have been reached, renumbering
-                                  // will not be possible
-  if (std::find (new_indices.begin(), new_indices.end(), DH::invalid_dof_index)
-      !=
-      new_indices.end())
-    Assert (false, ExcRenumberingIncomplete());
-  Assert (next_free_number == n_dofs,
-         ExcRenumberingIncomplete());
+                                    // test for all indices
+                                    // numbered. this mostly tests
+                                    // whether the
+                                    // front-marching-algorithm (which
+                                    // Cuthill-McKee actually is) has
+                                    // reached all points. it should
+                                    // usually do so, but might not for
+                                    // two reasons:
+                                    //
+                                    // - The algorithm above has a bug, or
+                                    // - The domain is not connected and
+                                    // consists of separate parts.
+                                    //
+                                    // In any case, if not all DoFs
+                                    // have been reached, renumbering
+                                    // will not be possible
+    if (std::find (new_indices.begin(), new_indices.end(), DH::invalid_dof_index)
+       !=
+       new_indices.end())
+      Assert (false, ExcRenumberingIncomplete());
+    Assert (next_free_number == n_dofs,
+           ExcRenumberingIncomplete());
 #endif
 
-  if (reversed_numbering)
-    for (std::vector<unsigned int>::iterator i=new_indices.begin();
-        i!=new_indices.end(); ++i)
-      *i = n_dofs-*i-1;
-}
+    if (reversed_numbering)
+      for (std::vector<unsigned int>::iterator i=new_indices.begin();
+          i!=new_indices.end(); ++i)
+       *i = n_dofs-*i-1;
+  }
 
 
-template <int dim>
-void DoFRenumbering::Cuthill_McKee (
-  MGDoFHandler<dim>               &dof_handler,
-  const unsigned int               level,
-  const bool                       reversed_numbering,
-  const std::vector<unsigned int> &starting_indices)
-{
+  template <int dim>
+  void Cuthill_McKee (MGDoFHandler<dim>               &dof_handler,
+                     const unsigned int               level,
+                     const bool                       reversed_numbering,
+                     const std::vector<unsigned int> &starting_indices)
+  {
 //TODO: we should be doing the same here as in the other compute_CMK function to preserve some memory
   
-                                  // make the connection graph
-  SparsityPattern sparsity (dof_handler.n_dofs(level),
-                           dof_handler.max_couplings_between_dofs());
-  MGTools::make_sparsity_pattern (dof_handler, sparsity, level);
+                                    // make the connection graph
+    SparsityPattern sparsity (dof_handler.n_dofs(level),
+                             dof_handler.max_couplings_between_dofs());
+    MGTools::make_sparsity_pattern (dof_handler, sparsity, level);
     
-  const unsigned int n_dofs = sparsity.n_rows();
-                                  // store the new dof numbers; invalid_dof_index means
-                                  // that no new number was chosen yet
-  std::vector<unsigned int> new_indices(n_dofs, DoFHandler<dim>::invalid_dof_index);
-  
-                                  // store the indices of the dofs renumbered
-                                  // in the last round. Default to starting
-                                  // points
-  std::vector<unsigned 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]==DoFHandler<dim>::invalid_dof_index) ||
-       (last_round_dofs[i]>=n_dofs))
-      last_round_dofs[i] = DoFHandler<dim>::invalid_dof_index;
-  
-  std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
-                 std::bind2nd(std::equal_to<unsigned int>(),
-                              DoFHandler<dim>::invalid_dof_index));
-  
-                                  // now if no valid points remain:
-                                  // find dof with lowest coordination
-                                  // number
-  
-  if (last_round_dofs.size() == 0)
-    {
-      unsigned int starting_point   = DoFHandler<dim>::invalid_dof_index;
-      unsigned int min_coordination = n_dofs;
-      for (unsigned 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] == SparsityPattern::invalid_entry)
-             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);
-    };
+    const unsigned int n_dofs = sparsity.n_rows();
+                                    // store the new dof numbers; invalid_dof_index means
+                                    // that no new number was chosen yet
+    std::vector<unsigned int> new_indices(n_dofs, DoFHandler<dim>::invalid_dof_index);
+  
+                                    // store the indices of the dofs renumbered
+                                    // in the last round. Default to starting
+                                    // points
+    std::vector<unsigned 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]==DoFHandler<dim>::invalid_dof_index) ||
+         (last_round_dofs[i]>=n_dofs))
+       last_round_dofs[i] = DoFHandler<dim>::invalid_dof_index;
+  
+    std::remove_if (last_round_dofs.begin(), last_round_dofs.end(),
+                   std::bind2nd(std::equal_to<unsigned int>(),
+                                DoFHandler<dim>::invalid_dof_index));
+  
+                                    // now if no valid points remain:
+                                    // find dof with lowest coordination
+                                    // number
+  
+    if (last_round_dofs.size() == 0)
+      {
+       unsigned int starting_point   = DoFHandler<dim>::invalid_dof_index;
+       unsigned int min_coordination = n_dofs;
+       for (unsigned 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] == SparsityPattern::invalid_entry)
+               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
-  unsigned int next_free_number = 0;
+                                    // store next free dof index
+    unsigned int next_free_number = 0;
 
-                                  // enumerate the first round dofs
-  for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
-    new_indices[last_round_dofs[i]] = next_free_number++;
+                                    // enumerate the first round dofs
+    for (unsigned int i=0; i!=last_round_dofs.size(); ++i)
+      new_indices[last_round_dofs[i]] = next_free_number++;
 
-  bool all_dofs_renumbered = false;
+    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
-      std::vector<unsigned 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] == SparsityPattern::invalid_entry)
-           break;
-         else
-           next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
-      
-                                      // sort dof numbers
-      std::sort (next_round_dofs.begin(), next_round_dofs.end());
-
-                                      // delete multiple entries
-      std::vector<unsigned int>::iterator end_sorted;
-      end_sorted = std::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_indices[next_round_dofs[s]] != DoFHandler<dim>::invalid_dof_index)
-         next_round_dofs.erase (next_round_dofs.begin() + 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
-      std::multimap<unsigned int, int> dofs_by_coordination;
-      
-                                      // find coordination number for
-                                      // each of these dofs
-      for (std::vector<unsigned 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)
+                                    // 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
+       std::vector<unsigned 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] == SparsityPattern::invalid_entry)
              break;
            else
-             ++coordination;
-
-                                          // insert this dof at its
-                                          // coordination number
-         const std::pair<const unsigned int, int> new_entry (coordination, *s);
-         dofs_by_coordination.insert (new_entry);
-       };
+             next_round_dofs.push_back (sparsity.get_column_numbers()[j]);
       
-                                      ////
-      std::multimap<unsigned int, int>::iterator i;
-      for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
-       new_indices[i->second] = next_free_number++;
-
-                                      // after that: copy this round's
-                                      // dofs for the next round
-      last_round_dofs = next_round_dofs;
-    };
+                                        // sort dof numbers
+       std::sort (next_round_dofs.begin(), next_round_dofs.end());
+
+                                        // delete multiple entries
+       std::vector<unsigned int>::iterator end_sorted;
+       end_sorted = std::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_indices[next_round_dofs[s]] != DoFHandler<dim>::invalid_dof_index)
+           next_round_dofs.erase (next_round_dofs.begin() + 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
+       std::multimap<unsigned int, int> dofs_by_coordination;
+      
+                                        // find coordination number for
+                                        // each of these dofs
+       for (std::vector<unsigned 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] == SparsityPattern::invalid_entry)
+               break;
+             else
+               ++coordination;
+
+                                            // insert this dof at its
+                                            // coordination number
+           const std::pair<const unsigned int, int> new_entry (coordination, *s);
+           dofs_by_coordination.insert (new_entry);
+         };
+      
+                                        ////
+       std::multimap<unsigned int, int>::iterator i;
+       for (i = dofs_by_coordination.begin(); i!=dofs_by_coordination.end(); ++i) 
+         new_indices[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
+                                    //  test for all indices numbered
 //TODO: Test fails. Do not check before unification.
-  if (std::find (new_indices.begin(), new_indices.end(),
-                DoFHandler<dim>::invalid_dof_index)
-      !=
-      new_indices.end())
-    Assert (false, ExcRenumberingIncomplete());
-  Assert (next_free_number == n_dofs,
-         ExcRenumberingIncomplete());
+    if (std::find (new_indices.begin(), new_indices.end(),
+                  DoFHandler<dim>::invalid_dof_index)
+       !=
+       new_indices.end())
+      Assert (false, ExcRenumberingIncomplete());
+    Assert (next_free_number == n_dofs,
+           ExcRenumberingIncomplete());
 #endif
 
-  if (reversed_numbering)
-    for (std::vector<unsigned int>::iterator i=new_indices.begin(); i!=new_indices.end(); ++i)
-      *i = n_dofs-*i;
+    if (reversed_numbering)
+      for (std::vector<unsigned int>::iterator i=new_indices.begin(); i!=new_indices.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_indices);
-}
+                                    // actually perform renumbering;
+                                    // this is dimension specific and
+                                    // thus needs an own function
+    dof_handler.renumber_dofs (level, new_indices);
+  }
 
 
-template <int dim>
-void
-DoFRenumbering::
-component_wise (DoFHandler<dim>                 &dof_handler,
-                const std::vector<unsigned int> &component_order_arg)
-{
-  std::vector<unsigned int> renumbering (dof_handler.n_dofs(),
-                                        DoFHandler<dim>::invalid_dof_index);
-
-  typedef
-    WrapDoFIterator<typename DoFHandler<dim>::active_cell_iterator> ITERATOR;
-
-  typename DoFHandler<dim>::active_cell_iterator
-    istart = dof_handler.begin_active();
-  ITERATOR start = istart;
-  const typename DoFHandler<dim>::cell_iterator
-    end = dof_handler.end();
-
-  unsigned int result =
-  compute_component_wise<dim, ITERATOR,
-    typename DoFHandler<dim>::cell_iterator>(renumbering,
-                                            start, end,
-                                            component_order_arg);
+  template <int dim>
+  void
+
+  component_wise (DoFHandler<dim>                 &dof_handler,
+                 const std::vector<unsigned int> &component_order_arg)
+  {
+    std::vector<unsigned int> renumbering (dof_handler.n_dofs(),
+                                          DoFHandler<dim>::invalid_dof_index);
+
+    typedef
+      internal::WrapDoFIterator<typename DoFHandler<dim>::active_cell_iterator> ITERATOR;
+
+    typename DoFHandler<dim>::active_cell_iterator
+      istart = dof_handler.begin_active();
+    ITERATOR start = istart;
+    const typename DoFHandler<dim>::cell_iterator
+      end = dof_handler.end();
 
-  if (result == 0) return;
+    unsigned int result =
+      compute_component_wise<dim, ITERATOR,
+      typename DoFHandler<dim>::cell_iterator>(renumbering,
+                                              start, end,
+                                              component_order_arg);
+
+    if (result == 0) return;
   
-  Assert (result == dof_handler.n_dofs(),
-         ExcRenumberingIncomplete());
+    Assert (result == dof_handler.n_dofs(),
+           ExcRenumberingIncomplete());
   
-  dof_handler.renumber_dofs (renumbering);
-}
+    dof_handler.renumber_dofs (renumbering);
+  }
 
 
 
-template <int dim>
-void DoFRenumbering::
-component_wise (MGDoFHandler<dim> &dof_handler,
-                const unsigned int level,
-                const std::vector<unsigned int> &component_order_arg)
-{
-  std::vector<unsigned int> renumbering (dof_handler.n_dofs(level),
-                                        DoFHandler<dim>::invalid_dof_index);
+  template <int dim>
+  void 
+  component_wise (MGDoFHandler<dim> &dof_handler,
+                 const unsigned int level,
+                 const std::vector<unsigned int> &component_order_arg)
+  {
+    std::vector<unsigned int> renumbering (dof_handler.n_dofs(level),
+                                          DoFHandler<dim>::invalid_dof_index);
 
-  typedef
-    WrapMGDoFIterator<typename MGDoFHandler<dim>::cell_iterator> ITERATOR;
+    typedef
+      internal::WrapMGDoFIterator<typename MGDoFHandler<dim>::cell_iterator> ITERATOR;
   
-  typename MGDoFHandler<dim>::cell_iterator
-    istart =dof_handler.begin(level);
-  ITERATOR start = istart;
-  typename MGDoFHandler<dim>::cell_iterator
-    iend = dof_handler.end(level);
-  const ITERATOR end = iend;
+    typename MGDoFHandler<dim>::cell_iterator
+      istart =dof_handler.begin(level);
+    ITERATOR start = istart;
+    typename MGDoFHandler<dim>::cell_iterator
+      iend = dof_handler.end(level);
+    const ITERATOR end = iend;
 
-  const unsigned int result =
-    compute_component_wise<dim, ITERATOR, ITERATOR>(
-      renumbering, start, end, component_order_arg);
+    const unsigned int result =
+      compute_component_wise<dim, ITERATOR, ITERATOR>(
+       renumbering, start, end, component_order_arg);
 
-  if (result == 0) return;
+    if (result == 0) return;
   
-  Assert (result == dof_handler.n_dofs(level),
-         ExcRenumberingIncomplete());
+    Assert (result == dof_handler.n_dofs(level),
+           ExcRenumberingIncomplete());
   
-  if (renumbering.size()!=0)
-    dof_handler.renumber_dofs (level, renumbering);
-}
+    if (renumbering.size()!=0)
+      dof_handler.renumber_dofs (level, renumbering);
+  }
 
 
 
-template <int dim>
-void DoFRenumbering::
-component_wise (MGDoFHandler<dim> &dof_handler,
-                const std::vector<unsigned int> &component_order_arg)
-{
-  Threads::ThreadGroup<> threads;
+  template <int dim>
+  void 
+  component_wise (MGDoFHandler<dim> &dof_handler,
+                 const std::vector<unsigned int> &component_order_arg)
+  {
+    Threads::ThreadGroup<> threads;
 
-  void (*non_mg_part) (DoFHandler<dim> &, const std::vector<unsigned int> &)
-    = &DoFRenumbering::component_wise<dim>;
-  void (*mg_part) (MGDoFHandler<dim> &, unsigned int, const std::vector<unsigned int> &)
-    = &DoFRenumbering::component_wise<dim>;
+    void (*non_mg_part) (DoFHandler<dim> &, const std::vector<unsigned int> &)
+      = &component_wise<dim>;
+    void (*mg_part) (MGDoFHandler<dim> &, unsigned int, const std::vector<unsigned int> &)
+      = &component_wise<dim>;
   
-  threads += Threads::spawn (non_mg_part) (static_cast<DoFHandler<dim>&> (dof_handler),
-                                          component_order_arg);
-  for (unsigned int level=0; level<dof_handler.get_tria().n_levels(); ++level)
-    threads += Threads::spawn (mg_part) (dof_handler, level, component_order_arg);
+    threads += Threads::spawn (non_mg_part) (static_cast<DoFHandler<dim>&> (dof_handler),
+                                            component_order_arg);
+    for (unsigned int level=0; level<dof_handler.get_tria().n_levels(); ++level)
+      threads += Threads::spawn (mg_part) (dof_handler, level, component_order_arg);
 
-  threads.join_all();
-}
+    threads.join_all();
+  }
 
 
 
 
-template <int dim, class ITERATOR, class ENDITERATOR>
-unsigned int
-DoFRenumbering::
-compute_component_wise (std::vector<unsigned int>& new_indices,
-                        ITERATOR& start,
-                        const ENDITERATOR& end,
-                        const std::vector<unsigned int> &component_order_arg)
-{
+  template <int dim, class ITERATOR, class ENDITERATOR>
+  unsigned int
+  compute_component_wise (std::vector<unsigned int>& new_indices,
+                         ITERATOR& start,
+                         const ENDITERATOR& end,
+                         const std::vector<unsigned int> &component_order_arg)
+  {
 //TODO: Modify for hp
-  const FiniteElement<dim>& fe = start->get_fe();
-  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+    const FiniteElement<dim>& fe = start->get_fe();
+    const unsigned int dofs_per_cell = fe.dofs_per_cell;
 
 //    Assert (new_indices.size() ==  start->dof_handler().n_dofs(),
 //       ExcDimensionMismatch(new_indices.size(),
 //                            start->dof_handler.n_dofs()));
 
-                                  // do nothing if the FE has only
-                                  // one component
-  if (fe.n_components() == 1)
-    {
-      new_indices.resize(0);
-      return 0;
-    }
+                                    // do nothing if the FE has only
+                                    // one component
+    if (fe.n_components() == 1)
+      {
+       new_indices.resize(0);
+       return 0;
+      }
 
-                                  // Copy last argument into a
-                                  // writable vector.
-  std::vector<unsigned int> component_order (component_order_arg);
-                                  // If the last argument was an
-                                  // empty vector, set up things to
-                                  // store components in the order
-                                  // found in the system.
-  if (component_order.size() == 0)
-    for (unsigned int i=0; i<fe.n_components(); ++i)
-      component_order.push_back (i);
-
-  Assert (component_order.size() == fe.n_components(),
-         ExcDimensionMismatch(component_order.size(), fe.n_components()));
-  
-  for (unsigned int i=0; i< component_order.size(); ++i)
-    Assert(component_order[i] < fe.n_components(),
-          ExcIndexRange(component_order[i], 0, fe.n_components()));
-  
-                                  // vector to hold the dof indices on
-                                  // the cell we visit at a time
-  std::vector<unsigned int> local_dof_indices(dofs_per_cell);
-
-                                   // prebuilt list to which component
-                                  // a given dof on a cell
-                                  // should go. note that we get into
-                                  // trouble here if the shape
-                                  // function is not primitive, since
-                                  // then there is no single vector
-                                  // component to which it
-                                  // belongs. in this case, assign it
-                                  // to the first vector component to
-                                  // which it belongs
-  std::vector<unsigned int> component_list (dofs_per_cell);
-  for (unsigned int i=0; i<component_list.size(); ++i)
-    if (fe.is_primitive(i))
-      component_list[i]
-       = component_order[fe.system_to_component_index(i).first];
-    else
+                                    // Copy last argument into a
+                                    // writable vector.
+    std::vector<unsigned int> component_order (component_order_arg);
+                                    // If the last argument was an
+                                    // empty vector, set up things to
+                                    // store components in the order
+                                    // found in the system.
+    if (component_order.size() == 0)
+      for (unsigned int i=0; i<fe.n_components(); ++i)
+       component_order.push_back (i);
+
+    Assert (component_order.size() == fe.n_components(),
+           ExcDimensionMismatch(component_order.size(), fe.n_components()));
+  
+    for (unsigned int i=0; i< component_order.size(); ++i)
+      Assert(component_order[i] < fe.n_components(),
+            ExcIndexRange(component_order[i], 0, fe.n_components()));
+  
+                                    // vector to hold the dof indices on
+                                    // the cell we visit at a time
+    std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+
+                                    // prebuilt list to which component
+                                    // a given dof on a cell
+                                    // should go. note that we get into
+                                    // trouble here if the shape
+                                    // function is not primitive, since
+                                    // then there is no single vector
+                                    // component to which it
+                                    // belongs. in this case, assign it
+                                    // to the first vector component to
+                                    // which it belongs
+    std::vector<unsigned int> component_list (dofs_per_cell);
+    for (unsigned int i=0; i<component_list.size(); ++i)
+      if (fe.is_primitive(i))
+       component_list[i]
+         = component_order[fe.system_to_component_index(i).first];
+      else
+       {
+         const unsigned int base_element =
+           fe.system_to_base_index(i).first.first;
+         const unsigned int base_multiplicity =
+           fe.system_to_base_index(i).first.second;
+                                          // sum up the number of
+                                          // components all the
+                                          // elements before that have
+         unsigned int c=0;
+         for (unsigned int b=0; b<base_element; ++b)
+           c += fe.base_element(b).n_components() *
+                fe.element_multiplicity(b);
+         for (unsigned int m=0; m<base_multiplicity; ++m)
+           c += fe.base_element(base_element).n_components();
+                                          // then associate this degree
+                                          // of freedom with this
+                                          // component
+         component_list[i] = component_order[c];
+       };
+  
+                                    // set up a map where for each
+                                    // component the respective degrees
+                                    // of freedom are collected.
+                                    //
+                                    // note that this map is sorted by
+                                    // component but that within each
+                                    // component it is NOT sorted by
+                                    // dof index. note also that some
+                                    // dof indices are entered
+                                    // multiply, so we will have to
+                                    // take care of that
+    std::vector<std::vector<unsigned int> >
+      component_to_dof_map (fe.n_components());
+    for (;start!=end;++start)
       {
-        const unsigned int base_element =
-          fe.system_to_base_index(i).first.first;
-        const unsigned int base_multiplicity =
-          fe.system_to_base_index(i).first.second;
-                                         // sum up the number of
-                                         // components all the
-                                         // elements before that have
-        unsigned int c=0;
-        for (unsigned int b=0; b<base_element; ++b)
-          c += fe.base_element(b).n_components() *
-               fe.element_multiplicity(b);
-        for (unsigned int m=0; m<base_multiplicity; ++m)
-          c += fe.base_element(base_element).n_components();
-                                         // then associate this degree
-                                         // of freedom with this
-                                         // component
-        component_list[i] = component_order[c];
+                                        // on each cell: get dof indices
+                                        // and insert them into the global
+                                        // list using their component
+       start.get_dof_indices (local_dof_indices);
+       for (unsigned int i=0; i<dofs_per_cell; ++i)
+         component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
       };
   
-                                  // set up a map where for each
-                                  // component the respective degrees
-                                  // of freedom are collected.
-                                  //
-                                  // note that this map is sorted by
-                                  // component but that within each
-                                  // component it is NOT sorted by
-                                  // dof index. note also that some
-                                  // dof indices are entered
-                                  // multiply, so we will have to
-                                  // take care of that
-  std::vector<std::vector<unsigned int> >
-    component_to_dof_map (fe.n_components());
-  for (;start!=end;++start)
-    {
-                                      // on each cell: get dof indices
-                                      // and insert them into the global
-                                      // list using their component
-      start.get_dof_indices (local_dof_indices);
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
-    };
-  
-                                  // now we've got all indices sorted
-                                  // into buckets labelled with their
-                                  // target component number. we've
-                                  // only got to traverse this list
-                                  // and assign the new indices
-                                  //
-                                  // however, we first want to sort
-                                  // the indices entered into the
-                                  // buckets to preserve the order
-                                  // within each component and during
-                                  // this also remove duplicate
-                                  // entries
-                                  //
-                                  // note that we no longer have to
-                                  // care about non-primitive shape
-                                  // functions since the buckets
-                                  // corresponding to the second and
-                                  // following vector components of a
-                                  // non-primitive FE will simply be
-                                  // empty, everything being shoved
-                                  // into the first one. The same
-                                  // holds if several components were
-                                  // joined into a single target.
-  for (unsigned int component=0; component<fe.n_components(); ++component)
-    {
-      std::sort (component_to_dof_map[component].begin(),
-                component_to_dof_map[component].end());
-      component_to_dof_map[component]
-        .erase (std::unique (component_to_dof_map[component].begin(),
-                             component_to_dof_map[component].end()),
-                component_to_dof_map[component].end());
-    };
+                                    // now we've got all indices sorted
+                                    // into buckets labelled with their
+                                    // target component number. we've
+                                    // only got to traverse this list
+                                    // and assign the new indices
+                                    //
+                                    // however, we first want to sort
+                                    // the indices entered into the
+                                    // buckets to preserve the order
+                                    // within each component and during
+                                    // this also remove duplicate
+                                    // entries
+                                    //
+                                    // note that we no longer have to
+                                    // care about non-primitive shape
+                                    // functions since the buckets
+                                    // corresponding to the second and
+                                    // following vector components of a
+                                    // non-primitive FE will simply be
+                                    // empty, everything being shoved
+                                    // into the first one. The same
+                                    // holds if several components were
+                                    // joined into a single target.
+    for (unsigned int component=0; component<fe.n_components(); ++component)
+      {
+       std::sort (component_to_dof_map[component].begin(),
+                  component_to_dof_map[component].end());
+       component_to_dof_map[component]
+         .erase (std::unique (component_to_dof_map[component].begin(),
+                              component_to_dof_map[component].end()),
+                 component_to_dof_map[component].end());
+      };
 
-                                   // now concatenate all the
-                                   // components in the order the user
-                                   // desired to see
-  unsigned int next_free_index = 0;
-  for (unsigned int c=0; c<fe.n_components(); ++c)
-    {
-      const unsigned int component = c;
+                                    // now concatenate all the
+                                    // components in the order the user
+                                    // desired to see
+    unsigned int next_free_index = 0;
+    for (unsigned int c=0; c<fe.n_components(); ++c)
+      {
+       const unsigned int component = c;
       
-      const typename std::vector<unsigned int>::const_iterator
-       begin_of_component = component_to_dof_map[component].begin(),
-       end_of_component   = component_to_dof_map[component].end();
+       const typename std::vector<unsigned int>::const_iterator
+         begin_of_component = component_to_dof_map[component].begin(),
+         end_of_component   = component_to_dof_map[component].end();
             
-      for (typename std::vector<unsigned int>::const_iterator
-             dof_index = begin_of_component;
-          dof_index != end_of_component; ++dof_index)
-       new_indices[*dof_index] = next_free_index++;
-    };
+       for (typename std::vector<unsigned int>::const_iterator
+              dof_index = begin_of_component;
+            dof_index != end_of_component; ++dof_index)
+         new_indices[*dof_index] = next_free_index++;
+      };
 
-  return next_free_index;
-}
+    return next_free_index;
+  }
 
 
 
-template <class DH>
-void
-DoFRenumbering::sort_selected_dofs_back (
-  DH&                      dof_handler,
-  const std::vector<bool>& selected_dofs)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
+  template <class DH>
+  void
+  sort_selected_dofs_back (DH&                      dof_handler,
+                          const std::vector<bool>& selected_dofs)
+  {
+    std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                         DH::invalid_dof_index);
+    compute_sort_selected_dofs_back(renumbering, dof_handler, selected_dofs);
 
-  dof_handler.renumber_dofs(renumbering);
-}
+    dof_handler.renumber_dofs(renumbering);
+  }
 
 
 
-template <class DH>
-void
-DoFRenumbering::compute_sort_selected_dofs_back (
-  std::vector<unsigned int>& new_indices,
-  const DH&                  dof_handler,
-  const std::vector<bool>&   selected_dofs)
-{
-  const unsigned int n_dofs = dof_handler.n_dofs();
-  Assert (selected_dofs.size() == n_dofs,
-         ExcDimensionMismatch (selected_dofs.size(), n_dofs));
-
-                                  // re-sort the dofs according to
-                                  // their selection state
-  Assert (new_indices.size() == n_dofs,
-         ExcDimensionMismatch(new_indices.size(), n_dofs));
-  
-  const unsigned int   n_selected_dofs = std::count (selected_dofs.begin(),
-                                                    selected_dofs.end(),
-                                                    false);
-  
-  unsigned int next_unselected = 0;
-  unsigned int next_selected   = n_selected_dofs;
-  for (unsigned int i=0; i<n_dofs; ++i)
-    if (selected_dofs[i] == false)
-      {
-       new_indices[i] = next_unselected;
-       ++next_unselected;
-      }
-    else
-      {
-       new_indices[i] = next_selected;
-       ++next_selected;
-      };
-  Assert (next_unselected == n_selected_dofs, ExcInternalError());
-  Assert (next_selected == n_dofs, ExcInternalError());
-}
+  template <class DH>
+  void
+  compute_sort_selected_dofs_back (std::vector<unsigned int>& new_indices,
+                                  const DH&                  dof_handler,
+                                  const std::vector<bool>&   selected_dofs)
+  {
+    const unsigned int n_dofs = dof_handler.n_dofs();
+    Assert (selected_dofs.size() == n_dofs,
+           ExcDimensionMismatch (selected_dofs.size(), n_dofs));
 
+                                    // re-sort the dofs according to
+                                    // their selection state
+    Assert (new_indices.size() == n_dofs,
+           ExcDimensionMismatch(new_indices.size(), n_dofs));
+  
+    const unsigned int   n_selected_dofs = std::count (selected_dofs.begin(),
+                                                      selected_dofs.end(),
+                                                      false);
+  
+    unsigned int next_unselected = 0;
+    unsigned int next_selected   = n_selected_dofs;
+    for (unsigned int i=0; i<n_dofs; ++i)
+      if (selected_dofs[i] == false)
+       {
+         new_indices[i] = next_unselected;
+         ++next_unselected;
+       }
+      else
+       {
+         new_indices[i] = next_selected;
+         ++next_selected;
+       };
+    Assert (next_unselected == n_selected_dofs, ExcInternalError());
+    Assert (next_selected == n_dofs, ExcInternalError());
+  }
 
 
-template <class DH>
-void
-DoFRenumbering::cell_wise_dg (
-  DH& dof,
-  const std::vector<typename DH::cell_iterator>& cells)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs());
-  std::vector<unsigned int> reverse(dof.n_dofs());
-  compute_cell_wise_dg(renumbering, reverse, dof, cells);
+
+  template <class DH>
+  void
+  cell_wise_dg (DH& dof,
+               const std::vector<typename DH::cell_iterator>& cells)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs());
+    std::vector<unsigned int> reverse(dof.n_dofs());
+    compute_cell_wise_dg(renumbering, reverse, dof, cells);
   
-  dof.renumber_dofs(renumbering);
-}
+    dof.renumber_dofs(renumbering);
+  }
 
+  
 
-template <class DH>
-void
-DoFRenumbering::cell_wise (
-  DH& dof,
-  const std::vector<typename DH::cell_iterator>& cells)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs());
-  std::vector<unsigned int> reverse(dof.n_dofs());
-  compute_cell_wise(renumbering, reverse, dof, cells);
+  template <class DH>
+  void
+  cell_wise (
+    DH& dof,
+    const std::vector<typename DH::cell_iterator>& cells)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs());
+    std::vector<unsigned int> reverse(dof.n_dofs());
+    compute_cell_wise(renumbering, reverse, dof, cells);
   
-  dof.renumber_dofs(renumbering);
-}
+    dof.renumber_dofs(renumbering);
+  }
 
 
 //TODO: Discuss if we cannot replace this function by the next
-template <class DH>
-void
-DoFRenumbering::compute_cell_wise_dg (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const DH& dof,
-  const typename std::vector<typename DH::cell_iterator>& cells)
-{
-  Assert(cells.size() == dof.get_tria().n_active_cells(),
-        ExcDimensionMismatch(cells.size(),
-                             dof.get_tria().n_active_cells()));
+  template <class DH>
+  void
+  compute_cell_wise_dg (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const DH& dof,
+    const typename std::vector<typename DH::cell_iterator>& cells)
+  {
+    Assert(cells.size() == dof.get_tria().n_active_cells(),
+          ExcDimensionMismatch(cells.size(),
+                               dof.get_tria().n_active_cells()));
 
-  unsigned int n_global_dofs = dof.n_dofs();
+    unsigned int n_global_dofs = dof.n_dofs();
 
-  // Actually, we compute the inverse of the reordering vector, called reverse here.
-  // Later, irs inverse is computed into new_indices, which is the return argument.
+                                    // Actually, we compute the
+                                    // inverse of the reordering
+                                    // vector, called reverse here.
+                                    // Later, irs inverse is computed
+                                    // into new_indices, which is the
+                                    // return argument.
   
-  Assert(new_indices.size() == n_global_dofs,
-        ExcDimensionMismatch(new_indices.size(), n_global_dofs));
-  Assert(reverse.size() == n_global_dofs,
-        ExcDimensionMismatch(reverse.size(), n_global_dofs));
+    Assert(new_indices.size() == n_global_dofs,
+          ExcDimensionMismatch(new_indices.size(), n_global_dofs));
+    Assert(reverse.size() == n_global_dofs,
+          ExcDimensionMismatch(reverse.size(), n_global_dofs));
   
-  std::vector<unsigned int> cell_dofs;
+    std::vector<unsigned int> cell_dofs;
 
-  unsigned int global_index = 0;
+    unsigned int global_index = 0;
   
-  typename std::vector<typename DH::cell_iterator>::const_iterator cell;
+    typename std::vector<typename DH::cell_iterator>::const_iterator cell;
 
-  for(cell = cells.begin(); cell != cells.end(); ++cell)
-    {
-      Assert((*cell)->get_fe().n_dofs_per_face()==0, ExcNotDGFEM());
-                                      // Determine the number of dofs
-                                      // on this cell and reinit the
-                                      // vector storing these
-                                      // numbers.
-      unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
-      cell_dofs.resize(n_cell_dofs);
+    for(cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+       Assert((*cell)->get_fe().n_dofs_per_face()==0, ExcNotDGFEM());
+                                        // Determine the number of dofs
+                                        // on this cell and reinit the
+                                        // vector storing these
+                                        // numbers.
+       unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
+       cell_dofs.resize(n_cell_dofs);
       
-      (*cell)->get_dof_indices(cell_dofs);
-
-                                      // Sort here to make sure that
-                                      // degrees of freedom inside a
-                                      // single cell are in the same
-                                      // order after renumbering.
-      std::sort(cell_dofs.begin(), cell_dofs.end());
-
-      for (unsigned int i=0;i<n_cell_dofs;++i)
-       {
-         reverse[global_index++] = cell_dofs[i];
-       }
-    }
-  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
+       (*cell)->get_dof_indices(cell_dofs);
+
+                                        // Sort here to make sure that
+                                        // degrees of freedom inside a
+                                        // single cell are in the same
+                                        // order after renumbering.
+       std::sort(cell_dofs.begin(), cell_dofs.end());
+
+       for (unsigned int i=0;i<n_cell_dofs;++i)
+         {
+           reverse[global_index++] = cell_dofs[i];
+         }
+      }
+    Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
 
-  for (unsigned int i=0;i<reverse.size(); ++i)
-    new_indices[reverse[i]] = i;
-}
+    for (unsigned int i=0;i<reverse.size(); ++i)
+      new_indices[reverse[i]] = i;
+  }
 
 
-template <class DH>
-void
-DoFRenumbering::compute_cell_wise (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const DH& dof,
-  const typename std::vector<typename DH::cell_iterator>& cells)
-{
-  Assert(cells.size() == dof.get_tria().n_active_cells(),
-        ExcDimensionMismatch(cells.size(),
-                             dof.get_tria().n_active_cells()));
+  template <class DH>
+  void
+  compute_cell_wise (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const DH& dof,
+    const typename std::vector<typename DH::cell_iterator>& cells)
+  {
+    Assert(cells.size() == dof.get_tria().n_active_cells(),
+          ExcDimensionMismatch(cells.size(),
+                               dof.get_tria().n_active_cells()));
 
-  unsigned int n_global_dofs = dof.n_dofs();
+    unsigned int n_global_dofs = dof.n_dofs();
   
-  // Actually, we compute the inverse of the reordering vector, called reverse here.
-  // Later, irs inverse is computed into new_indices, which is the return argument.
+                                    // Actually, we compute the
+                                    // inverse of the reordering
+                                    // vector, called reverse here.
+                                    // Later, irs inverse is computed
+                                    // into new_indices, which is the
+                                    // return argument.
   
-  Assert(new_indices.size() == n_global_dofs,
-        ExcDimensionMismatch(new_indices.size(), n_global_dofs));
-  Assert(reverse.size() == n_global_dofs,
-        ExcDimensionMismatch(reverse.size(), n_global_dofs));
+    Assert(new_indices.size() == n_global_dofs,
+          ExcDimensionMismatch(new_indices.size(), n_global_dofs));
+    Assert(reverse.size() == n_global_dofs,
+          ExcDimensionMismatch(reverse.size(), n_global_dofs));
 
-                                  // For continuous elements, we must
-                                  // make sure, that each dof is
-                                  // reordered only once.
-  std::vector<bool> already_sorted(n_global_dofs, false);
-  std::vector<unsigned int> cell_dofs;
+                                    // For continuous elements, we must
+                                    // make sure, that each dof is
+                                    // reordered only once.
+    std::vector<bool> already_sorted(n_global_dofs, false);
+    std::vector<unsigned int> cell_dofs;
 
-  unsigned int global_index = 0;
+    unsigned int global_index = 0;
   
-  typename std::vector<typename DH::cell_iterator>::const_iterator cell;
+    typename std::vector<typename DH::cell_iterator>::const_iterator cell;
 
-  for(cell = cells.begin(); cell != cells.end(); ++cell)
-    {
-                                      // Determine the number of dofs
-                                      // on this cell and reinit the
-                                      // vector storing these
-                                      // numbers.
-      unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
-      cell_dofs.resize(n_cell_dofs);
+    for(cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+                                        // Determine the number of dofs
+                                        // on this cell and reinit the
+                                        // vector storing these
+                                        // numbers.
+       unsigned int n_cell_dofs = (*cell)->get_fe().n_dofs_per_cell();
+       cell_dofs.resize(n_cell_dofs);
       
-      (*cell)->get_dof_indices(cell_dofs);
-
-                                      // Sort here to make sure that
-                                      // degrees of freedom inside a
-                                      // single cell are in the same
-                                      // order after renumbering.
-      std::sort(cell_dofs.begin(), cell_dofs.end());
-
-      for (unsigned int i=0;i<n_cell_dofs;++i)
-       {
-         if (!already_sorted[cell_dofs[i]])
-           {
-             already_sorted[cell_dofs[i]] = true;
-             reverse[global_index++] = cell_dofs[i];
-           }
-       }
-    }
-  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
+       (*cell)->get_dof_indices(cell_dofs);
+
+                                        // Sort here to make sure that
+                                        // degrees of freedom inside a
+                                        // single cell are in the same
+                                        // order after renumbering.
+       std::sort(cell_dofs.begin(), cell_dofs.end());
+
+       for (unsigned int i=0;i<n_cell_dofs;++i)
+         {
+           if (!already_sorted[cell_dofs[i]])
+             {
+               already_sorted[cell_dofs[i]] = true;
+               reverse[global_index++] = cell_dofs[i];
+             }
+         }
+      }
+    Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
 
-  for (unsigned int i=0;i<reverse.size(); ++i)
-    new_indices[reverse[i]] = i;
-}
+    for (unsigned int i=0;i<reverse.size(); ++i)
+      new_indices[reverse[i]] = i;
+  }
 
+  
 
-template <int dim>
-void DoFRenumbering::cell_wise_dg (
-  MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs(level));
-  std::vector<unsigned int> reverse(dof.n_dofs(level));
+  template <int dim>
+  void cell_wise_dg (
+    MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs(level));
+    std::vector<unsigned int> reverse(dof.n_dofs(level));
   
-  compute_cell_wise_dg(renumbering, reverse, dof, level, cells);
-  dof.renumber_dofs(level, renumbering);
-}
+    compute_cell_wise_dg(renumbering, reverse, dof, level, cells);
+    dof.renumber_dofs(level, renumbering);
+  }
 
+  
 
-template <int dim>
-void DoFRenumbering::cell_wise (
-  MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs(level));
-  std::vector<unsigned int> reverse(dof.n_dofs(level));
+  template <int dim>
+  void cell_wise (
+    MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs(level));
+    std::vector<unsigned int> reverse(dof.n_dofs(level));
   
-  compute_cell_wise(renumbering, reverse, dof, level, cells);
-  dof.renumber_dofs(level, renumbering);
-}
+    compute_cell_wise(renumbering, reverse, dof, level, cells);
+    dof.renumber_dofs(level, renumbering);
+  }
 
+  
 
-template <int dim>
-void DoFRenumbering::compute_cell_wise_dg (
-  std::vector<unsigned int>& new_order,
-  std::vector<unsigned int>& reverse,
-  const MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
-{
-  Assert(cells.size() == dof.get_tria().n_cells(level),
-        ExcDimensionMismatch(cells.size(),
-                             dof.get_tria().n_cells(level)));
-  switch (dim)
-    {
-      case 3:
-           Assert(dof.get_fe().n_dofs_per_quad()==0,
-                  ExcNotDGFEM());
-      case 2:
-           Assert(dof.get_fe().n_dofs_per_line()==0,
-                  ExcNotDGFEM());
-      default:
-           Assert(dof.get_fe().n_dofs_per_vertex()==0,
-                  ExcNotDGFEM());
-    }
+  template <int dim>
+  void compute_cell_wise_dg (
+    std::vector<unsigned int>& new_order,
+    std::vector<unsigned int>& reverse,
+    const MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
+  {
+    Assert(cells.size() == dof.get_tria().n_cells(level),
+          ExcDimensionMismatch(cells.size(),
+                               dof.get_tria().n_cells(level)));
+    switch (dim)
+      {
+       case 3:
+             Assert(dof.get_fe().n_dofs_per_quad()==0,
+                    ExcNotDGFEM());
+       case 2:
+             Assert(dof.get_fe().n_dofs_per_line()==0,
+                    ExcNotDGFEM());
+       default:
+             Assert(dof.get_fe().n_dofs_per_vertex()==0,
+                    ExcNotDGFEM());
+      }
   
-  Assert (new_order.size() == dof.n_dofs(level),
-         ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
-  Assert (reverse.size() == dof.n_dofs(level),
-         ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
+    Assert (new_order.size() == dof.n_dofs(level),
+           ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
+    Assert (reverse.size() == dof.n_dofs(level),
+           ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
   
-  unsigned int n_global_dofs = dof.n_dofs(level);
-  unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
+    unsigned int n_global_dofs = dof.n_dofs(level);
+    unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
   
-  std::vector<unsigned int> cell_dofs(n_cell_dofs);
+    std::vector<unsigned int> cell_dofs(n_cell_dofs);
   
-  unsigned int global_index = 0;
+    unsigned int global_index = 0;
   
-  typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
+    typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
 
-  for(cell = cells.begin(); cell != cells.end(); ++cell)
-    {
-      Assert ((*cell)->level() == (int) level, ExcInternalError());
+    for(cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+       Assert ((*cell)->level() == (int) level, ExcInternalError());
 
-      (*cell)->get_mg_dof_indices(cell_dofs);
-      std::sort(cell_dofs.begin(), cell_dofs.end());
+       (*cell)->get_mg_dof_indices(cell_dofs);
+       std::sort(cell_dofs.begin(), cell_dofs.end());
 
-      for (unsigned int i=0;i<n_cell_dofs;++i)
-       {
-         reverse[global_index++] = cell_dofs[i];
-       }
-    }
-  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
+       for (unsigned int i=0;i<n_cell_dofs;++i)
+         {
+           reverse[global_index++] = cell_dofs[i];
+         }
+      }
+    Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
 
-  for (unsigned int i=0;i<new_order.size(); ++i)
-    new_order[reverse[i]] = i;  
-}
+    for (unsigned int i=0;i<new_order.size(); ++i)
+      new_order[reverse[i]] = i;  
+  }
 
 
 
-template <int dim>
-void DoFRenumbering::compute_cell_wise (
-  std::vector<unsigned int>& new_order,
-  std::vector<unsigned int>& reverse,
-  const MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
-{
-  Assert(cells.size() == dof.get_tria().n_cells(level),
-        ExcDimensionMismatch(cells.size(),
-                             dof.get_tria().n_cells(level)));
-  Assert (new_order.size() == dof.n_dofs(level),
-         ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
-  Assert (reverse.size() == dof.n_dofs(level),
-         ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
+  template <int dim>
+  void compute_cell_wise (
+    std::vector<unsigned int>& new_order,
+    std::vector<unsigned int>& reverse,
+    const MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const typename std::vector<typename MGDoFHandler<dim>::cell_iterator>& cells)
+  {
+    Assert(cells.size() == dof.get_tria().n_cells(level),
+          ExcDimensionMismatch(cells.size(),
+                               dof.get_tria().n_cells(level)));
+    Assert (new_order.size() == dof.n_dofs(level),
+           ExcDimensionMismatch(new_order.size(), dof.n_dofs(level)));
+    Assert (reverse.size() == dof.n_dofs(level),
+           ExcDimensionMismatch(reverse.size(), dof.n_dofs(level)));
   
-  unsigned int n_global_dofs = dof.n_dofs(level);
-  unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
+    unsigned int n_global_dofs = dof.n_dofs(level);
+    unsigned int n_cell_dofs = dof.get_fe().n_dofs_per_cell();
   
-  std::vector<bool> already_sorted(n_global_dofs, false);
-  std::vector<unsigned int> cell_dofs(n_cell_dofs);
+    std::vector<bool> already_sorted(n_global_dofs, false);
+    std::vector<unsigned int> cell_dofs(n_cell_dofs);
   
-  unsigned int global_index = 0;
+    unsigned int global_index = 0;
   
-  typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
+    typename std::vector<typename MGDoFHandler<dim>::cell_iterator>::const_iterator cell;
 
-  for(cell = cells.begin(); cell != cells.end(); ++cell)
-    {
-      Assert ((*cell)->level() == (int) level, ExcInternalError());
+    for(cell = cells.begin(); cell != cells.end(); ++cell)
+      {
+       Assert ((*cell)->level() == (int) level, ExcInternalError());
 
-      (*cell)->get_mg_dof_indices(cell_dofs);
-      std::sort(cell_dofs.begin(), cell_dofs.end());
+       (*cell)->get_mg_dof_indices(cell_dofs);
+       std::sort(cell_dofs.begin(), cell_dofs.end());
       
-      for (unsigned int i=0;i<n_cell_dofs;++i)
-       {
-         if (!already_sorted[cell_dofs[i]])
-           {
-             already_sorted[cell_dofs[i]] = true;
-             reverse[global_index++] = cell_dofs[i];
-           }
-       }
-    }
-  Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
+       for (unsigned int i=0;i<n_cell_dofs;++i)
+         {
+           if (!already_sorted[cell_dofs[i]])
+             {
+               already_sorted[cell_dofs[i]] = true;
+               reverse[global_index++] = cell_dofs[i];
+             }
+         }
+      }
+    Assert(global_index == n_global_dofs, ExcRenumberingIncomplete());
 
-  for (unsigned int i=0;i<new_order.size(); ++i)
-    new_order[reverse[i]] = i;  
-}
+    for (unsigned int i=0;i<new_order.size(); ++i)
+      new_order[reverse[i]] = i;  
+  }
 
 
 
-template <class DH, int dim>
-void
-DoFRenumbering::downstream_dg (DH& dof, const Point<dim>& direction)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs());
-  compute_downstream_dg(renumbering, dof, direction);
+  template <class DH, int dim>
+  void
+  downstream_dg (DH& dof, const Point<dim>& direction)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs());
+    compute_downstream_dg(renumbering, dof, direction);
 
-  dof.renumber_dofs(renumbering);
-}
+    dof.renumber_dofs(renumbering);
+  }
 
 
 
-template <class DH, int dim>
-void
-DoFRenumbering::downstream (DH& dof, const Point<dim>& direction)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs());
-  std::vector<unsigned int> reverse(dof.n_dofs());
-  compute_downstream(renumbering, reverse, dof, direction);
+  template <class DH, int dim>
+  void
+  downstream (DH& dof, const Point<dim>& direction)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs());
+    std::vector<unsigned int> reverse(dof.n_dofs());
+    compute_downstream(renumbering, reverse, dof, direction);
 
-  dof.renumber_dofs(renumbering);
-}
+    dof.renumber_dofs(renumbering);
+  }
 
 
 
-template <class DH, int dim>
-void
-DoFRenumbering::compute_downstream_dg (
-  std::vector<unsigned int>& new_indices,
-  const DH& dof,
-  const Point<dim>& direction)
-{
-  std::vector<typename DH::cell_iterator>
-    ordered_cells(dof.get_tria().n_active_cells());
-  const CompareDownstream<dim> comparator(direction);
+  template <class DH, int dim>
+  void
+  compute_downstream_dg (
+    std::vector<unsigned int>& new_indices,
+    const DH& dof,
+    const Point<dim>& direction)
+  {
+    std::vector<typename DH::cell_iterator>
+      ordered_cells(dof.get_tria().n_active_cells());
+    const internal::CompareDownstream<dim> comparator(direction);
   
-  typename DH::active_cell_iterator begin = dof.begin_active();
-  typename DH::active_cell_iterator end = dof.end();
+    typename DH::active_cell_iterator begin = dof.begin_active();
+    typename DH::active_cell_iterator end = dof.end();
   
-  copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
 
-  std::vector<unsigned int> reverse(new_indices.size());
-  compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
-}
+    std::vector<unsigned int> reverse(new_indices.size());
+    compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
+  }
 
+  
 
-template <class DH, int dim>
-void
-DoFRenumbering::compute_downstream_dg (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const DH& dof,
-  const Point<dim>& direction)
-{
-  std::vector<typename DH::cell_iterator>
-    ordered_cells(dof.get_tria().n_active_cells());
-  const CompareDownstream<dim> comparator(direction);
+  template <class DH, int dim>
+  void
+  compute_downstream_dg (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const DH& dof,
+    const Point<dim>& direction)
+  {
+    std::vector<typename DH::cell_iterator>
+      ordered_cells(dof.get_tria().n_active_cells());
+    const internal::CompareDownstream<dim> comparator(direction);
   
-  typename DH::active_cell_iterator begin = dof.begin_active();
-  typename DH::active_cell_iterator end = dof.end();
+    typename DH::active_cell_iterator begin = dof.begin_active();
+    typename DH::active_cell_iterator end = dof.end();
   
-  copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
   
-  compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
-}
+    compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
+  }
 
 
-template <class DH, int dim>
-void
-DoFRenumbering::compute_downstream (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const DH& dof,
-  const Point<dim>& direction)
-{
-  std::vector<typename DH::cell_iterator>
-    ordered_cells(dof.get_tria().n_active_cells());
-  const CompareDownstream<dim> comparator(direction);
+  template <class DH, int dim>
+  void
+  compute_downstream (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const DH& dof,
+    const Point<dim>& direction)
+  {
+    std::vector<typename DH::cell_iterator>
+      ordered_cells(dof.get_tria().n_active_cells());
+    const internal::CompareDownstream<dim> comparator(direction);
   
-  typename DH::active_cell_iterator begin = dof.begin_active();
-  typename DH::active_cell_iterator end = dof.end();
+    typename DH::active_cell_iterator begin = dof.begin_active();
+    typename DH::active_cell_iterator end = dof.end();
   
-  copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
   
-  compute_cell_wise(new_indices, reverse, dof, ordered_cells);
-}
+    compute_cell_wise(new_indices, reverse, dof, ordered_cells);
+  }
 
 
-template <int dim>
-void DoFRenumbering::downstream_dg (MGDoFHandler<dim>& dof,
-                                   const unsigned int level,
-                                   const Point<dim>& direction)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs(level));
-  std::vector<unsigned int> reverse(dof.n_dofs(level));
-  compute_downstream_dg(renumbering, reverse, dof, level, direction);
   
-  dof.renumber_dofs(level, renumbering);
-}
+  template <int dim>
+  void downstream_dg (MGDoFHandler<dim>& dof,
+                     const unsigned int level,
+                     const Point<dim>& direction)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs(level));
+    std::vector<unsigned int> reverse(dof.n_dofs(level));
+    compute_downstream_dg(renumbering, reverse, dof, level, direction);
+  
+    dof.renumber_dofs(level, renumbering);
+  }
 
+  
 
-template <int dim>
-void DoFRenumbering::downstream (MGDoFHandler<dim>& dof,
-                                const unsigned int level,
-                                const Point<dim>& direction)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs(level));
-  std::vector<unsigned int> reverse(dof.n_dofs(level));
-  compute_downstream(renumbering, reverse, dof, level, direction);
+  template <int dim>
+  void downstream (MGDoFHandler<dim>& dof,
+                  const unsigned int level,
+                  const Point<dim>& direction)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs(level));
+    std::vector<unsigned int> reverse(dof.n_dofs(level));
+    compute_downstream(renumbering, reverse, dof, level, direction);
   
-  dof.renumber_dofs(level, renumbering);
-}
+    dof.renumber_dofs(level, renumbering);
+  }
 
+  
 
-template <int dim>
-void
-DoFRenumbering::compute_downstream_dg (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const Point<dim>& direction)
-{
-  std::vector<typename MGDoFHandler<dim>::cell_iterator>
-    ordered_cells(dof.get_tria().n_cells(level));
-  const CompareDownstream<dim> comparator(direction);
+  template <int dim>
+  void
+  compute_downstream_dg (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const Point<dim>& direction)
+  {
+    std::vector<typename MGDoFHandler<dim>::cell_iterator>
+      ordered_cells(dof.get_tria().n_cells(level));
+    const internal::CompareDownstream<dim> comparator(direction);
   
-  typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
-  typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
+    typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
+    typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
   
-  std::copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    std::copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
   
-  compute_cell_wise_dg(new_indices, reverse, dof, level, ordered_cells);
-}
+    compute_cell_wise_dg(new_indices, reverse, dof, level, ordered_cells);
+  }
 
+  
 
-template <int dim>
-void
-DoFRenumbering::compute_downstream (
-  std::vector<unsigned int>& new_indices,
-  std::vector<unsigned int>& reverse,
-  const MGDoFHandler<dim>& dof,
-  const unsigned int level,
-  const Point<dim>& direction)
-{
-  std::vector<typename MGDoFHandler<dim>::cell_iterator>
-    ordered_cells(dof.get_tria().n_cells(level));
-  const CompareDownstream<dim> comparator(direction);
+  template <int dim>
+  void
+  compute_downstream (
+    std::vector<unsigned int>& new_indices,
+    std::vector<unsigned int>& reverse,
+    const MGDoFHandler<dim>& dof,
+    const unsigned int level,
+    const Point<dim>& direction)
+  {
+    std::vector<typename MGDoFHandler<dim>::cell_iterator>
+      ordered_cells(dof.get_tria().n_cells(level));
+    const internal::CompareDownstream<dim> comparator(direction);
   
-  typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
-  typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
+    typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
+    typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
   
-  std::copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    std::copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
   
-  compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
-}
+    compute_cell_wise(new_indices, reverse, dof, level, ordered_cells);
+  }
 
 
 
 /**
  * Provide comparator for DoFCellAccessors
  */
-
-template <int dim>
-struct ClockCells
-{
-                                    /**
-                                     * Center of rotation.
-                                     */
-    const Point<dim>& center;
-                                    /**
-                                     * Revert sorting order.
-                                     */
-    bool counter;
+  namespace internal
+  {
+    template <int dim>
+    struct ClockCells
+    {
+                                        /**
+                                         * Center of rotation.
+                                         */
+       const Point<dim>& center;
+                                        /**
+                                         * Revert sorting order.
+                                         */
+       bool counter;
     
-                                    /**
-                                     * Constructor.
-                                     */
-    ClockCells (const Point<dim>& center, bool counter) :
-                   center(center),
-                   counter(counter)
-      {}
-                                    /**
-                                     * Comparison operator
-                                     */
-    template <class DHCellIterator>
-    bool operator () (const DHCellIterator& c1,
-                     const DHCellIterator& c2) const
-      {
+                                        /**
+                                         * Constructor.
+                                         */
+       ClockCells (const Point<dim>& center, bool counter) :
+                       center(center),
+                       counter(counter)
+         {}
+                                        /**
+                                         * Comparison operator
+                                         */
+       template <class DHCellIterator>
+       bool operator () (const DHCellIterator& c1,
+                         const DHCellIterator& c2) const
+         {
        
-       const Point<dim> v1 = c1->center() - center;
-       const Point<dim> v2 = c2->center() - center;
-       const double s1 = std::atan2(v1(0), v1(1));
-       const double s2 = std::atan2(v2(0), v2(1));
-       return ( counter ? (s1>s2) : (s2>s1));
-      }
-};
+           const Point<dim> v1 = c1->center() - center;
+           const Point<dim> v2 = c2->center() - center;
+           const double s1 = std::atan2(v1(0), v1(1));
+           const double s2 = std::atan2(v2(0), v2(1));
+           return ( counter ? (s1>s2) : (s2>s1));
+         }
+    };
+  }
 
+  
 
-template <class DH, int dim>
-void
-DoFRenumbering::clockwise_dg (
-  DH& dof,
-  const Point<dim>& center,
-  const bool counter)
-{
-  std::vector<unsigned int> renumbering(dof.n_dofs());
-  compute_clockwise_dg(renumbering, dof, center, counter);
+  template <class DH, int dim>
+  void
+  clockwise_dg (
+    DH& dof,
+    const Point<dim>& center,
+    const bool counter)
+  {
+    std::vector<unsigned int> renumbering(dof.n_dofs());
+    compute_clockwise_dg(renumbering, dof, center, counter);
 
-  dof.renumber_dofs(renumbering);
-}
+    dof.renumber_dofs(renumbering);
+  }
 
 
 
-template <class DH, int dim>
-void
-DoFRenumbering::compute_clockwise_dg (
-  std::vector<unsigned int>& new_indices,
-  const DH& dof,
-  const Point<dim>& center,
-  const bool counter)
-{
-  std::vector<typename DH::cell_iterator>
-    ordered_cells(dof.get_tria().n_active_cells());
-  ClockCells<dim> comparator(center, counter);
+  template <class DH, int dim>
+  void
+  compute_clockwise_dg (
+    std::vector<unsigned int>& new_indices,
+    const DH& dof,
+    const Point<dim>& center,
+    const bool counter)
+  {
+    std::vector<typename DH::cell_iterator>
+      ordered_cells(dof.get_tria().n_active_cells());
+    internal::ClockCells<dim> comparator(center, counter);
   
-  typename DH::active_cell_iterator begin = dof.begin_active();
-  typename DH::active_cell_iterator end = dof.end();
+    typename DH::active_cell_iterator begin = dof.begin_active();
+    typename DH::active_cell_iterator end = dof.end();
   
-  std::copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    std::copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
 
-  std::vector<unsigned int> reverse(new_indices.size());
-  compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
-}
+    std::vector<unsigned int> reverse(new_indices.size());
+    compute_cell_wise_dg(new_indices, reverse, dof, ordered_cells);
+  }
 
+  
 
-template <int dim>
-void DoFRenumbering::clockwise_dg (MGDoFHandler<dim>& dof,
-                                  const unsigned int level,
-                                  const Point<dim>& center,
-                                  const bool counter)
-{
-  std::vector<typename MGDoFHandler<dim>::cell_iterator>
-    ordered_cells(dof.get_tria().n_cells(level));
-  ClockCells<dim> comparator(center, counter);
+  template <int dim>
+  void clockwise_dg (MGDoFHandler<dim>& dof,
+                    const unsigned int level,
+                    const Point<dim>& center,
+                    const bool counter)
+  {
+    std::vector<typename MGDoFHandler<dim>::cell_iterator>
+      ordered_cells(dof.get_tria().n_cells(level));
+    internal::ClockCells<dim> comparator(center, counter);
   
-  typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
-  typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
+    typename MGDoFHandler<dim>::cell_iterator begin = dof.begin(level);
+    typename MGDoFHandler<dim>::cell_iterator end = dof.end(level);
   
-  std::copy (begin, end, ordered_cells.begin());
-  std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
+    std::copy (begin, end, ordered_cells.begin());
+    std::sort (ordered_cells.begin(), ordered_cells.end(), comparator);
 
-  cell_wise_dg(dof, level, ordered_cells);
-}
+    cell_wise_dg(dof, level, ordered_cells);
+  }
 
 
 
-template <class DH>
-void
-DoFRenumbering::random (DH& dof_handler)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DH::invalid_dof_index);
-  compute_random(renumbering, dof_handler);
+  template <class DH>
+  void
+  random (DH& dof_handler)
+  {
+    std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                         DH::invalid_dof_index);
+    compute_random(renumbering, dof_handler);
 
-  dof_handler.renumber_dofs(renumbering);
-}
+    dof_handler.renumber_dofs(renumbering);
+  }
 
+  
 
-template <class DH>
-void
-DoFRenumbering::compute_random (
-  std::vector<unsigned int>& new_indices,
-  const DH&                  dof_handler)
-{
-  const unsigned int n_dofs = dof_handler.n_dofs();
-  Assert(new_indices.size() == n_dofs,
-        ExcDimensionMismatch(new_indices.size(), n_dofs));
+  template <class DH>
+  void
+  compute_random (
+    std::vector<unsigned int>& new_indices,
+    const DH&                  dof_handler)
+  {
+    const unsigned int n_dofs = dof_handler.n_dofs();
+    Assert(new_indices.size() == n_dofs,
+          ExcDimensionMismatch(new_indices.size(), n_dofs));
   
-  for (unsigned i=0; i<n_dofs; ++i)
-    new_indices[i] = i;
+    for (unsigned i=0; i<n_dofs; ++i)
+      new_indices[i] = i;
   
-  std::random_shuffle (new_indices.begin(), new_indices.end());
-}
+    std::random_shuffle (new_indices.begin(), new_indices.end());
+  }
 
 
 
-template <int dim>
-void
-DoFRenumbering::subdomain_wise (DoFHandler<dim> &dof_handler)
-{
-  std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
-                                       DoFHandler<dim>::invalid_dof_index);
-  compute_subdomain_wise(renumbering, dof_handler);
+  template <int dim>
+  void
+  subdomain_wise (DoFHandler<dim> &dof_handler)
+  {
+    std::vector<unsigned int> renumbering(dof_handler.n_dofs(),
+                                         DoFHandler<dim>::invalid_dof_index);
+    compute_subdomain_wise(renumbering, dof_handler);
 
-  dof_handler.renumber_dofs(renumbering);
-}
+    dof_handler.renumber_dofs(renumbering);
+  }
 
+  
 
-template <int dim>
-void
-DoFRenumbering::
-compute_subdomain_wise (std::vector<unsigned int> &new_dof_indices,
-                        const DoFHandler<dim>     &dof_handler)
-{
-  const unsigned int n_dofs = dof_handler.n_dofs();
-  Assert (new_dof_indices.size() == n_dofs,
-         ExcDimensionMismatch (new_dof_indices.size(), n_dofs));
-
-                                   // first get the association of each dof
-                                   // with a subdomain and determine the total
-                                   // number of subdomain ids used
-  std::vector<unsigned int> subdomain_association (n_dofs);
-  DoFTools::get_subdomain_association (dof_handler,
-                                       subdomain_association);
-  const unsigned int n_subdomains
-    = *std::max_element (subdomain_association.begin(),
-                         subdomain_association.end()) + 1;
-  
-                                   // then renumber the subdomains by first
-                                   // looking at those belonging to subdomain
-                                   // 0, then those of subdomain 1, etc. note
-                                   // that the algorithm is stable, i.e. if
-                                   // two dofs i,j have i<j and belong to the
-                                   // same subdomain, then they will be in
-                                   // this order also after reordering
-  std::fill (new_dof_indices.begin(), new_dof_indices.end(),
-             numbers::invalid_unsigned_int);
-  unsigned int next_free_index = 0;
-  for (unsigned int subdomain=0; subdomain<n_subdomains; ++subdomain)
-    for (unsigned int i=0; i<n_dofs; ++i)
-      if (subdomain_association[i] == subdomain)
-        {
-          Assert (new_dof_indices[i] == numbers::invalid_unsigned_int,
-                  ExcInternalError());
-          new_dof_indices[i] = next_free_index;
-          ++next_free_index;
-        }
-
-                                   // we should have numbered all dofs
-  Assert (next_free_index == n_dofs, ExcInternalError());
-  Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
-                     numbers::invalid_unsigned_int)
-          == new_dof_indices.end(),
-          ExcInternalError());
-}
+  template <int dim>
+  void
+  compute_subdomain_wise (std::vector<unsigned int> &new_dof_indices,
+                         const DoFHandler<dim>     &dof_handler)
+  {
+    const unsigned int n_dofs = dof_handler.n_dofs();
+    Assert (new_dof_indices.size() == n_dofs,
+           ExcDimensionMismatch (new_dof_indices.size(), n_dofs));
+
+                                    // first get the association of each dof
+                                    // with a subdomain and determine the total
+                                    // number of subdomain ids used
+    std::vector<unsigned int> subdomain_association (n_dofs);
+    DoFTools::get_subdomain_association (dof_handler,
+                                        subdomain_association);
+    const unsigned int n_subdomains
+      = *std::max_element (subdomain_association.begin(),
+                          subdomain_association.end()) + 1;
+  
+                                    // then renumber the subdomains by first
+                                    // looking at those belonging to subdomain
+                                    // 0, then those of subdomain 1, etc. note
+                                    // that the algorithm is stable, i.e. if
+                                    // two dofs i,j have i<j and belong to the
+                                    // same subdomain, then they will be in
+                                    // this order also after reordering
+    std::fill (new_dof_indices.begin(), new_dof_indices.end(),
+              numbers::invalid_unsigned_int);
+    unsigned int next_free_index = 0;
+    for (unsigned int subdomain=0; subdomain<n_subdomains; ++subdomain)
+      for (unsigned int i=0; i<n_dofs; ++i)
+       if (subdomain_association[i] == subdomain)
+         {
+           Assert (new_dof_indices[i] == numbers::invalid_unsigned_int,
+                   ExcInternalError());
+           new_dof_indices[i] = next_free_index;
+           ++next_free_index;
+         }
+
+                                    // we should have numbered all dofs
+    Assert (next_free_index == n_dofs, ExcInternalError());
+    Assert (std::find (new_dof_indices.begin(), new_dof_indices.end(),
+                      numbers::invalid_unsigned_int)
+           == new_dof_indices.end(),
+           ExcInternalError());
+  }
 
 
 
 
 // explicit instantiations
-//TODO[WB]: also implement the following boost functions for hp DoFHandlers etc.
-template
-void
-DoFRenumbering::boost::Cuthill_McKee (DoFHandler<deal_II_dimension> &, bool, bool);
-
-template
-void
-DoFRenumbering::boost::compute_Cuthill_McKee (std::vector<unsigned int> &,
-                                             const DoFHandler<deal_II_dimension> &, bool, bool);
-
-template
-void
-DoFRenumbering::boost::king_ordering (DoFHandler<deal_II_dimension> &, bool, bool);
-
-template
-void
-DoFRenumbering::boost::compute_king_ordering (std::vector<unsigned int> &,
-                                             const DoFHandler<deal_II_dimension> &, bool, bool);
-
-template
-void
-DoFRenumbering::boost::minimum_degree (DoFHandler<deal_II_dimension> &, bool, bool);
-
-template
-void
-DoFRenumbering::boost::compute_minimum_degree (std::vector<unsigned int> &,
-                                              const DoFHandler<deal_II_dimension> &, bool, bool);
-
+  namespace boost
+  {
+//TODO[WB]: also implement the following boost for hp DoFHandlers etc.
+    template
+    void
+    Cuthill_McKee (DoFHandler<deal_II_dimension> &, bool, bool);
+
+    template
+    void
+    compute_Cuthill_McKee (std::vector<unsigned int> &,
+                          const DoFHandler<deal_II_dimension> &, bool, bool);
+
+    template
+    void
+    king_ordering (DoFHandler<deal_II_dimension> &, bool, bool);
+
+    template
+    void
+    compute_king_ordering (std::vector<unsigned int> &,
+                          const DoFHandler<deal_II_dimension> &, bool, bool);
+
+    template
+    void
+    minimum_degree (DoFHandler<deal_II_dimension> &, bool, bool);
+
+    template
+    void
+    compute_minimum_degree (std::vector<unsigned int> &,
+                           const DoFHandler<deal_II_dimension> &, bool, bool);
+  }
+  
 
 // non-boost functions:
 
-template
-void DoFRenumbering::Cuthill_McKee<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&,
- const bool,
- const bool,
- const std::vector<unsigned int>&);
-
-template
-void
-DoFRenumbering::compute_Cuthill_McKee<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&,
- const bool,
- const bool,
- const std::vector<unsigned int>&);
-
-
-template
-void DoFRenumbering::Cuthill_McKee<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const bool,
- const bool,
- const std::vector<unsigned int>&);
-
-template
-void
-DoFRenumbering::compute_Cuthill_McKee<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const bool,
- const bool,
- const std::vector<unsigned int>&);
-
-
-template
-void DoFRenumbering::component_wise<deal_II_dimension>
-(DoFHandler<deal_II_dimension>&,
- const std::vector<unsigned int>&);
+  template
+  void Cuthill_McKee<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&,
  const bool,
  const bool,
  const std::vector<unsigned int>&);
+
+  template
+  void
+  compute_Cuthill_McKee<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&,
  const bool,
  const bool,
  const std::vector<unsigned int>&);
+
+
+  template
+  void Cuthill_McKee<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const bool,
  const bool,
  const std::vector<unsigned int>&);
+
+  template
+  void
+  compute_Cuthill_McKee<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const bool,
  const bool,
  const std::vector<unsigned int>&);
+
+
+  template
+  void component_wise<deal_II_dimension>
+  (DoFHandler<deal_II_dimension>&,
  const std::vector<unsigned int>&);
 
 //  template
 //  void
-//  DoFRenumbering::compute_component_wise<deal_II_dimension>
+//  compute_component_wise<deal_II_dimension>
 //  (std::vector<unsigned int>&,
 //   const DoFHandler<deal_II_dimension>&,
 //   const std::vector<unsigned int>&);
 
-template
-void DoFRenumbering::component_wise<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&,
- unsigned int,
- const std::vector<unsigned int>&);
+  template
+  void component_wise<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&,
  unsigned int,
  const std::vector<unsigned int>&);
 
-template
-void DoFRenumbering::component_wise<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&,
- const std::vector<unsigned int>&);
+  template
+  void component_wise<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&,
  const std::vector<unsigned int>&);
 
 // DG renumbering for DoFHandler
 
-template
-void
-DoFRenumbering::cell_wise_dg<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&,
- const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template
-void
-DoFRenumbering::compute_cell_wise_dg<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&,
- const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template
-void
-DoFRenumbering::downstream_dg<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template
-void
-DoFRenumbering::compute_downstream_dg<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-template
-void
-DoFRenumbering::compute_downstream_dg<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template void
-DoFRenumbering::cell_wise<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&,
- const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::compute_cell_wise<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&,
- const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::downstream<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
-
-template void
-DoFRenumbering::compute_downstream<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
-
-template
-void
-DoFRenumbering::clockwise_dg<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&, bool);
-
-template
-void
-DoFRenumbering::compute_clockwise_dg<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&, const DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&, const bool);
+  template
+  void
+  cell_wise_dg<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&,
  const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template
+  void
+  compute_cell_wise_dg<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&,
  const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template
+  void
+  downstream_dg<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template
+  void
+  compute_downstream_dg<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+  template
+  void
+  compute_downstream_dg<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template void
+  cell_wise<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&,
  const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  compute_cell_wise<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&,
  const std::vector<DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  downstream<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
+
+  template void
+  compute_downstream<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,std::vector<unsigned int>&,
  const DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&);
+
+  template
+  void
+  clockwise_dg<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension>&, const Point<deal_II_dimension>&, bool);
+
+  template
+  void
+  compute_clockwise_dg<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&, const DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&, const bool);
 
 // Renumbering for hp::DoFHandler
 
-template void
-DoFRenumbering::cell_wise_dg<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::compute_cell_wise_dg<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::cell_wise<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::compute_cell_wise<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
-
-template void
-DoFRenumbering::downstream_dg<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template void
-DoFRenumbering::compute_downstream_dg<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-template
-void
-DoFRenumbering::compute_downstream_dg<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template void
-DoFRenumbering::downstream<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template void
-DoFRenumbering::compute_downstream<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&);
-
-template
-void
-DoFRenumbering::clockwise_dg<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&, bool);
-
-template
-void
-DoFRenumbering::compute_clockwise_dg<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension>&,
- const Point<deal_II_dimension>&,
- const bool);
+  template void
+  cell_wise_dg<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  compute_cell_wise_dg<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  cell_wise<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  compute_cell_wise<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const std::vector<hp::DoFHandler<deal_II_dimension>::cell_iterator>&);
+
+  template void
+  downstream_dg<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template void
+  compute_downstream_dg<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+  template
+  void
+  compute_downstream_dg<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template void
+  downstream<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template void
+  compute_downstream<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&);
+
+  template
+  void
+  clockwise_dg<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&, bool);
+
+  template
+  void
+  compute_clockwise_dg<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
  const hp::DoFHandler<deal_II_dimension>&,
  const Point<deal_II_dimension>&,
  const bool);
 
 // MG
 
-template
-void DoFRenumbering::downstream_dg<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&, const unsigned int,
- const Point<deal_II_dimension>&);
+  template
+  void downstream_dg<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&, const unsigned int,
  const Point<deal_II_dimension>&);
 
-template
-void DoFRenumbering::downstream<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&, const unsigned int,
- const Point<deal_II_dimension>&);
+  template
+  void downstream<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&, const unsigned int,
  const Point<deal_II_dimension>&);
 
-template
-void DoFRenumbering::clockwise_dg<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&,
- const unsigned int,
- const Point<deal_II_dimension>&, bool);
+  template
+  void clockwise_dg<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&,
  const unsigned int,
  const Point<deal_II_dimension>&, bool);
 
 // Generic numbering schemes
 
-template
-void DoFRenumbering::random<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension> &);
-
-template
-void
-DoFRenumbering::compute_random<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension> &);
-
-template
-void DoFRenumbering::sort_selected_dofs_back<DoFHandler<deal_II_dimension> >
-(DoFHandler<deal_II_dimension> &,
- const std::vector<bool> &);
-
-template
-void
-DoFRenumbering::compute_sort_selected_dofs_back<DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension> &,
- const std::vector<bool> &);
-
-template
-void DoFRenumbering::random<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension> &);
-
-template
-void
-DoFRenumbering::compute_random<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension> &);
-
-template
-void DoFRenumbering::sort_selected_dofs_back<hp::DoFHandler<deal_II_dimension> >
-(hp::DoFHandler<deal_II_dimension> &,
- const std::vector<bool> &);
-
-template
-void
-DoFRenumbering::compute_sort_selected_dofs_back<hp::DoFHandler<deal_II_dimension> >
-(std::vector<unsigned int>&,
- const hp::DoFHandler<deal_II_dimension> &,
- const std::vector<bool> &);
-
-template
-void DoFRenumbering::subdomain_wise<deal_II_dimension>
-(DoFHandler<deal_II_dimension> &);
-
-template
-void
-DoFRenumbering::compute_subdomain_wise<deal_II_dimension>
-(std::vector<unsigned int>&,
- const DoFHandler<deal_II_dimension> &);
-
-template
-void DoFRenumbering::Cuthill_McKee<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&,
- const unsigned int,
- const bool,
- const std::vector<unsigned int>&);
-
-template
-void DoFRenumbering::cell_wise_dg<deal_II_dimension>
-(MGDoFHandler<deal_II_dimension>&,
- const unsigned int,
- const std::vector<MGDoFHandler<deal_II_dimension>::cell_iterator>&);
-template
-void DoFRenumbering::compute_cell_wise_dg<deal_II_dimension>
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const MGDoFHandler<deal_II_dimension>&, const unsigned int,
- const std::vector<MGDoFHandler<deal_II_dimension>::cell_iterator>&);
-template
-void DoFRenumbering::compute_downstream_dg<deal_II_dimension>
-(std::vector<unsigned int>&, std::vector<unsigned int>&,
- const MGDoFHandler<deal_II_dimension>&, const unsigned int,
- const Point<deal_II_dimension>&);
+  template
+  void random<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension> &);
+
+  template
+  void
+  compute_random<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
+   const DoFHandler<deal_II_dimension> &);
+
+  template
+  void sort_selected_dofs_back<DoFHandler<deal_II_dimension> >
+  (DoFHandler<deal_II_dimension> &,
+   const std::vector<bool> &);
+
+  template
+  void
+  compute_sort_selected_dofs_back<DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
+   const DoFHandler<deal_II_dimension> &,
+   const std::vector<bool> &);
+
+  template
+  void random<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension> &);
+
+  template
+  void
+  compute_random<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
+   const hp::DoFHandler<deal_II_dimension> &);
+
+  template
+  void sort_selected_dofs_back<hp::DoFHandler<deal_II_dimension> >
+  (hp::DoFHandler<deal_II_dimension> &,
+   const std::vector<bool> &);
+
+  template
+  void
+  compute_sort_selected_dofs_back<hp::DoFHandler<deal_II_dimension> >
+  (std::vector<unsigned int>&,
+   const hp::DoFHandler<deal_II_dimension> &,
+   const std::vector<bool> &);
+
+  template
+  void subdomain_wise<deal_II_dimension>
+  (DoFHandler<deal_II_dimension> &);
+
+  template
+  void
+  compute_subdomain_wise<deal_II_dimension>
+  (std::vector<unsigned int>&,
+   const DoFHandler<deal_II_dimension> &);
+
+  template
+  void Cuthill_McKee<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&,
+   const unsigned int,
+   const bool,
+   const std::vector<unsigned int>&);
+
+  template
+  void cell_wise_dg<deal_II_dimension>
+  (MGDoFHandler<deal_II_dimension>&,
+   const unsigned int,
+   const std::vector<MGDoFHandler<deal_II_dimension>::cell_iterator>&);
+  template
+  void compute_cell_wise_dg<deal_II_dimension>
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
+   const MGDoFHandler<deal_II_dimension>&, const unsigned int,
+   const std::vector<MGDoFHandler<deal_II_dimension>::cell_iterator>&);
+  template
+  void compute_downstream_dg<deal_II_dimension>
+  (std::vector<unsigned int>&, std::vector<unsigned int>&,
+   const MGDoFHandler<deal_II_dimension>&, const unsigned int,
+   const Point<deal_II_dimension>&);
+
+}  // namespace DoFRenumbering
 
 DEAL_II_NAMESPACE_CLOSE
index bed76caf08e2cae04b98df2d1e5606435528cc01..a18014e31a5996e2b994c11b27d57dedd2bb667d 100644 (file)
@@ -283,6 +283,21 @@ constraints individually.
 <h3>deal.II</h3>
 
 <ol> 
+  <li> <p>New: There is now a namespace DoFRenumbering::boost that contains the
+  implementation of interfaces to three reordering strategies provided by the
+  Boost Graph Library, namely DoFRenumbering::boost::Cuthill_McKee,
+  DoFRenumbering::boost::king_ordering, and DoFRenumbering::boost::minimum_degree.
+  <br>
+  (WB 2008/03/07)
+  </p></li>
+
+  <li> <p>Changed: The DoFRenumbering class has been converted into a namespace.
+  Since all its members were static, this hasn't changed anything fundamental
+  and use of the functions therein should work in exactly the same way as before.
+  <br>
+  (WB 2008/03/07)
+  </p></li>
+
   <li> <p>New: The new TriaObjectAccessor::set_all_boundary_indicators function
   does not only set the boundary indicator of the current object but of all that
   bound it as well. For example, in 3d, if TriaObjectAccessor::set_boundary_indicator
@@ -293,8 +308,7 @@ constraints individually.
   (WB 2008/02/26)
   </p></li>
 
-
-  <li> <p>Improved: The tria pointer and the debug_output_grid function in GridIn
+  <li> <p>Improved: The tria pointer and the GridIn::debug_output_grid function
   are changed from private to protected to be accessible in derived classes.
   <br>
   (RH 2008/02/24)

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.