From: bangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Sun, 4 Nov 2012 23:40:43 +0000 (+0000)
Subject: Merge from upstream.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c771306ccfee86716d71cbc497baaf67237477ca;p=dealii-svn.git

Merge from upstream.

git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27372 0785d39b-7218-0410-832d-ea1e28bc413d
---

c771306ccfee86716d71cbc497baaf67237477ca
diff --cc deal.II/include/deal.II/dofs/dof_accessor.templates.h
index ae4b831bc8,cefac603cc..8643bd7694
--- a/deal.II/include/deal.II/dofs/dof_accessor.templates.h
+++ b/deal.II/include/deal.II/dofs/dof_accessor.templates.h
@@@ -736,42 -734,44 +736,44 @@@ namespace interna
                                               // leaves out is faces in 3d --
                                               // let's hope that this never is
                                               // a problem
- 	        Assert ((dim==structdim
- 		             ?
- 		             typename
- 		             dealii::DoFHandler<dim,spacedim>::
- 		             raw_cell_iterator (&dof_handler.get_tria(),
- 					                    obj_level,
- 					                    obj_index,
- 					                    &dof_handler)->used()
- 		             :
- 		             (structdim==1
- 		              ?
- 		              typename
- 		              dealii::DoFHandler<dim,spacedim>::
- 		              raw_line_iterator (&dof_handler.get_tria(),
- 					                     obj_level,
- 					                     obj_index,
- 					                     &dof_handler)->used()
- 		              :
- 		              true))
- 		            == true,
- 		            ExcMessage ("This cell is not active and therefore can't be "
- 				                "queried for its active FE indices"));
 -            Assert ((dim==structdim
 -                     ?
 -                     typename
 -		     dealii::internal::DoFHandler::
 -		     Iterators<dealii::DoFHandler<dim,spacedim> >::
 -                     raw_cell_iterator (&dof_handler.get_tria(),
 -                                        obj_level,
 -                                        obj_index,
 -                                        &dof_handler)->used()
 -                     :
 -                     (structdim==1
 -                      ?
 -                      typename
 -                      dealii::internal::DoFHandler::
 -		      Iterators<dealii::DoFHandler<dim,spacedim> >::
 -                      raw_line_iterator (&dof_handler.get_tria(),
 -                                         obj_level,
 -                                         obj_index,
 -                                         &dof_handler)->used()
 -                      :
 -                      true))
 -                    == true,
 -                    ExcMessage ("This cell is not active and therefore can't be "
 -                                "queried for its active FE indices"));
 -            return 1;
 -          }
 -
 -
 -
 -        template <int structdim, int dim, int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::DoFHandler<dim,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<structdim> &)
 -          {
++		Assert ((dim==structdim
++			 ?
++			 typename
++			 dealii::internal::DoFHandler::
++			 Iterators<dealii::DoFHandler<dim,spacedim> >::
++			 raw_cell_iterator (&dof_handler.get_tria(),
++					    obj_level,
++					    obj_index,
++					    &dof_handler)->used()
++			 :
++			 (structdim==1
++			  ?
++			  typename
++			  dealii::internal::DoFHandler::
++			  Iterators<dealii::DoFHandler<dim,spacedim> >::
++			  raw_line_iterator (&dof_handler.get_tria(),
++					     obj_level,
++					     obj_index,
++					     &dof_handler)->used()
++			  :
++			  true))
++			== true,
++			ExcMessage ("This cell is not active and therefore can't be "
++				    "queried for its active FE indices"));
 +	        return 1;
 +	      }
 +
 +
 +
 +	    template <int structdim, int dim, int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::DoFHandler<dim,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<structdim> &)
 +	      {
                                               // check that the object we look
                                               // at is in fact active. the
                                               // problem is that we have
@@@ -786,310 -786,312 +788,312 @@@
                                               // leaves out is faces in 3d --
                                               // let's hope that this never is
                                               // a problem
 -            Assert ((dim==structdim
 -                     ?
 -                     typename
 -		     dealii::internal::DoFHandler::
 -		     Iterators<dealii::DoFHandler<dim,spacedim> >::
 -                     raw_cell_iterator (&dof_handler.get_tria(),
 -                                        obj_level,
 -                                        obj_index,
 -                                        &dof_handler)->used()
 -                     :
 -                     (structdim==1
 -                      ?
 -                      typename
 -		      dealii::internal::DoFHandler::
 -		      Iterators<dealii::DoFHandler<dim,spacedim> >::
 -                      raw_line_iterator (&dof_handler.get_tria(),
 -                                         obj_level,
 -                                         obj_index,
 -                                         &dof_handler)->used()
 -                      :
 -                      true))
 -                    == true,
 -                    ExcMessage ("This cell is not active and therefore can't be "
 -                                "queried for its active FE indices"));
 -            Assert (n == 0, ExcIndexRange (n, 0, 1));
 -
 -            return dealii::DoFHandler<dim,spacedim>::default_fe_index;
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 -                                                                           obj_index,
 -                                                                           fe_index,
 -                                                                           obj_level);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 -                                                                             obj_index);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 -                                                                             obj_level,
 -                                                                             obj_index,
 -                                                                             n);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.fe_index_is_active(dof_handler,
 -                                                               obj_index,
 -                                                               fe_index,
 -                                                               obj_level);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                             const unsigned int ,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
 -                                                                 obj_index);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
 -                                                                 obj_level,
 -                                                                 obj_index,
 -                                                                 n);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 -                                                                           obj_index,
 -                                                                           fe_index,
 -                                                                           obj_level);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 -                                                                             obj_index);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 -                                                                             obj_level,
 -                                                                             obj_index,
 -                                                                             n);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.fe_index_is_active(dof_handler,
 -                                                               obj_index,
 -                                                               fe_index,
 -                                                               obj_level);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
 -                                                                 obj_index);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<1> &)
 -          {
 -            return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
 -                                                                 obj_level,
 -                                                                 obj_index,
 -                                                                 n);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.faces->quads.fe_index_is_active(dof_handler,
 -                                                               obj_index,
 -                                                               fe_index,
 -                                                               obj_level);
 -          }
 -
 -        template <int spacedim>
 -        static
 -        bool
 -        fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                            const unsigned int obj_level,
 -                            const unsigned int obj_index,
 -                            const unsigned int fe_index,
 -                            const dealii::internal::int2type<3> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 -                                                                           obj_index,
 -                                                                           fe_index,
 -                                                                           obj_level);
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int ,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.faces->quads.n_active_fe_indices (dof_handler,
 -                                                                 obj_index);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<2> &)
 -          {
 -            return dof_handler.faces->quads.nth_active_fe_index (dof_handler,
 -                                                                 obj_level,
 -                                                                 obj_index,
 -                                                                 n);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const dealii::internal::int2type<3> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 -                                                                             obj_index);
 -          }
 -
 -
 -
 -        template <int spacedim>
 -        static
 -        unsigned int
 -        nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 -                             const unsigned int obj_level,
 -                             const unsigned int obj_index,
 -                             const unsigned int n,
 -                             const dealii::internal::int2type<3> &)
 -          {
 -            return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 -                                                                             obj_level,
 -                                                                             obj_index,
 -                                                                             n);
 -          }
 +	        Assert ((dim==structdim
- 		             ?
- 		             typename
- 		             dealii::DoFHandler<dim,spacedim>::
- 		             raw_cell_iterator (&dof_handler.get_tria(),
- 					                    obj_level,
- 					                    obj_index,
- 					                    &dof_handler)->used()
- 		             :
- 		             (structdim==1
- 		              ?
- 		              typename
- 		              dealii::DoFHandler<dim,spacedim>::
- 		              raw_line_iterator (&dof_handler.get_tria(),
- 					                     obj_level,
- 					                     obj_index,
- 				   	                     &dof_handler)->used()
- 		              :
- 		              true))
- 		            == true,
- 		            ExcMessage ("This cell is not active and therefore can't be "
- 			  	                "queried for its active FE indices"));
++			 ?
++			 typename
++			 dealii::internal::DoFHandler::
++			 Iterators<dealii::DoFHandler<dim,spacedim> >::
++			 raw_cell_iterator (&dof_handler.get_tria(),
++					    obj_level,
++					    obj_index,
++					    &dof_handler)->used()
++			 :
++			 (structdim==1
++			  ?
++			  typename
++			  dealii::internal::DoFHandler::
++			  Iterators<dealii::DoFHandler<dim,spacedim> >::
++			  raw_line_iterator (&dof_handler.get_tria(),
++					     obj_level,
++					     obj_index,
++					     &dof_handler)->used()
++			  :
++			  true))
++			== true,
++			ExcMessage ("This cell is not active and therefore can't be "
++				    "queried for its active FE indices"));
 +	        Assert (n == 0, ExcIndexRange (n, 0, 1));
 +
 +	        return dealii::DoFHandler<dim,spacedim>::default_fe_index;
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +      	bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 +					 				                                       obj_index,
 +									                                       fe_index,
 +									                                       obj_level);
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 +									                                         obj_index);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<1,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 +									                                         obj_level,
 +									                                         obj_index,
 +									                                         n);
 +	      }
 +
 +
 +	    template <int spacedim>
 +  	    static
 +     	bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.fe_index_is_active(dof_handler,
 +							                                   obj_index,
 +							                                   fe_index,
 +							                                   obj_level);
 +	      }
 +
 +
 +  	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                 const unsigned int ,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
 +								                                 obj_index);
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
 +								                                 obj_level,
 +								                                 obj_index,
 +								                                 n);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 +									                                       obj_index,
 +									                                       fe_index,
 +									                                       obj_level);
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 +									                                         obj_index);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<2,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 +									                                         obj_level,
 +									                                         obj_index,
 +									                                         n);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.fe_index_is_active(dof_handler,
 +							                                   obj_index,
 +							                                   fe_index,
 +							                                   obj_level);
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.n_active_fe_indices (dof_handler,
 +								                                 obj_index);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<1> &)
 +	      {
 +	        return dof_handler.faces->lines.nth_active_fe_index (dof_handler,
 +								                                 obj_level,
 +							 	                                 obj_index,
 +								                                 n);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.faces->quads.fe_index_is_active(dof_handler,
 +							                                   obj_index,
 +							                                   fe_index,
 +							                                   obj_level);
 +	      }
 +
 +	    template <int spacedim>
 +	    static
 +	    bool
 +	    fe_index_is_active (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                const unsigned int obj_level,
 +			                const unsigned int obj_index,
 +			                const unsigned int fe_index,
 +			                const dealii::internal::int2type<3> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.fe_index_is_active(dof_handler,
 +									                                       obj_index,
 +									                                       fe_index,
 +									                                       obj_level);
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int ,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.faces->quads.n_active_fe_indices (dof_handler,
 +				  			  	                                 obj_index);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<2> &)
 +	      {
 +	        return dof_handler.faces->quads.nth_active_fe_index (dof_handler,
 +								                                 obj_level,
 +							 	                                 obj_index,
 +								                                 n);
 +	      }
 +
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    n_active_fe_indices (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const dealii::internal::int2type<3> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.n_active_fe_indices (dof_handler,
 +									                                         obj_index);
 +	      }
 +
 +
 +
 +  	    template <int spacedim>
 +	    static
 +	    unsigned int
 +	    nth_active_fe_index (const dealii::hp::DoFHandler<3,spacedim> &dof_handler,
 +			                 const unsigned int obj_level,
 +			                 const unsigned int obj_index,
 +			                 const unsigned int n,
 +			                 const dealii::internal::int2type<3> &)
 +	      {
 +	        return dof_handler.levels[obj_level]->dof_object.nth_active_fe_index (dof_handler,
 +									                                         obj_level,
 +									                                         obj_index,
 +									                                         n);
 +	      }
  
                                           /**
                                            * Set the @p local_index-th
@@@ -1736,72 -1699,15 +1740,72 @@@ namespace interna
                                         // applies, if the face_rotation or
                                         // face_orientation is non-standard
        for (unsigned int quad=0; quad<6; ++quad)
 -        for (unsigned int d=0; d<dofs_per_quad; ++d)
 -          *next++ = accessor.quad(quad)->dof_index(accessor.get_fe(fe_index).
 -                                                   adjust_quad_dof_index_for_face_orientation(d,
 -                                                                                              accessor.face_orientation(quad),
 -                                                                                              accessor.face_flip(quad),
 -                                                                                              accessor.face_rotation(quad)),
 -                                                   fe_index);
 +	    for (unsigned int d=0; d<dofs_per_quad; ++d)
 +	      *next++ = accessor.quad(quad)->dof_index(accessor.get_fe(fe_index).
 +						                           adjust_quad_dof_index_for_face_orientation(d,
 +											                                                  accessor.face_orientation(quad),
 +											                                                  accessor.face_flip(quad),
 +											                                                  accessor.face_rotation(quad)),
 +						                           fe_index);
        for (unsigned int d=0; d<dofs_per_hex; ++d)
 -        *next++ = accessor.dof_index(d,fe_index);
 +	    *next++ = accessor.dof_index(d,fe_index);
 +    }
-     
++
 +    template<class DH>
 +    void get_mg_dof_indices (const dealii::DoFAccessor<1, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
 +      const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
 +      std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-       
++
 +      for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +          *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-       
++
 +      for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +        *next++ = accessor.mg_dof_index (level, dof);
-       
++
 +      Assert (next == dof_indices.end (), ExcInternalError ());
 +    }
-     
++
 +    template<class DH>
 +    void get_mg_dof_indices (const dealii::DoFAccessor<2, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
 +      const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
 +      std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-       
++
 +      for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +          *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-       
++
 +      for (unsigned int line = 0; line < GeometryInfo<2>::lines_per_cell; ++line)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +          *next++ = accessor.line (line)->mg_dof_index (level, dof);
-       
++
 +      for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +        *next++ = accessor.mg_dof_index (level, dof);
-       
++
 +      Assert (next == dof_indices.end (), ExcInternalError ());
 +    }
-     
++
 +    template<class DH>
 +    void get_mg_dof_indices (const dealii::DoFAccessor<3, DH>& accessor, const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) {
 +      const FiniteElement<DH::dimension, DH::space_dimension>& fe = accessor.get_dof_handler ().get_fe ()[fe_index];
 +      std::vector<types::global_dof_index>::iterator next = dof_indices.begin ();
-       
++
 +      for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +          *next++ = accessor.mg_vertex_dof_index (level, vertex, dof);
-       
++
 +      for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +          *next++ = accessor.line (line)->mg_dof_index (level, dof);
-       
++
 +      for (unsigned int quad = 0; quad < GeometryInfo<3>::quads_per_cell; ++quad)
 +        for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +          *next++ = accessor.quad (quad)->mg_dof_index (level, dof);
-       
++
 +      for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof)
 +        *next++ = accessor.mg_dof_index (level, dof);
-       
++
 +      Assert (next == dof_indices.end (), ExcInternalError ());
      }
    }
  }
@@@ -1871,34 -1777,9 +1875,34 @@@ DoFAccessor<structdim,DH>::get_dof_indi
    dealii::internal::DoFAccessor::get_dof_indices (*this, dof_indices, fe_index);
  }
  
 +template<int structdim, class DH>
 +inline
 +void DoFAccessor<structdim, DH>::get_mg_dof_indices (const int level, std::vector<types::global_dof_index>& dof_indices, const unsigned int fe_index) const {
 +  Assert (this->dof_handler != 0, ExcInvalidObject ());
 +  Assert (&this->dof_handler->get_fe () != 0, ExcInvalidObject ());
-   
++
 +  switch (structdim) {
 +    case 1: {
 +	  Assert (dof_indices.size () == 2 * this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + this->dof_handler->get_fe ()[fe_index].dofs_per_line, ExcVectorDoesNotMatch ());
 +	  break;
 +    }
-     
++
 +    case 2: {
 +	  Assert (dof_indices.size () == 4 * (this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + this->dof_handler->get_fe ()[fe_index].dofs_per_line) + this->dof_handler->get_fe ()[fe_index].dofs_per_quad, ExcVectorDoesNotMatch ());
 +	  break;
 +    }
-     
+ 
 +    case 3: {
 +	  Assert (dof_indices.size () == 8 * this->dof_handler->get_fe ()[fe_index].dofs_per_vertex + 12 * this->dof_handler->get_fe ()[fe_index].dofs_per_line + 6 * this->dof_handler->get_fe ()[fe_index].dofs_per_quad + this->dof_handler->get_fe ()[fe_index].dofs_per_hex, ExcVectorDoesNotMatch ());
 +	  break;
 +    }
-     
+ 
 +    default:
 +	  Assert (false, ExcNotImplemented ());
 +  }
-   
+ 
 +  internal::DoFAccessor::get_mg_dof_indices (*this, level, dof_indices, fe_index);
 +}
  
  template <int structdim, class DH>
  inline
diff --cc deal.II/include/deal.II/lac/constraint_matrix.h
index 8fabf6a013,a79738f20e..af4f3ad83d
--- a/deal.II/include/deal.II/lac/constraint_matrix.h
+++ b/deal.II/include/deal.II/lac/constraint_matrix.h
@@@ -381,10 -381,11 +381,11 @@@ class ConstraintMatrix : public Subscri
                                        * constraint line <i>i</i>, according
                                        * to the discussion in the general
                                        * class description.
-                                       * Note: the line needs to be added with
+ 				      *
+                                       * @note the line needs to be added with
                                        * one of the add_line() calls first.
                                        */
 -    void set_inhomogeneity (const unsigned int line,
 +    void set_inhomogeneity (const types::global_dof_index line,
                              const double       value);
  
                                       /**
diff --cc deal.II/include/deal.II/numerics/vector_tools.templates.h
index c7c30bc364,ae40ad9d46..10882b2e97
--- a/deal.II/include/deal.II/numerics/vector_tools.templates.h
+++ b/deal.II/include/deal.II/numerics/vector_tools.templates.h
@@@ -402,23 -401,24 +401,24 @@@ namespace VectorTool
                                         // that is actually wholly on
                                         // the boundary, not only by
                                         // one line or one vertex
-       typename dealii::DoFHandler<dim,spacedim>::active_face_iterator
-         face = dof_handler.begin_active_face(),
-         endf = dof_handler.end_face();
+       typename dealii::DoFHandler<dim,spacedim>::active_cell_iterator
+         cell = dof_handler.begin_active(),
+         endc = dof_handler.end();
 -      std::vector<unsigned int> face_dof_indices (fe.dofs_per_face);
 +      std::vector<types::global_dof_index> face_dof_indices (fe.dofs_per_face);
-       for (; face!=endf; ++face)
-         if (face->at_boundary())
-           {
-             face->get_dof_indices (face_dof_indices);
-             for (unsigned int i=0; i<fe.dofs_per_face; ++i)
-                                                // enter zero boundary values
-                                                // for all boundary nodes
-                                                //
-                                                // we need not care about
-                                                // vector valued elements here,
-                                                // since we set all components
-               boundary_values[face_dof_indices[i]] = 0.;
-           }
+       for (; cell!=endc; ++cell)
+ 	for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+ 	  if (cell->at_boundary(f))
+ 	    {
+ 	      cell->face(f)->get_dof_indices (face_dof_indices);
+ 	      for (unsigned int i=0; i<fe.dofs_per_face; ++i)
+ 						 // enter zero boundary values
+ 						 // for all boundary nodes
+ 						 //
+ 						 // we need not care about
+ 						 // vector valued elements here,
+ 						 // since we set all components
+ 		boundary_values[face_dof_indices[i]] = 0.;
+ 	    }
      }
    }
  
@@@ -1939,10 -1929,9 +1929,9 @@@
    interpolate_boundary_values (const DH                 &dof,
                                 const types::boundary_id            boundary_component,
                                 const Function<DH::space_dimension>           &boundary_function,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const ComponentMask       &component_mask)
    {
-     Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
      interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                                  dof, boundary_component,
                                  boundary_function, boundary_values, component_mask);
@@@ -1954,10 -1943,9 +1943,9 @@@
    void
    interpolate_boundary_values (const DH                 &dof,
                                 const typename FunctionMap<DH::space_dimension>::type &function_map,
 -                               std::map<unsigned int,double> &boundary_values,
 +                               std::map<types::global_dof_index,double> &boundary_values,
                                 const ComponentMask       &component_mask)
    {
-     Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
      interpolate_boundary_values(StaticMappingQ1<DH::dimension,DH::space_dimension>::mapping,
                                  dof, function_map,
                                  boundary_values, component_mask);
@@@ -2263,10 -2249,9 +2249,9 @@@
    project_boundary_values (const DoFHandler<dim,spacedim>    &dof,
                             const typename FunctionMap<spacedim>::type &boundary_functions,
                             const Quadrature<dim-1>  &q,
 -                           std::map<unsigned int,double> &boundary_values,
 +                           std::map<types::global_dof_index,double> &boundary_values,
                             std::vector<unsigned int> component_mapping)
    {
-     Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
      project_boundary_values(StaticMappingQ1<dim,spacedim>::mapping, dof, boundary_functions, q,
                              boundary_values, component_mapping);
    }
diff --cc deal.II/source/dofs/dof_handler.cc
index 0a68adaa5c,afa53fe64e..10ea886f82
--- a/deal.II/source/dofs/dof_handler.cc
+++ b/deal.II/source/dofs/dof_handler.cc
@@@ -257,427 -257,99 +257,427 @@@ namespace interna
                                            * was selected when calling
                                            * @p distribute_dofs the last time.
                                            */
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<1,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<1,spacedim>::invalid_dof_index);
 -
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels
 -                  .push_back (new internal::DoFHandler::DoFLevel<1>);
 -
 -                dof_handler.levels.back()->dof_object.dofs
 -                  .resize (dof_handler.tria->n_raw_cells(i) *
 -                           dof_handler.selected_fe->dofs_per_line,
 -                           DoFHandler<1,spacedim>::invalid_dof_index);
 -
 -                dof_handler.levels.back()->cell_dof_indices_cache
 -                  .resize (dof_handler.tria->n_raw_cells(i) *
 -                           dof_handler.selected_fe->dofs_per_cell,
 -                           DoFHandler<1,spacedim>::invalid_dof_index);
 -              }
 -          }
 -
 -
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<2,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<2,spacedim>::invalid_dof_index);
 -
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 +	    template <int spacedim>
 +	    static
 +	    void reserve_space (DoFHandler<1,spacedim> &dof_handler)
 +	      {
 +	        dof_handler.vertex_dofs
 +	          .resize(dof_handler.tria->n_vertices() *
 +		              dof_handler.selected_fe->dofs_per_vertex,
 +		              DoFHandler<1,spacedim>::invalid_dof_index);
 +
 +	        for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +	          {
- 		        dof_handler.levels
- 		          .push_back (new internal::DoFHandler::DoFLevel<1>);
- 
- 		        dof_handler.levels.back()->dof_object.dofs
- 		          .resize (dof_handler.tria->n_raw_lines(i) *
- 			               dof_handler.selected_fe->dofs_per_line,
- 			               DoFHandler<1,spacedim>::invalid_dof_index);
- 
- 		        dof_handler.levels.back()->cell_dof_indices_cache
- 		          .resize (dof_handler.tria->n_raw_lines(i) *
- 			               dof_handler.selected_fe->dofs_per_cell,
- 			               DoFHandler<1,spacedim>::invalid_dof_index);
- 	          }
++		    dof_handler.levels
++		      .push_back (new internal::DoFHandler::DoFLevel<1>);
++
++		    dof_handler.levels.back()->dof_object.dofs
++		      .resize (dof_handler.tria->n_raw_cells(i) *
++			       dof_handler.selected_fe->dofs_per_line,
++			       DoFHandler<1,spacedim>::invalid_dof_index);
++
++		    dof_handler.levels.back()->cell_dof_indices_cache
++		      .resize (dof_handler.tria->n_raw_cells(i) *
++			       dof_handler.selected_fe->dofs_per_cell,
++			       DoFHandler<1,spacedim>::invalid_dof_index);
++		  }
 +	      }
 +
 +
 +	    template <int spacedim>
 +	    static
 +	    void reserve_space (DoFHandler<2,spacedim> &dof_handler)
 +	      {
 +	        dof_handler.vertex_dofs
 +	          .resize(dof_handler.tria->n_vertices() *
 +		              dof_handler.selected_fe->dofs_per_vertex,
 +		              DoFHandler<2,spacedim>::invalid_dof_index);
 +
 +	        for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +	          {
 +		        dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<2>);
  
- 		        dof_handler.levels.back()->dof_object.dofs
- 		          .resize (dof_handler.tria->n_raw_quads(i) *
- 			               dof_handler.selected_fe->dofs_per_quad,
- 			               DoFHandler<2,spacedim>::invalid_dof_index);
+                 dof_handler.levels.back()->dof_object.dofs
+                   .resize (dof_handler.tria->n_raw_cells(i) *
+                            dof_handler.selected_fe->dofs_per_quad,
+                            DoFHandler<2,spacedim>::invalid_dof_index);
  
- 		        dof_handler.levels.back()->cell_dof_indices_cache
- 		          .resize (dof_handler.tria->n_raw_quads(i) *
- 			               dof_handler.selected_fe->dofs_per_cell,
- 			               DoFHandler<2,spacedim>::invalid_dof_index);
- 	          }
+                 dof_handler.levels.back()->cell_dof_indices_cache
+                   .resize (dof_handler.tria->n_raw_cells(i) *
+                            dof_handler.selected_fe->dofs_per_cell,
+                            DoFHandler<2,spacedim>::invalid_dof_index);
+               }
  
 -            dof_handler.faces = new internal::DoFHandler::DoFFaces<2>;
 -            dof_handler.faces->lines.dofs
 -              .resize (dof_handler.tria->n_raw_lines() *
 -                       dof_handler.selected_fe->dofs_per_line,
 -                       DoFHandler<2,spacedim>::invalid_dof_index);
 -          }
 +	        dof_handler.faces = new internal::DoFHandler::DoFFaces<2>;
 +	        dof_handler.faces->lines.dofs
 +	          .resize (dof_handler.tria->n_raw_lines() *
 +		               dof_handler.selected_fe->dofs_per_line,
 +		               DoFHandler<2,spacedim>::invalid_dof_index);
 +	      }
  
  
 -        template <int spacedim>
 -        static
 -        void reserve_space (DoFHandler<3,spacedim> &dof_handler)
 -          {
 -            dof_handler.vertex_dofs
 -              .resize(dof_handler.tria->n_vertices() *
 -                      dof_handler.selected_fe->dofs_per_vertex,
 -                      DoFHandler<3,spacedim>::invalid_dof_index);
 +	    template <int spacedim>
 +	    static
 +	    void reserve_space (DoFHandler<3,spacedim> &dof_handler)
 +	      {
 +	        dof_handler.vertex_dofs
 +	          .resize(dof_handler.tria->n_vertices() *
 +		              dof_handler.selected_fe->dofs_per_vertex,
 +		              DoFHandler<3,spacedim>::invalid_dof_index);
  
 -            for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 -              {
 -                dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>);
 +	        for (unsigned int i=0; i<dof_handler.tria->n_levels(); ++i)
 +	          {
 +		        dof_handler.levels.push_back (new internal::DoFHandler::DoFLevel<3>);
  
- 		        dof_handler.levels.back()->dof_object.dofs
- 		          .resize (dof_handler.tria->n_raw_hexs(i) *
- 			               dof_handler.selected_fe->dofs_per_hex,
- 			               DoFHandler<3,spacedim>::invalid_dof_index);
+                 dof_handler.levels.back()->dof_object.dofs
+                   .resize (dof_handler.tria->n_raw_cells(i) *
+                            dof_handler.selected_fe->dofs_per_hex,
+                            DoFHandler<3,spacedim>::invalid_dof_index);
  
- 		        dof_handler.levels.back()->cell_dof_indices_cache
- 		          .resize (dof_handler.tria->n_raw_hexs(i) *
- 			               dof_handler.selected_fe->dofs_per_cell,
- 			               DoFHandler<3,spacedim>::invalid_dof_index);
- 	          }
- 	        dof_handler.faces = new internal::DoFHandler::DoFFaces<3>;
+                 dof_handler.levels.back()->cell_dof_indices_cache
+                   .resize (dof_handler.tria->n_raw_cells(i) *
+                            dof_handler.selected_fe->dofs_per_cell,
+                            DoFHandler<3,spacedim>::invalid_dof_index);
+               }
+             dof_handler.faces = new internal::DoFHandler::DoFFaces<3>;
  
 -            dof_handler.faces->lines.dofs
 -              .resize (dof_handler.tria->n_raw_lines() *
 -                       dof_handler.selected_fe->dofs_per_line,
 -                       DoFHandler<3,spacedim>::invalid_dof_index);
 -            dof_handler.faces->quads.dofs
 -              .resize (dof_handler.tria->n_raw_quads() *
 -                       dof_handler.selected_fe->dofs_per_quad,
 -                       DoFHandler<3,spacedim>::invalid_dof_index);
 -          }
 +	        dof_handler.faces->lines.dofs
 +	          .resize (dof_handler.tria->n_raw_lines() *
 +		               dof_handler.selected_fe->dofs_per_line,
 +		               DoFHandler<3,spacedim>::invalid_dof_index);
 +	        dof_handler.faces->quads.dofs
 +	          .resize (dof_handler.tria->n_raw_quads() *
 +		               dof_handler.selected_fe->dofs_per_quad,
 +		               DoFHandler<3,spacedim>::invalid_dof_index);
 +	      }
 +
 +	    template<int spacedim>
 +	    static
 +	    void reserve_space_mg (DoFHandler<1, spacedim>& dof_handler) {
 +	      Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +	      dof_handler.clear_mg_space ();
 +
 +	      const dealii::Triangulation<1, spacedim>& tria = dof_handler.get_tria ();
 +	      const unsigned int& dofs_per_line = dof_handler.get_fe ().dofs_per_line;
 +          const unsigned int& n_levels = tria.n_levels ();
 +
 +	      for (unsigned int i = 0; i < n_levels; ++i) {
 +	        dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<1>);
 +	        dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<unsigned int> (tria.n_raw_lines (i) * dofs_per_line, DoFHandler<1>::invalid_dof_index);
 +	      }
 +
 +	      const unsigned int& n_vertices = tria.n_vertices ();
 +
 +	      dof_handler.mg_vertex_dofs.resize (n_vertices);
 +
 +	      std::vector<unsigned int> max_level (n_vertices, 0);
 +	      std::vector<unsigned int> min_level (n_vertices, n_levels);
 +
 +	      for (typename dealii::Triangulation<1, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +	        const unsigned int level = cell->level ();
 +
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) {
 +	          const unsigned int vertex_index = cell->vertex_index (vertex);
 +
 +	          if (min_level[vertex_index] > level)
 +	            min_level[vertex_index] = level;
 +
 +	          if (max_level[vertex_index] < level)
 +	            max_level[vertex_index] = level;
 +	        }
 +	      }
 +
 +	      for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +	        if (tria.vertex_used (vertex)) {
 +	          Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], dof_handler.get_fe ().dofs_per_vertex);
 +	        }
 +
 +	        else {
 +	          Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] == 0, ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +	        }
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void reserve_space_mg (DoFHandler<2, spacedim>& dof_handler) {
 +	      Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +	      dof_handler.clear_mg_space ();
 +
 +	      const dealii::FiniteElement<2, spacedim>& fe = dof_handler.get_fe ();
 +	      const dealii::Triangulation<2, spacedim>& tria = dof_handler.get_tria ();
 +	      const unsigned int& n_levels = tria.n_levels ();
 +
 +	      for (unsigned int i = 0; i < n_levels; ++i) {
 +	        dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<2>);
 +	        dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<unsigned int> (tria.n_raw_quads (i) * fe.dofs_per_quad, DoFHandler<2>::invalid_dof_index);
 +	      }
 +
 +	      dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<2>;
 +	      dof_handler.mg_faces->lines.dofs = std::vector<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<2>::invalid_dof_index);
 +
 +	      const unsigned int& n_vertices = tria.n_vertices ();
 +
 +	      dof_handler.mg_vertex_dofs.resize (n_vertices);
 +
 +	      std::vector<unsigned int> max_level (n_vertices, 0);
 +	      std::vector<unsigned int> min_level (n_vertices, n_levels);
 +
 +	      for (typename dealii::Triangulation<2, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +	        const unsigned int level = cell->level ();
 +
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex) {
 +	          const unsigned int vertex_index = cell->vertex_index (vertex);
 +
 +	          if (min_level[vertex_index] > level)
 +	            min_level[vertex_index] = level;
 +
 +   	          if (max_level[vertex_index] < level)
 +	            max_level[vertex_index] = level;
 +	        }
 +	      }
 +
 +	      for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +	        if (tria.vertex_used (vertex)) {
 +	          Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
 +	        }
 +
 +	        else {
 +	          Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] == 0, ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +	        }
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void reserve_space_mg (DoFHandler<3, spacedim>& dof_handler) {
 +	      Assert (dof_handler.get_tria ().n_levels () > 0, ExcMessage ("Invalid triangulation"));
 +	      dof_handler.clear_mg_space ();
 +
 +	      const dealii::FiniteElement<3, spacedim>& fe = dof_handler.get_fe ();
 +	      const dealii::Triangulation<3, spacedim>& tria = dof_handler.get_tria ();
 +	      const unsigned int& n_levels = tria.n_levels ();
 +
 +	      for (unsigned int i = 0; i < n_levels; ++i) {
 +	        dof_handler.mg_levels.push_back (new internal::DoFHandler::DoFLevel<3>);
 +	        dof_handler.mg_levels.back ()->dof_object.dofs = std::vector<unsigned int> (tria.n_raw_hexs (i) * fe.dofs_per_hex, DoFHandler<3>::invalid_dof_index);
 +	      }
 +
 +	      dof_handler.mg_faces = new internal::DoFHandler::DoFFaces<3>;
 +	      dof_handler.mg_faces->lines.dofs = std::vector<unsigned int> (tria.n_raw_lines () * fe.dofs_per_line, DoFHandler<3>::invalid_dof_index);
 +	      dof_handler.mg_faces->quads.dofs = std::vector<unsigned int> (tria.n_raw_quads () * fe.dofs_per_quad, DoFHandler<3>::invalid_dof_index);
 +
 +	      const unsigned int& n_vertices = tria.n_vertices ();
 +
 +	      dof_handler.mg_vertex_dofs.resize (n_vertices);
 +
 +	      std::vector<unsigned int> max_level (n_vertices, 0);
 +	      std::vector<unsigned int> min_level (n_vertices, n_levels);
 +
 +	      for (typename dealii::Triangulation<3, spacedim>::cell_iterator cell = tria.begin (); cell != tria.end (); ++cell) {
 +	        const unsigned int level = cell->level ();
 +
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex) {
 +	          const unsigned int vertex_index = cell->vertex_index (vertex);
 +
 +	          if (min_level[vertex_index] > level)
 +	            min_level[vertex_index] = level;
 +
 +	          if (max_level[vertex_index] < level)
 +	            max_level[vertex_index] = level;
 +	        }
 +	      }
 +
 +	      for (unsigned int vertex = 0; vertex < n_vertices; ++vertex)
 +	        if (tria.vertex_used (vertex)) {
 +	          Assert (min_level[vertex] < n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] >= min_level[vertex], ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (min_level[vertex], max_level[vertex], fe.dofs_per_vertex);
 +	        }
 +
 +	        else {
 +	          Assert (min_level[vertex] == n_levels, ExcInternalError ());
 +	          Assert (max_level[vertex] == 0, ExcInternalError ());
 +	          dof_handler.mg_vertex_dofs[vertex].init (1, 0, 0);
 +	        }
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<1, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +	  	  const FiniteElement<1, spacedim>& fe = cell->get_fe ();
 +
 +	      if (fe.dofs_per_vertex > 0)
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<1>::vertices_per_cell; ++vertex) {
 +	          typename DoFHandler<1, spacedim>::cell_iterator neighbor = cell->neighbor (vertex);
 +
 +	          if (neighbor.state () == IteratorState::valid)
 +	            if (neighbor->user_flag_set () && (neighbor->level () == cell->level ())) {
 +	              if (vertex == 0)
 +	                for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +	                  cell->set_mg_vertex_dof_index (cell->level (), 0, dof, neighbor->mg_vertex_dof_index (cell->level (), 1, dof));
 +
 +	              else
 +	                for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +	                  cell->set_mg_vertex_dof_index (cell->level (), 1, dof, neighbor->mg_vertex_dof_index (cell->level (), 0, dof));
 +
 +	              continue;
 +	            }
 +
 +	          for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +	            cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
 +	        }
 +
 +	      if (fe.dofs_per_line > 0)
 +	        for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +	          cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +
 +	      cell->set_user_flag ();
 +	      return next_free_dof;
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<2, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +	      const FiniteElement<2, spacedim>& fe = cell->get_fe ();
 +
 +	      if (fe.dofs_per_vertex > 0)
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<2>::vertices_per_cell; ++vertex)
 +	          if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<2>::invalid_dof_index)
 +	            for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +	              cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
 +
 +	      if (fe.dofs_per_line > 0)
 +	        for (unsigned int face = 0; face < GeometryInfo<2>::faces_per_cell; ++face) {
 +	          typename DoFHandler<2, spacedim>::line_iterator line = cell->line (face);
 +
 +	          if (line->mg_dof_index (cell->level (), 0) == DoFHandler<2>::invalid_dof_index)
 +	            for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +	              line->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +	        }
 +
 +	      if (fe.dofs_per_quad > 0)
 +	        for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +	          cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +
 +	      cell->set_user_flag ();
 +	      return next_free_dof;
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index distribute_dofs_on_cell (typename DoFHandler<3, spacedim>::cell_iterator& cell, types::global_dof_index next_free_dof) {
 +	      const FiniteElement<3, spacedim>& fe = cell->get_fe ();
 +
 +	      if (fe.dofs_per_vertex > 0)
 +	        for (unsigned int vertex = 0; vertex < GeometryInfo<3>::vertices_per_cell; ++vertex)
 +	          if (cell->mg_vertex_dof_index (cell->level (), vertex, 0) == DoFHandler<3>::invalid_dof_index)
 +	            for (unsigned int dof = 0; dof < fe.dofs_per_vertex; ++dof)
 +	              cell->set_mg_vertex_dof_index (cell->level (), vertex, dof, next_free_dof++);
 +
 +	      if (fe.dofs_per_line > 0)
 +	        for (unsigned int line = 0; line < GeometryInfo<3>::lines_per_cell; ++line) {
 +	          typename DoFHandler<3, spacedim>::line_iterator line_it = cell->line (line);
 +
 +	          if (line_it->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index)
 +	            for (unsigned int dof = 0; dof < fe.dofs_per_line; ++dof)
 +	              line_it->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +	        }
 +
 +	      if (fe.dofs_per_quad > 0)
 +	        for (unsigned int face = 0; face < GeometryInfo<3>::quads_per_cell; ++face) {
 +	          typename DoFHandler<3, spacedim>::quad_iterator quad = cell->quad (face);
 +
 +	          if (quad->mg_dof_index (cell->level (), 0) == DoFHandler<3>::invalid_dof_index)
 +	            for (unsigned int dof = 0; dof < fe.dofs_per_quad; ++dof)
 +	              quad->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +	        }
 +
 +	      if (fe.dofs_per_hex > 0)
 +	        for (unsigned int dof = 0; dof < fe.dofs_per_hex; ++dof)
 +	          cell->set_mg_dof_index (cell->level (), dof, next_free_dof++);
 +
 +	      cell->set_user_flag ();
 +	      return next_free_dof;
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +	      return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +	      return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) {
 +	      return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<1>) {
 +	      return mg_faces.lines.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<2>) {
 +	      return mg_faces.quads.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    types::global_dof_index get_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const int2type<3>) {
 +	      return mg_level.dof_object.get_dof_index (dof_handler, obj_index, fe_index, local_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<1, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<1>& mg_level, internal::DoFHandler::DoFFaces<1>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +	      mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>&, internal::DoFHandler::DoFFaces<2>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +	      mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<2, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<2>& mg_level, internal::DoFHandler::DoFFaces<2>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) {
 +	      mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<1>) {
 +	      mg_faces.lines.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>&, internal::DoFHandler::DoFFaces<3>& mg_faces, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<2>) {
 +	      mg_faces.quads.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
 +
 +	    template<int spacedim>
 +	    static
 +	    void set_dof_index (const DoFHandler<3, spacedim>& dof_handler, internal::DoFHandler::DoFLevel<3>& mg_level, internal::DoFHandler::DoFFaces<3>&, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index, const int2type<3>) {
 +	      mg_level.dof_object.set_dof_index (dof_handler, obj_index, fe_index, local_index, global_index);
 +	    }
      };
    }
  }
@@@ -1571,13 -565,17 +895,17 @@@ types::global_dof_index DoFHandler<dim,
    std::set<int> boundary_dofs;
  
    const unsigned int dofs_per_face = get_fe().dofs_per_face;
 -  std::vector<unsigned int> dofs_on_face(dofs_per_face);
 +  std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
  
-                                    // loop over all faces to check
-                                    // whether they are at a
-                                    // boundary. note that we need not
-                                    // take special care of single
-                                    // lines (using
+                                    // loop over all faces of all cells
+                                    // and see whether they are at a
+                                    // boundary. note (i) that we visit
+                                    // interior faces twice (which we
+                                    // don't care about) but exterior
+                                    // faces only once as is
+                                    // appropriate, and (ii) that we
+                                    // need not take special care of
+                                    // single lines (using
                                     // @p{cell->has_boundary_lines}),
                                     // since we do not support
                                     // boundaries of dimension dim-2,
@@@ -1607,17 -606,25 +936,25 @@@ DoFHandler<dim,spacedim>::n_boundary_do
    std::set<int> boundary_dofs;
  
    const unsigned int dofs_per_face = get_fe().dofs_per_face;
 -  std::vector<unsigned int> dofs_on_face(dofs_per_face);
 +  std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
-   active_face_iterator face = begin_active_face (),
- 		               endf = end_face();
-   for (; face!=endf; ++face)
-     if (boundary_indicators.find(face->boundary_indicator()) !=
- 	    boundary_indicators.end())
-       {
- 	    face->get_dof_indices (dofs_on_face);
- 	    for (unsigned int i=0; i<dofs_per_face; ++i)
- 	      boundary_dofs.insert(dofs_on_face[i]);
-       }
+ 
+                                    // same as in the previous
+                                    // function, but with an additional
+                                    // check for the boundary indicator
+   active_cell_iterator cell = begin_active (),
+ 		       endc = end();
+   for (; cell!=endc; ++cell)
+     for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if (cell->at_boundary(f)
+ 	  &&
+ 	  (boundary_indicators.find(cell->face(f)->boundary_indicator()) !=
+ 	   boundary_indicators.end()))
+ 	{
+ 	  cell->face(f)->get_dof_indices (dofs_on_face);
+ 	  for (unsigned int i=0; i<dofs_per_face; ++i)
+ 	    boundary_dofs.insert(dofs_on_face[i]);
+ 	}
+ 
    return boundary_dofs.size();
  }
  
@@@ -1633,19 -640,27 +970,27 @@@ DoFHandler<dim,spacedim>::n_boundary_do
    std::set<int> boundary_dofs;
  
    const unsigned int dofs_per_face = get_fe().dofs_per_face;
 -  std::vector<unsigned int> dofs_on_face(dofs_per_face);
 +  std::vector<types::global_dof_index> dofs_on_face(dofs_per_face);
-   active_face_iterator face = begin_active_face (),
- 		               endf = end_face();
-   for (; face!=endf; ++face)
-     if (std::find (boundary_indicators.begin(),
- 		           boundary_indicators.end(),
- 		           face->boundary_indicator()) !=
- 	    boundary_indicators.end())
-       {
- 	    face->get_dof_indices (dofs_on_face);
- 	    for (unsigned int i=0; i<dofs_per_face; ++i)
- 	      boundary_dofs.insert(dofs_on_face[i]);
-       }
+ 
+                                    // same as in the previous
+                                    // function, but with a different
+                                    // check for the boundary indicator
+   active_cell_iterator cell = begin_active (),
+ 		       endc = end();
+   for (; cell!=endc; ++cell)
+     for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
+       if (cell->at_boundary(f)
+ 	  &&
+ 	  (std::find (boundary_indicators.begin(),
+                    boundary_indicators.end(),
+ 		      cell->face(f)->boundary_indicator()) !=
+ 	   boundary_indicators.end()))
+ 	{
+ 	  cell->face(f)->get_dof_indices (dofs_on_face);
+ 	  for (unsigned int i=0; i<dofs_per_face; ++i)
+ 	    boundary_dofs.insert(dofs_on_face[i]);
+ 	}
+ 
    return boundary_dofs.size();
  }
  
@@@ -1882,76 -840,12 +1227,78 @@@ void DoFHandler<dim,spacedim>::clear_sp
    delete faces;
    faces = 0;
  
 -  std::vector<unsigned int> tmp;
 +  std::vector<types::global_dof_index> tmp;
    std::swap (vertex_dofs, tmp);
+ 
+   number_cache.clear ();
  }
  
 +template<int dim, int spacedim>
 +template<int structdim>
 +types::global_dof_index DoFHandler<dim, spacedim>::get_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index) const {
 +  return internal::DoFHandler::Implementation::get_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, internal::int2type<structdim> ());
 +}
 +
 +template<int dim, int spacedim>
 +template<int structdim>
 +void DoFHandler<dim, spacedim>::set_dof_index (const unsigned int obj_level, const unsigned int obj_index, const unsigned int fe_index, const unsigned int local_index, const types::global_dof_index global_index) const {
 +  internal::DoFHandler::Implementation::set_dof_index (*this, *this->mg_levels[obj_level], *this->mg_faces, obj_index, fe_index, local_index, global_index, internal::int2type<structdim> ());
 +}
 +
 +template<int dim, int spacedim>
 +DoFHandler<dim, spacedim>::MGVertexDoFs::MGVertexDoFs (): coarsest_level (numbers::invalid_unsigned_int), finest_level (0), indices (0), indices_offset (0) {
 +}
 +
 +template<int dim, int spacedim>
 +DoFHandler<dim, spacedim>::MGVertexDoFs::~MGVertexDoFs () {
 +  delete[] indices;
 +  delete[] indices_offset;
 +}
 +
 +template<int dim, int spacedim>
 +void DoFHandler<dim, spacedim>::MGVertexDoFs::init (const unsigned int cl, const unsigned int fl, const unsigned int dofs_per_vertex) {
 +  if (indices != 0) {
 +    delete[] indices;
 +    indices = 0;
 +  }
 +
 +  if (indices_offset != 0) {
 +    delete[] indices_offset;
 +    indices_offset = 0;
 +  }
 +
 +  coarsest_level = cl;
 +  finest_level = fl;
 +
 +  if (cl > fl)
 +    return;
 +
 +  const unsigned int n_levels = finest_level - coarsest_level + 1;
 +  const unsigned int n_indices = n_levels * dofs_per_vertex;
 +
 +  indices = new unsigned int[n_indices];
 +  Assert (indices != 0, ExcNoMemory ());
 +
 +  for (unsigned int i = 0; i < n_indices; ++i)
 +    indices[i] = DoFHandler<dim, spacedim>::invalid_dof_index;
 +
 +  indices_offset = new unsigned int[n_levels];
 +  Assert (indices != 0, ExcNoMemory ());
 +
 +  for (unsigned int i = 0; i < n_levels; ++i)
 +    indices_offset[i] = i * dofs_per_vertex;
 +}
 +
 +template<int dim, int spacedim>
 +unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_coarsest_level () const {
 +  return coarsest_level;
 +}
 +
 +template<int dim, int spacedim>
 +unsigned int DoFHandler<dim, spacedim>::MGVertexDoFs::get_finest_level () const {
 +  return finest_level;
 +}
 +
  
  /*-------------- Explicit Instantiations -------------------------------*/
  #include "dof_handler.inst"
diff --cc deal.II/source/dofs/dof_renumbering.cc
index e969f3aabd,05f978abeb..f3edb6c3ab
--- a/deal.II/source/dofs/dof_renumbering.cc
+++ b/deal.II/source/dofs/dof_renumbering.cc
@@@ -15,7 -15,7 +15,8 @@@
  #include <deal.II/base/thread_management.h>
  #include <deal.II/base/utilities.h>
  #include <deal.II/base/quadrature_lib.h>
 +#include <deal.II/base/types.h>
+ #include <deal.II/base/template_constraints.h>
  
  #include <deal.II/lac/sparsity_pattern.h>
  #include <deal.II/lac/sparsity_tools.h>
@@@ -128,8 -128,7 +129,7 @@@ namespace DoFRenumberin
  
    namespace boost
    {
- #ifndef DEAL_II_BOOST_GRAPH_COMPILER_BUG
 -    namespace types
 +    namespace boosttypes
      {
        using namespace ::boost;
        using namespace std;
@@@ -205,19 -203,9 +204,9 @@@
                             const bool       reversed_numbering,
                             const bool       use_constraints)
      {
- #ifdef DEAL_II_BOOST_GRAPH_COMPILER_BUG
-       (void)new_dof_indices;
-       (void)dof_handler;
-       (void)reversed_numbering;
-       (void)use_constraints;
-       Assert (false,
-               ExcMessage ("Due to a bug in your compiler, this function triggers an internal "
-                           "compiler error and has been disabled. If you need to use the "
-                           "function, you need to upgrade to a newer version of the compiler."));
- #else
 -      types::Graph
 +      boosttypes::Graph
          graph(dof_handler.n_dofs());
 -      types::property_map<types::Graph,types::vertex_degree_t>::type
 +      boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
          graph_degree;
  
        internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
@@@ -273,19 -260,9 +261,9 @@@
                             const bool       reversed_numbering,
                             const bool       use_constraints)
      {
- #ifdef DEAL_II_BOOST_GRAPH_COMPILER_BUG
-       (void)new_dof_indices;
-       (void)dof_handler;
-       (void)reversed_numbering;
-       (void)use_constraints;
-       Assert (false,
-               ExcMessage ("Due to a bug in your compiler, this function triggers an internal "
-                           "compiler error and has been disabled. If you need to use the "
-                           "function, you need to upgrade to a newer version of the compiler."));
- #else
 -      types::Graph
 +      boosttypes::Graph
          graph(dof_handler.n_dofs());
 -      types::property_map<types::Graph,types::vertex_degree_t>::type
 +      boosttypes::property_map<boosttypes::Graph,boosttypes::vertex_degree_t>::type
          graph_degree;
  
        internal::create_graph (dof_handler, use_constraints, graph, graph_degree);
diff --cc deal.II/source/dofs/dof_tools.inst.in
index 9d6f34d4f8,b0f889f3f6..ba94bf15cd
--- a/deal.II/source/dofs/dof_tools.inst.in
+++ b/deal.II/source/dofs/dof_tools.inst.in
@@@ -776,6 -764,6 +776,14 @@@ DoFTools::extract_level_dofs<MGDoFHandl
   const BlockMask&,
   std::vector<bool>&);
  
++template
++void
++DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension, deal_II_dimension+1> >
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
++ const BlockMask&,
++ std::vector<bool>&);
++
  #endif
  
  
@@@ -800,9 -788,9 +808,17 @@@ DoFTools::extract_level_dofs<DoFHandler
   const ComponentMask &,
   std::vector<bool>&);
  
+ template
+ void
 -DoFTools::extract_level_dofs<1,3>
++DoFTools::extract_level_dofs<DoFHandler<1,3> >
++(const unsigned int level,
++ const DoFHandler<1,3>&,
++ const BlockMask&,
++ std::vector<bool>&);
++
 +template
 +void
 +DoFTools::extract_level_dofs<MGDoFHandler<1, 3> >
  (const unsigned int level,
   const MGDoFHandler<1,3>&,
   const BlockMask &,
diff --cc deal.II/source/grid/grid_tools.cc
index d67815fbe8,321023958a..d5a2255cc8
--- a/deal.II/source/grid/grid_tools.cc
+++ b/deal.II/source/grid/grid_tools.cc
@@@ -717,128 -701,128 +701,128 @@@ namespace GridTool
          for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_cell; v++)
            if (cell->vertex_index(v) == vertex)
              {
--              // OK, we found a cell that contains
--              // the given vertex. We add it
--              // to the list.
++					       // OK, we found a cell that contains
++					       // the given vertex. We add it
++					       // to the list.
                adjacent_cells.insert(cell);
  
--              // as explained above, in 2+d we need to check whether
--              // this vertex is on a face behind which there is a
--              // (possibly) coarser neighbor. if this is the case,
--              // then we need to also add this neighbor
++					       // as explained above, in 2+d we need to check whether
++					       // this vertex is on a face behind which there is a
++					       // (possibly) coarser neighbor. if this is the case,
++					       // then we need to also add this neighbor
                if (dim >= 2)
                  for (unsigned int vface = 0; vface < dim; vface++)
                    {
                      const unsigned int face =
--                        GeometryInfo<dim>::vertex_to_face[v][vface];
++		      GeometryInfo<dim>::vertex_to_face[v][vface];
  
                      if (!cell->at_boundary(face)
                          &&
                          cell->neighbor(face)->active())
                        {
--                        // there is a (possibly) coarser cell behind a
--                        // face to which the vertex belongs. the
--                        // vertex we are looking at is then either a
--                        // vertex of that coarser neighbor, or it is a
--                        // hanging node on one of the faces of that
--                        // cell. in either case, it is adjacent to the
--                        // vertex, so add it to the list as well (if
--                        // the cell was already in the list then the
--                        // std::set makes sure that we get it only
--                        // once)
++							 // there is a (possibly) coarser cell behind a
++							 // face to which the vertex belongs. the
++							 // vertex we are looking at is then either a
++							 // vertex of that coarser neighbor, or it is a
++							 // hanging node on one of the faces of that
++							 // cell. in either case, it is adjacent to the
++							 // vertex, so add it to the list as well (if
++							 // the cell was already in the list then the
++							 // std::set makes sure that we get it only
++							 // once)
                          adjacent_cells.insert (cell->neighbor(face));
                        }
                    }
  
--              // in any case, we have found a cell, so go to the next cell
++					       // in any case, we have found a cell, so go to the next cell
                goto next_cell;
              }
  
--        // in 3d also loop over the edges
++					 // in 3d also loop over the edges
          if (dim >= 3)
            {
              for (unsigned int e=0; e<GeometryInfo<dim>::lines_per_cell; ++e)
                if (cell->line(e)->has_children())
--                // the only place where this vertex could have been
--                // hiding is on the mid-edge point of the edge we
--                // are looking at
++						 // the only place where this vertex could have been
++						 // hiding is on the mid-edge point of the edge we
++						 // are looking at
                  if (cell->line(e)->child(0)->vertex_index(1) == vertex)
                    {
                      adjacent_cells.insert(cell);
  
--                    // jump out of this tangle of nested loops
++						     // jump out of this tangle of nested loops
                      goto next_cell;
                    }
            }
  
--        // in more than 3d we would probably have to do the same as
--        // above also for even lower-dimensional objects
++					 // in more than 3d we would probably have to do the same as
++					 // above also for even lower-dimensional objects
          Assert (dim <= 3, ExcNotImplemented());
  
--        // move on to the next cell if we have found the
--        // vertex on the current one
++					 // move on to the next cell if we have found the
++					 // vertex on the current one
          next_cell:
          ;
        }
  
--    // if this was an active vertex then there needs to have been
--    // at least one cell to which it is adjacent!
++				     // if this was an active vertex then there needs to have been
++				     // at least one cell to which it is adjacent!
      Assert (adjacent_cells.size() > 0, ExcInternalError());
  
--    // return the result as a vector, rather than the set we built above
++				     // return the result as a vector, rather than the set we built above
      return
--    std::vector<typename Container<dim,spacedim>::active_cell_iterator>
--    (adjacent_cells.begin(), adjacent_cells.end());
++      std::vector<typename Container<dim,spacedim>::active_cell_iterator>
++      (adjacent_cells.begin(), adjacent_cells.end());
    }
  
  
  
    namespace
    {
--	template <int dim, template<int, int> class Container, int spacedim>
--	void find_active_cell_around_point_internal(const Container<dim,spacedim>& container,
--	  std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells,
--	  std::set<typename Container<dim,spacedim>::active_cell_iterator>& adjacent_cells)
++    template <int dim, template<int, int> class Container, int spacedim>
++    void find_active_cell_around_point_internal(const Container<dim,spacedim>& container,
++						std::set<typename Container<dim,spacedim>::active_cell_iterator>& searched_cells,
++						std::set<typename Container<dim,spacedim>::active_cell_iterator>& adjacent_cells)
++    {
++      typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
++
++				       // update the searched cells
++      searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
++				       // now we to collect all neighbors
++				       // of the cells in adjacent_cells we
++				       // have not yet searched.
++      std::set<cell_iterator> adjacent_cells_new;
++
++      typename std::set<cell_iterator>::const_iterator
++	cell = adjacent_cells.begin(),
++	endc = adjacent_cells.end();
++      for(; cell != endc; ++cell)
  	{
--	  typedef typename Container<dim,spacedim>::active_cell_iterator cell_iterator;
--
--					// update the searched cells
--	  searched_cells.insert(adjacent_cells.begin(), adjacent_cells.end());
--					// now we to collect all neighbors
--					// of the cells in adjacent_cells we
--					// have not yet searched.
--	  std::set<cell_iterator> adjacent_cells_new;
--
--	  typename std::set<cell_iterator>::const_iterator
--	  cell = adjacent_cells.begin(),
--	  endc = adjacent_cells.end();
--	  for(; cell != endc; ++cell)
--	  {
--		std::vector<cell_iterator> active_neighbors;
--		get_active_neighbors<Container<dim, spacedim> >(*cell, active_neighbors);
--		for (unsigned int i=0; i<active_neighbors.size(); ++i)
--		  if(searched_cells.find(active_neighbors[i]) == searched_cells.end())
--			adjacent_cells_new.insert(active_neighbors[i]);
--	  }
--	  adjacent_cells.clear();
--	  adjacent_cells.insert(adjacent_cells_new.begin(), adjacent_cells_new.end());
--	  if (adjacent_cells.size() == 0)
--      {
--              // we haven't found any other cell that would be a
--              // neighbor of a previously found cell, but we know
--              // that we haven't checked all cells yet. that means
--              // that the domain is disconnected. in that case,
--              // choose the first previously untouched cell we
--              // can find
--        cell_iterator it = container.begin_active();
--		for ( ; it!=container.end();++it)
--		  if(searched_cells.find(it) == searched_cells.end())
--		  {
--			adjacent_cells.insert(it);
--			 break;
--		  }
--       }
++	  std::vector<cell_iterator> active_neighbors;
++	  get_active_neighbors<Container<dim, spacedim> >(*cell, active_neighbors);
++	  for (unsigned int i=0; i<active_neighbors.size(); ++i)
++	    if(searched_cells.find(active_neighbors[i]) == searched_cells.end())
++	      adjacent_cells_new.insert(active_neighbors[i]);
++	}
++      adjacent_cells.clear();
++      adjacent_cells.insert(adjacent_cells_new.begin(), adjacent_cells_new.end());
++      if (adjacent_cells.size() == 0)
++	{
++					   // we haven't found any other cell that would be a
++					   // neighbor of a previously found cell, but we know
++					   // that we haven't checked all cells yet. that means
++					   // that the domain is disconnected. in that case,
++					   // choose the first previously untouched cell we
++					   // can find
++	  cell_iterator it = container.begin_active();
++	  for ( ; it!=container.end();++it)
++	    if(searched_cells.find(it) == searched_cells.end())
++	      {
++		adjacent_cells.insert(it);
++		break;
++	      }
  	}
++    }
    }
  
    template <int dim, template<int, int> class Container, int spacedim>
@@@ -874,80 -858,80 +858,80 @@@
      std::vector<cell_iterator> adjacent_cells_tmp =
        find_cells_adjacent_to_vertex(container, vertex);
  
--									// Make sure that we have found
--									// at least one cell adjacent to vertex.
--	  Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
--
--									// Copy all the cells into a std::set
--	 std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
--	 std::set<cell_iterator> searched_cells;
--
--									// Determine the maximal number of cells
--									// in the grid.
--									// As long as we have not found
--									// the cell and have not searched
--									// every cell in the triangulation,
--									// we keep on looking.
--	  const unsigned int n_cells =get_tria(container).n_cells();
--	  bool found = false;
--	  unsigned int cells_searched = 0;
--	  while (!found && cells_searched < n_cells)
++				     // Make sure that we have found
++				     // at least one cell adjacent to vertex.
++    Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
++
++				     // Copy all the cells into a std::set
++    std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
++    std::set<cell_iterator> searched_cells;
++
++				     // Determine the maximal number of cells
++				     // in the grid.
++				     // As long as we have not found
++				     // the cell and have not searched
++				     // every cell in the triangulation,
++				     // we keep on looking.
++    const unsigned int n_cells =get_tria(container).n_cells();
++    bool found = false;
++    unsigned int cells_searched = 0;
++    while (!found && cells_searched < n_cells)
++      {
++	typename std::set<cell_iterator>::const_iterator
++	  cell = adjacent_cells.begin(),
++	  endc = adjacent_cells.end();
++	for(; cell != endc; ++cell)
  	  {
--	    typename std::set<cell_iterator>::const_iterator
--		cell = adjacent_cells.begin(),
--		endc = adjacent_cells.end();
--		for(; cell != endc; ++cell)
--		  {
--			try
--			  {
--				const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
--
--												// calculate the infinity norm of
--												// the distance vector to the unit cell.
--				const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
--
--												// We compare if the point is inside the
--												// unit cell (or at least not too far
--												// outside). If it is, it is also checked
--												// that the cell has a more refined state
--				if (dist < best_distance ||
--					(dist == best_distance && (*cell)->level() > best_level))
--				  {
--					found 		  = true;
--					best_distance = dist;
--					best_level    = (*cell)->level();
--					best_cell     = std::make_pair(*cell, p_cell);
--				  }
--			  }
--			catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
--			  {
--							// ok, the transformation
--							// failed presumably
--							// because the point we
--							// are looking for lies
--							// outside the current
--							// cell. this means that
--							// the current cell can't
--							// be the cell around the
--							// point, so just ignore
--							// this cell and move on
--							// to the next
--			  }
--		  }
--							//udpate the number of cells searched
--		  cells_searched += adjacent_cells.size();
--							// if we have not found the cell in
--							// question and have not yet searched every
--							// cell, we expand our search to
--							// all the not already searched neighbors of
--							// the cells in adjacent_cells. This is
--							// what find_active_cell_around_poin_internal
--							// is for.
--		  if(!found && cells_searched < n_cells)
++	    try
++	      {
++		const Point<dim> p_cell = mapping.transform_real_to_unit_cell(*cell, p);
++
++						 // calculate the infinity norm of
++						 // the distance vector to the unit cell.
++		const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
++
++						 // We compare if the point is inside the
++						 // unit cell (or at least not too far
++						 // outside). If it is, it is also checked
++						 // that the cell has a more refined state
++		if (dist < best_distance ||
++		    (dist == best_distance && (*cell)->level() > best_level))
  		  {
--			find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++		    found 		  = true;
++		    best_distance = dist;
++		    best_level    = (*cell)->level();
++		    best_cell     = std::make_pair(*cell, p_cell);
  		  }
++	      }
++	    catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
++	      {
++						 // ok, the transformation
++						 // failed presumably
++						 // because the point we
++						 // are looking for lies
++						 // outside the current
++						 // cell. this means that
++						 // the current cell can't
++						 // be the cell around the
++						 // point, so just ignore
++						 // this cell and move on
++						 // to the next
++	      }
++	  }
++					 //udpate the number of cells searched
++	cells_searched += adjacent_cells.size();
++					 // if we have not found the cell in
++					 // question and have not yet searched every
++					 // cell, we expand our search to
++					 // all the not already searched neighbors of
++					 // the cells in adjacent_cells. This is
++					 // what find_active_cell_around_poin_internal
++					 // is for.
++	if(!found && cells_searched < n_cells)
++	  {
++	    find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
  	  }
++      }
  
      Assert (best_cell.first.state() == IteratorState::valid,
              ExcPointNotFound<spacedim>(p));
@@@ -993,84 -977,84 +977,84 @@@
                                           // all adjacent cells
          unsigned int vertex = find_closest_vertex(container, p);
  
--		std::vector<cell_iterator> adjacent_cells_tmp =
--        find_cells_adjacent_to_vertex(container, vertex);
--
--									// Make sure that we have found
--									// at least one cell adjacent to vertex.
--	  Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
--
--									// Copy all the cells into a std::set
--	 std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
--	 std::set<cell_iterator> searched_cells;
--
--									// Determine the maximal number of cells
--									// in the grid.
--									// As long as we have not found
--									// the cell and have not searched
--									// every cell in the triangulation,
--									// we keep on looking.
--	  const unsigned int n_cells =get_tria(container).n_cells();
--	  bool found = false;
--	  unsigned int cells_searched = 0;
--	  while (!found && cells_searched < n_cells)
++	std::vector<cell_iterator> adjacent_cells_tmp =
++	  find_cells_adjacent_to_vertex(container, vertex);
++
++					 // Make sure that we have found
++					 // at least one cell adjacent to vertex.
++	Assert(adjacent_cells_tmp.size()>0, ExcInternalError());
++
++					 // Copy all the cells into a std::set
++	std::set<cell_iterator> adjacent_cells(adjacent_cells_tmp.begin(), adjacent_cells_tmp.end());
++	std::set<cell_iterator> searched_cells;
++
++					 // Determine the maximal number of cells
++					 // in the grid.
++					 // As long as we have not found
++					 // the cell and have not searched
++					 // every cell in the triangulation,
++					 // we keep on looking.
++	const unsigned int n_cells =get_tria(container).n_cells();
++	bool found = false;
++	unsigned int cells_searched = 0;
++	while (!found && cells_searched < n_cells)
  	  {
  	    typename std::set<cell_iterator>::const_iterator
--		cell = adjacent_cells.begin(),
--		endc = adjacent_cells.end();
--		for(; cell != endc; ++cell)
++	      cell = adjacent_cells.begin(),
++	      endc = adjacent_cells.end();
++	    for(; cell != endc; ++cell)
++	      {
++		try
  		  {
--			try
--			  {
--				const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
--
--
--				// calculate the infinity norm of
--				// the distance vector to the unit cell.
--				const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
--
--												// We compare if the point is inside the
--												// unit cell (or at least not too far
--												// outside). If it is, it is also checked
--												// that the cell has a more refined state
--				if (dist < best_distance ||
--					(dist == best_distance && (*cell)->level() > best_level))
--				  {
--					found 		  = true;
--					best_distance = dist;
--					best_level    = (*cell)->level();
--					best_cell     = std::make_pair(*cell, p_cell);
--				  }
--			  }
--			catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
--			  {
--							// ok, the transformation
--							// failed presumably
--							// because the point we
--							// are looking for lies
--							// outside the current
--							// cell. this means that
--							// the current cell can't
--							// be the cell around the
--							// point, so just ignore
--							// this cell and move on
--							// to the next
--			  }
++		    const Point<dim> p_cell = mapping[(*cell)->active_fe_index()].transform_real_to_unit_cell(*cell, p);
++
++
++						     // calculate the infinity norm of
++						     // the distance vector to the unit cell.
++		    const double dist = GeometryInfo<dim>::distance_to_unit_cell(p_cell);
++
++						     // We compare if the point is inside the
++						     // unit cell (or at least not too far
++						     // outside). If it is, it is also checked
++						     // that the cell has a more refined state
++		    if (dist < best_distance ||
++			(dist == best_distance && (*cell)->level() > best_level))
++		      {
++			found 		  = true;
++			best_distance = dist;
++			best_level    = (*cell)->level();
++			best_cell     = std::make_pair(*cell, p_cell);
++		      }
  		  }
--		  	  							//udpate the number of cells searched
--		  cells_searched += adjacent_cells.size();
--							// if we have not found the cell in
--							// question and have not yet searched every
--							// cell, we expand our search to
--							// all the not already searched neighbors of
--							// the cells in adjacent_cells.
--		  if(!found && cells_searched < n_cells)
++		catch (typename MappingQ1<dim,spacedim>::ExcTransformationFailed &)
  		  {
--			find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++						     // ok, the transformation
++						     // failed presumably
++						     // because the point we
++						     // are looking for lies
++						     // outside the current
++						     // cell. this means that
++						     // the current cell can't
++						     // be the cell around the
++						     // point, so just ignore
++						     // this cell and move on
++						     // to the next
  		  }
++	      }
++					     //udpate the number of cells searched
++	    cells_searched += adjacent_cells.size();
++					     // if we have not found the cell in
++					     // question and have not yet searched every
++					     // cell, we expand our search to
++					     // all the not already searched neighbors of
++					     // the cells in adjacent_cells.
++	    if(!found && cells_searched < n_cells)
++	      {
++		find_active_cell_around_point_internal(container, searched_cells, adjacent_cells);
++	      }
  
  	  }
--	  }
++      }
  
      Assert (best_cell.first.state() == IteratorState::valid,
              ExcPointNotFound<spacedim>(p));
@@@ -2177,45 -2161,45 +2161,45 @@@
                                     const CellIterator &cell2,
                                     const int          direction,
                                     const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
--                                                      &offset)
++				   &offset)
    {
--                                   // An orthogonal equality test for
--                                   // cells:
--                                   // Two cells are equal if the
--                                   // corresponding vertices of the
--                                   // boundary faces can be transformed
--                                   // into each other parallel to the
--                                   // given direction.
--                                   // It is assumed that
--                                   // cell1->face(2*direction) and
--                                   // cell2->face(2*direction+1) are the
--                                   // corresponding boundary faces.
--                                   // The optional argument offset will
--                                   // be added to each vertex of cell1
++				     // An orthogonal equality test for
++				     // cells:
++				     // Two cells are equal if the
++				     // corresponding vertices of the
++				     // boundary faces can be transformed
++				     // into each other parallel to the
++				     // given direction.
++				     // It is assumed that
++				     // cell1->face(2*direction) and
++				     // cell2->face(2*direction+1) are the
++				     // corresponding boundary faces.
++				     // The optional argument offset will
++				     // be added to each vertex of cell1
  
      static const int dim = CellIterator::AccessorType::dimension;
      static const int space_dim = CellIterator::AccessorType::space_dimension;
      Assert(dim == space_dim,
             ExcNotImplemented());
--                                   // ... otherwise we would need two
--                                   // directions: One for selecting the
--                                   // faces, thus living in dim,
--                                   // the other direction for selecting the
--                                   // direction for the orthogonal
--                                   // projection, thus living in
--                                   // space_dim.
++				     // ... otherwise we would need two
++				     // directions: One for selecting the
++				     // faces, thus living in dim,
++				     // the other direction for selecting the
++				     // direction for the orthogonal
++				     // projection, thus living in
++				     // space_dim.
  
      for (int i = 0; i < space_dim; ++i) {
--                                   // Only compare
--                                   // coordinate-components != direction:
++				       // Only compare
++				       // coordinate-components != direction:
        if (i == direction)
          continue;
  
        for (unsigned int v = 0; v < GeometryInfo<dim>::vertices_per_face; ++v)
          if (fabs (  cell1->face(2*direction)->vertex(v)(i)
--                  + offset[i]
--                  - cell2->face(2*direction+1)->vertex(v)(i))
--                    > 1.e-10)
++		    + offset[i]
++		    - cell2->face(2*direction+1)->vertex(v)(i))
++	    > 1.e-10)
            return false;
      }
      return true;
@@@ -2229,7 -2213,7 +2213,7 @@@
                                 const types::boundary_id                  boundary_component,
                                 int                                         direction,
                                 const dealii::Tensor<1,CellIterator::AccessorType::space_dimension>
--                                 &offset)
++			       &offset)
    {
      static const int space_dim = CellIterator::AccessorType::space_dimension;
      Assert (0<=direction && direction<space_dim,
@@@ -2238,11 -2222,11 +2222,11 @@@
      std::set<CellIterator> cells1;
      std::set<CellIterator> cells2;
  
--                                   // collect boundary cells, i.e. cells
--                                   // with face(2*direction), resp.
--                                   // face(2*direction+1) being on the
--                                   // boundary and having the correct
--                                   // color:
++				     // collect boundary cells, i.e. cells
++				     // with face(2*direction), resp.
++				     // face(2*direction+1) being on the
++				     // boundary and having the correct
++				     // color:
      CellIterator cell = begin;
      for (; cell!= end; ++cell) {
        if (cell->face(2*direction)->at_boundary() &&
@@@ -2260,16 -2244,16 +2244,16 @@@
  
      std::map<CellIterator, CellIterator> matched_cells;
  
--                                   // Match with a complexity of O(n^2).
--                                   // This could be improved...
++				     // Match with a complexity of O(n^2).
++				     // This could be improved...
      typedef typename std::set<CellIterator>::const_iterator SetIterator;
      for (SetIterator it1 = cells1.begin(); it1 != cells1.end(); ++it1) {
        for (SetIterator it2 = cells2.begin(); it2 != cells2.end(); ++it2) {
          if (orthogonal_equality(*it1, *it2, direction, offset)) {
--                                   // We have a match, so insert the
--                                   // matching pairs and remove the matched
--                                   // cell in cells2 to speed up the
--                                   // matching:
++					   // We have a match, so insert the
++					   // matching pairs and remove the matched
++					   // cell in cells2 to speed up the
++					   // matching:
            matched_cells[*it1] = *it2;
            cells2.erase(it2);
            break;
@@@ -2291,14 -2275,14 +2275,14 @@@
                                 const types::boundary_id boundary_component,
                                 int                        direction,
                                 const dealii::Tensor<1,DH::space_dimension>
--                                 &offset)
++			       &offset)
    {
      return collect_periodic_cell_pairs<typename DH::cell_iterator>
--                                      (dof_handler.begin_active(),
--                                       dof_handler.end(),
--                                       boundary_component,
--                                       direction,
--                                       offset);
++      (dof_handler.begin_active(),
++       dof_handler.end(),
++       boundary_component,
++       direction,
++       offset);
    }
  
  
diff --cc deal.II/source/hp/dof_handler.cc
index 41abf8114c,1e162dcfb9..8847f900c1
--- a/deal.II/source/hp/dof_handler.cc
+++ b/deal.II/source/hp/dof_handler.cc
@@@ -3123,8 -2142,17 +2143,17 @@@ types::global_dof_index DoFHandler<1,3>
    template<int dim, int spacedim>
    void
    DoFHandler<dim,spacedim>::
 -  compute_line_dof_identities (std::vector<unsigned int> &new_dof_indices) const
 +  compute_line_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
    {
+ 				     // we will mark lines that we
+ 				     // have already treated, so first
+ 				     // save and clear the user flags
+ 				     // on lines and later restore
+ 				     // them
+     std::vector<bool> user_flags;
+     this->get_tria().save_user_flags_line(user_flags);
+     const_cast<Triangulation<dim,spacedim> &>(this->get_tria()).clear_user_flags_line ();
+ 
                                       // An implementation of the
                                       // algorithm described in the hp
                                       // paper, including the
@@@ -3161,176 -2189,185 +2190,185 @@@
                                       // at most 2
      Table<2,std_cxx1x::shared_ptr<internal::hp::DoFIdentities> >
        line_dof_identities (finite_elements->size(),
 -                           finite_elements->size());
 +			               finite_elements->size());
  
-     for (line_iterator line=begin_line(); line!=end_line(); ++line)
-       {
+     for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
+       for (unsigned int l=0; l<GeometryInfo<dim>::lines_per_cell; ++l)
+ 	if (cell->line(l)->user_flag_set() == false)
+ 	  {
+ 	    const line_iterator line = cell->line(l);
+ 	    line->set_user_flag ();
+ 
  	    unsigned int unique_sets_of_dofs
  	      = line->n_active_fe_indices();
  
-                                          // do a first loop over all sets of
-                                          // dofs and do identity
-                                          // uniquification
+ 					     // do a first loop over all sets of
+ 					     // dofs and do identity
+ 					     // uniquification
  	    for (unsigned int f=0; f<line->n_active_fe_indices(); ++f)
  	      for (unsigned int g=f+1; g<line->n_active_fe_indices(); ++g)
- 	        {
- 	          const unsigned int fe_index_1 = line->nth_active_fe_index (f),
- 				                 fe_index_2 = line->nth_active_fe_index (g);
- 
- 	          if (((*finite_elements)[fe_index_1].dofs_per_line
- 		           ==
- 		           (*finite_elements)[fe_index_2].dofs_per_line)
- 		          &&
- 		          ((*finite_elements)[fe_index_1].dofs_per_line > 0))
- 		        {
- 		          internal::hp::ensure_existence_of_dof_identities<1>
- 		            ((*finite_elements)[fe_index_1],
- 		             (*finite_elements)[fe_index_2],
- 		             line_dof_identities[fe_index_1][fe_index_2]);
-                                                    // see if these sets of dofs
-                                                    // are identical. the first
-                                                    // condition for this is that
-                                                    // indeed there are n
-                                                    // identities
- 		          if (line_dof_identities[fe_index_1][fe_index_2]->size()
- 		              ==
- 		              (*finite_elements)[fe_index_1].dofs_per_line)
- 		            {
- 		              unsigned int i=0;
- 		              for (; i<(*finite_elements)[fe_index_1].dofs_per_line; ++i)
- 			            if (((*(line_dof_identities[fe_index_1][fe_index_2]))[i].first != i)
- 			                &&
- 			                ((*(line_dof_identities[fe_index_1][fe_index_2]))[i].second != i))
-                                                            // not an identity
- 			              break;
- 
- 		              if (i == (*finite_elements)[fe_index_1].dofs_per_line)
- 			            {
-                                                            // the dofs of these
-                                                            // two finite
-                                                            // elements are
-                                                            // identical. as a
-                                                            // safety check,
-                                                            // ensure that none
-                                                            // of the two FEs is
-                                                            // trying to dominate
-                                                            // the other, which
-                                                            // wouldn't make any
-                                                            // sense in this case
- 			              Assert ((*finite_elements)[fe_index_1].compare_for_face_domination
- 				                  ((*finite_elements)[fe_index_2])
- 				                  ==
- 				                  FiniteElementDomination::either_element_can_dominate,
- 				                  ExcInternalError());
- 
- 			              --unique_sets_of_dofs;
- 
- 			              for (unsigned int j=0; j<(*finite_elements)[fe_index_1].dofs_per_line; ++j)
- 			                {
- 			                  const types::global_dof_index master_dof_index
- 				                = line->dof_index (j, fe_index_1);
- 			                  const types::global_dof_index slave_dof_index
- 				                = line->dof_index (j, fe_index_2);
- 
-                                                                // if master dof
-                                                                // was already
-                                                                // constrained,
-                                                                // constrain to
-                                                                // that one,
-                                                                // otherwise
-                                                                // constrain
-                                                                // slave to
-                                                                // master
- 			                  if (new_dof_indices[master_dof_index] !=
- 				                  numbers::invalid_unsigned_int)
- 				                {
- 				                  Assert (new_dof_indices[new_dof_indices[master_dof_index]] ==
- 					                      numbers::invalid_unsigned_int,
- 					                      ExcInternalError());
- 
- 				                  new_dof_indices[slave_dof_index]
- 				                    = new_dof_indices[master_dof_index];
- 				                }
- 			                  else
- 				                {
- 				                  Assert ((new_dof_indices[master_dof_index] ==
- 					                       numbers::invalid_unsigned_int)
- 					                      ||
- 					                      (new_dof_indices[slave_dof_index] ==
- 					                       master_dof_index),
- 					                      ExcInternalError());
- 
- 				                  new_dof_indices[slave_dof_index] = master_dof_index;
- 				                }
- 			                }
- 			            }
- 		            }
- 		        }
- 	        }
- 
-                                          // if at this point, there is only
-                                          // one unique set of dofs left, then
-                                          // we have taken care of everything
-                                          // above. if there are two, then we
-                                          // need to deal with them here. if
-                                          // there are more, then we punt, as
-                                          // described in the paper (and
-                                          // mentioned above)
+ 		{
+ 		  const unsigned int fe_index_1 = line->nth_active_fe_index (f),
+ 				     fe_index_2 = line->nth_active_fe_index (g);
+ 
+ 		  if (((*finite_elements)[fe_index_1].dofs_per_line
+ 		       ==
+ 		       (*finite_elements)[fe_index_2].dofs_per_line)
+ 		      &&
+ 		      ((*finite_elements)[fe_index_1].dofs_per_line > 0))
+ 		    {
+ 		      internal::hp::ensure_existence_of_dof_identities<1>
+ 			((*finite_elements)[fe_index_1],
+ 			 (*finite_elements)[fe_index_2],
+ 			 line_dof_identities[fe_index_1][fe_index_2]);
+ 						       // see if these sets of dofs
+ 						       // are identical. the first
+ 						       // condition for this is that
+ 						       // indeed there are n
+ 						       // identities
+ 		      if (line_dof_identities[fe_index_1][fe_index_2]->size()
+ 			  ==
+ 			  (*finite_elements)[fe_index_1].dofs_per_line)
+ 			{
+ 			  unsigned int i=0;
+ 			  for (; i<(*finite_elements)[fe_index_1].dofs_per_line; ++i)
+ 			    if (((*(line_dof_identities[fe_index_1][fe_index_2]))[i].first != i)
+ 				&&
+ 				((*(line_dof_identities[fe_index_1][fe_index_2]))[i].second != i))
+ 							       // not an identity
+ 			      break;
+ 
+ 			  if (i == (*finite_elements)[fe_index_1].dofs_per_line)
+ 			    {
+ 							       // the dofs of these
+ 							       // two finite
+ 							       // elements are
+ 							       // identical. as a
+ 							       // safety check,
+ 							       // ensure that none
+ 							       // of the two FEs is
+ 							       // trying to dominate
+ 							       // the other, which
+ 							       // wouldn't make any
+ 							       // sense in this case
+ 			      Assert ((*finite_elements)[fe_index_1].compare_for_face_domination
+ 				      ((*finite_elements)[fe_index_2])
+ 				      ==
+ 				      FiniteElementDomination::either_element_can_dominate,
+ 				      ExcInternalError());
+ 
+ 			      --unique_sets_of_dofs;
+ 
+ 			      for (unsigned int j=0; j<(*finite_elements)[fe_index_1].dofs_per_line; ++j)
+ 				{
 -				  const unsigned int master_dof_index
++				  const types::global_dof_index master_dof_index
+ 				    = line->dof_index (j, fe_index_1);
 -				  const unsigned int slave_dof_index
++				  const types::global_dof_index slave_dof_index
+ 				    = line->dof_index (j, fe_index_2);
+ 
+ 								   // if master dof
+ 								   // was already
+ 								   // constrained,
+ 								   // constrain to
+ 								   // that one,
+ 								   // otherwise
+ 								   // constrain
+ 								   // slave to
+ 								   // master
+ 				  if (new_dof_indices[master_dof_index] !=
+ 				      numbers::invalid_unsigned_int)
+ 				    {
+ 				      Assert (new_dof_indices[new_dof_indices[master_dof_index]] ==
+ 					      numbers::invalid_unsigned_int,
+ 					      ExcInternalError());
+ 
+ 				      new_dof_indices[slave_dof_index]
+ 					= new_dof_indices[master_dof_index];
+ 				    }
+ 				  else
+ 				    {
+ 				      Assert ((new_dof_indices[master_dof_index] ==
+ 					       numbers::invalid_unsigned_int)
+ 					      ||
+ 					      (new_dof_indices[slave_dof_index] ==
+ 					       master_dof_index),
+ 					      ExcInternalError());
+ 
+ 				      new_dof_indices[slave_dof_index] = master_dof_index;
+ 				    }
+ 				}
+ 			    }
+ 			}
+ 		    }
+ 		}
+ 
+ 					     // if at this point, there is only
+ 					     // one unique set of dofs left, then
+ 					     // we have taken care of everything
+ 					     // above. if there are two, then we
+ 					     // need to deal with them here. if
+ 					     // there are more, then we punt, as
+ 					     // described in the paper (and
+ 					     // mentioned above)
  //TODO: The check for 'dim==2' was inserted by intuition. It fixes
  // the previous problems with step-27 in 3D. But an explanation
  // for this is still required, and what we do here is not what we
  // describe in the paper!.
  	    if ((unique_sets_of_dofs == 2) && (dim == 2))
  	      {
-                                              // find out which is the
-                                              // most dominating finite
-                                              // element of the ones that
-                                              // are used on this line
- 	        const unsigned int most_dominating_fe_index
- 	          = internal::hp::get_most_dominating_fe_index<dim,spacedim> (line);
- 
- 	        const unsigned int n_active_fe_indices
- 	          = line->n_active_fe_indices ();
- 
-                                              // loop over the indices of
-                                              // all the finite elements
-                                              // that are not dominating,
-                                              // and identify their dofs
-                                              // to the most dominating
-                                              // one
- 	        for (unsigned int f=0; f<n_active_fe_indices; ++f)
- 	          if (line->nth_active_fe_index (f) !=
- 		          most_dominating_fe_index)
- 		        {
- 		          const unsigned int
- 		            other_fe_index = line->nth_active_fe_index (f);
- 
- 		          internal::hp::ensure_existence_of_dof_identities<1>
- 		            ((*finite_elements)[most_dominating_fe_index],
- 		             (*finite_elements)[other_fe_index],
- 		             line_dof_identities[most_dominating_fe_index][other_fe_index]);
- 
- 		          internal::hp::DoFIdentities &identities
- 		            = *line_dof_identities[most_dominating_fe_index][other_fe_index];
- 		          for (unsigned int i=0; i<identities.size(); ++i)
- 		            {
- 		              const types::global_dof_index master_dof_index
- 			            = line->dof_index (identities[i].first, most_dominating_fe_index);
- 		              const types::global_dof_index slave_dof_index
- 			            = line->dof_index (identities[i].second, other_fe_index);
- 
- 		              Assert ((new_dof_indices[master_dof_index] ==
- 			                   numbers::invalid_unsigned_int)
- 			                  ||
- 			                  (new_dof_indices[slave_dof_index] ==
- 			                   master_dof_index),
- 			                  ExcInternalError());
- 
- 		              new_dof_indices[slave_dof_index] = master_dof_index;
- 		            }
- 		        }
+ 						 // find out which is the
+ 						 // most dominating finite
+ 						 // element of the ones that
+ 						 // are used on this line
+ 		const unsigned int most_dominating_fe_index
+ 		  = internal::hp::get_most_dominating_fe_index<dim,spacedim> (line);
+ 
+ 		const unsigned int n_active_fe_indices
+ 		  = line->n_active_fe_indices ();
+ 
+ 						 // loop over the indices of
+ 						 // all the finite elements
+ 						 // that are not dominating,
+ 						 // and identify their dofs
+ 						 // to the most dominating
+ 						 // one
+ 		for (unsigned int f=0; f<n_active_fe_indices; ++f)
+ 		  if (line->nth_active_fe_index (f) !=
+ 		      most_dominating_fe_index)
+ 		    {
+ 		      const unsigned int
+ 			other_fe_index = line->nth_active_fe_index (f);
+ 
+ 		      internal::hp::ensure_existence_of_dof_identities<1>
+ 			((*finite_elements)[most_dominating_fe_index],
+ 			 (*finite_elements)[other_fe_index],
+ 			 line_dof_identities[most_dominating_fe_index][other_fe_index]);
+ 
+ 		      internal::hp::DoFIdentities &identities
+ 			= *line_dof_identities[most_dominating_fe_index][other_fe_index];
+ 		      for (unsigned int i=0; i<identities.size(); ++i)
+ 			{
 -			  const unsigned int master_dof_index
++			  const types::global_dof_index master_dof_index
+ 			    = line->dof_index (identities[i].first, most_dominating_fe_index);
 -			  const unsigned int slave_dof_index
++			  const types::global_dof_index slave_dof_index
+ 			    = line->dof_index (identities[i].second, other_fe_index);
+ 
+ 			  Assert ((new_dof_indices[master_dof_index] ==
+ 				   numbers::invalid_unsigned_int)
+ 				  ||
+ 				  (new_dof_indices[slave_dof_index] ==
+ 				   master_dof_index),
+ 				  ExcInternalError());
+ 
+ 			  new_dof_indices[slave_dof_index] = master_dof_index;
+ 			}
+ 		    }
  	      }
-       }
+ 	  }
+ 
+ 				     // finally restore the user flags
+     const_cast<Triangulation<dim,spacedim> &>(this->get_tria())
+       .load_user_flags_line(user_flags);
    }
  
  
@@@ -3372,8 -2409,17 +2410,17 @@@
    template<int dim, int spacedim>
    void
    DoFHandler<dim,spacedim>::
 -  compute_quad_dof_identities (std::vector<unsigned int> &new_dof_indices) const
 +  compute_quad_dof_identities (std::vector<types::global_dof_index> &new_dof_indices) const
    {
+ 				     // we will mark quads that we
+ 				     // have already treated, so first
+ 				     // save and clear the user flags
+ 				     // on quads and later restore
+ 				     // them
+     std::vector<bool> user_flags;
+     this->get_tria().save_user_flags_quad(user_flags);
+     const_cast<Triangulation<dim,spacedim> &>(this->get_tria()).clear_user_flags_quad ();
+ 
                                       // An implementation of the
                                       // algorithm described in the hp
                                       // paper, including the
@@@ -3396,59 -2442,69 +2443,69 @@@
                                       // particularly frequent case
      Table<2,std_cxx1x::shared_ptr<internal::hp::DoFIdentities> >
        quad_dof_identities (finite_elements->size(),
 -                           finite_elements->size());
 +			               finite_elements->size());
  
-     for (quad_iterator quad=begin_quad(); quad!=end_quad(); ++quad)
-       if (quad->n_active_fe_indices() == 2)
- 	    {
-                                            // find out which is the
-                                            // most dominating finite
-                                            // element of the ones that
-                                            // are used on this quad
-           const unsigned int most_dominating_fe_index
-             = internal::hp::get_most_dominating_fe_index<dim,spacedim> (quad);
- 
- 	      const unsigned int n_active_fe_indices
- 	        = quad->n_active_fe_indices ();
- 
-                                            // loop over the indices of
-                                            // all the finite elements
-                                            // that are not dominating,
-                                            // and identify their dofs
-                                            // to the most dominating
-                                            // one
-           for (unsigned int f=0; f<n_active_fe_indices; ++f)
-             if (quad->nth_active_fe_index (f) !=
-                 most_dominating_fe_index)
-               {
-                 const unsigned int
-                   other_fe_index = quad->nth_active_fe_index (f);
- 
-                 internal::hp::ensure_existence_of_dof_identities<2>
-                   ((*finite_elements)[most_dominating_fe_index],
-                    (*finite_elements)[other_fe_index],
-                    quad_dof_identities[most_dominating_fe_index][other_fe_index]);
- 
-                 internal::hp::DoFIdentities &identities
-                   = *quad_dof_identities[most_dominating_fe_index][other_fe_index];
-                 for (unsigned int i=0; i<identities.size(); ++i)
-                   {
-                     const types::global_dof_index master_dof_index
-                       = quad->dof_index (identities[i].first, most_dominating_fe_index);
-                     const types::global_dof_index slave_dof_index
-                       = quad->dof_index (identities[i].second, other_fe_index);
- 
-                     Assert ((new_dof_indices[master_dof_index] ==
-                              numbers::invalid_unsigned_int)
-                             ||
-                             (new_dof_indices[slave_dof_index] ==
-                              master_dof_index),
-                             ExcInternalError());
- 
-                     new_dof_indices[slave_dof_index] = master_dof_index;
- 		        }
- 	        }
- 	    }
+     for (active_cell_iterator cell=begin_active(); cell!=end(); ++cell)
+       for (unsigned int q=0; q<GeometryInfo<dim>::quads_per_cell; ++q)
+ 	if ((cell->quad(q)->user_flag_set() == false)
+ 	    &&
+ 	    (cell->quad(q)->n_active_fe_indices() == 2))
+ 	  {
+ 	    const quad_iterator quad = cell->quad(q);
+ 	    quad->set_user_flag ();
+ 
+ 					     // find out which is the
+ 					     // most dominating finite
+ 					     // element of the ones that
+ 					     // are used on this quad
+ 	    const unsigned int most_dominating_fe_index
+ 	      = internal::hp::get_most_dominating_fe_index<dim,spacedim> (quad);
+ 
+ 	    const unsigned int n_active_fe_indices
+ 	      = quad->n_active_fe_indices ();
+ 
+ 					     // loop over the indices of
+ 					     // all the finite elements
+ 					     // that are not dominating,
+ 					     // and identify their dofs
+ 					     // to the most dominating
+ 					     // one
+ 	    for (unsigned int f=0; f<n_active_fe_indices; ++f)
+ 	      if (quad->nth_active_fe_index (f) !=
+ 		  most_dominating_fe_index)
+ 		{
+ 		  const unsigned int
+ 		    other_fe_index = quad->nth_active_fe_index (f);
+ 
+ 		  internal::hp::ensure_existence_of_dof_identities<2>
+ 		    ((*finite_elements)[most_dominating_fe_index],
+ 		     (*finite_elements)[other_fe_index],
+ 		     quad_dof_identities[most_dominating_fe_index][other_fe_index]);
+ 
+ 		  internal::hp::DoFIdentities &identities
+ 		    = *quad_dof_identities[most_dominating_fe_index][other_fe_index];
+ 		  for (unsigned int i=0; i<identities.size(); ++i)
+ 		    {
 -		      const unsigned int master_dof_index
++		      const types::global_dof_index master_dof_index
+ 			= quad->dof_index (identities[i].first, most_dominating_fe_index);
 -		      const unsigned int slave_dof_index
++		      const types::global_dof_index slave_dof_index
+ 			= quad->dof_index (identities[i].second, other_fe_index);
+ 
+ 		      Assert ((new_dof_indices[master_dof_index] ==
+ 			       numbers::invalid_unsigned_int)
+ 			      ||
+ 			      (new_dof_indices[slave_dof_index] ==
+ 			       master_dof_index),
+ 			      ExcInternalError());
+ 
+ 		      new_dof_indices[slave_dof_index] = master_dof_index;
+ 		    }
+ 		}
+ 	  }
+ 
+ 				     // finally restore the user flags
+     const_cast<Triangulation<dim,spacedim> &>(this->get_tria())
+       .load_user_flags_quad(user_flags);
    }
  
  
@@@ -3739,9 -2812,12 +2814,12 @@@
    template<>
    void
    DoFHandler<2,2>::
 -  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
 -                          internal::int2type<2>)
 +  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
 +			              internal::int2type<2>)
    {
+     const unsigned int dim = 2;
+     const unsigned int spacedim = 2;
+ 
      Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
  
      renumber_dofs_internal (new_numbers, internal::int2type<1>());
@@@ -3769,9 -2864,12 +2866,12 @@@
    template<>
    void
    DoFHandler<2,3>::
 -  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
 -                          internal::int2type<2>)
 +  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
 +			              internal::int2type<2>)
    {
+     const unsigned int dim = 2;
+     const unsigned int spacedim = 3;
+ 
      Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
  
      renumber_dofs_internal (new_numbers, internal::int2type<1>());
@@@ -3798,9 -2915,12 +2917,12 @@@
    template<>
    void
    DoFHandler<3,3>::
 -  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
 -                          internal::int2type<2>)
 +  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
 +			              internal::int2type<2>)
    {
+     const unsigned int dim = 3;
+     const unsigned int spacedim = 3;
+ 
      Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
  
      renumber_dofs_internal (new_numbers, internal::int2type<1>());
@@@ -3827,9 -2966,12 +2968,12 @@@
    template<>
    void
    DoFHandler<3,3>::
 -  renumber_dofs_internal (const std::vector<unsigned int> &new_numbers,
 -                          internal::int2type<3>)
 +  renumber_dofs_internal (const std::vector<types::global_dof_index> &new_numbers,
 +			              internal::int2type<3>)
    {
+     const unsigned int dim = 3;
+     const unsigned int spacedim = 3;
+ 
      Assert (new_numbers.size() == n_dofs(), ExcRenumberingIncomplete());
  
      renumber_dofs_internal (new_numbers, internal::int2type<2>());