]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Merge from upstream.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 4 Nov 2012 23:40:43 +0000 (23:40 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Sun, 4 Nov 2012 23:40:43 +0000 (23:40 +0000)
git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27372 0785d39b-7218-0410-832d-ea1e28bc413d

24 files changed:
1  2 
deal.II/examples/step-18/step-18.cc
deal.II/examples/step-42/step-42.cc
deal.II/examples/step-6/step-6.cc
deal.II/include/deal.II/dofs/dof_accessor.h
deal.II/include/deal.II/dofs/dof_accessor.templates.h
deal.II/include/deal.II/dofs/dof_handler.h
deal.II/include/deal.II/dofs/dof_tools.h
deal.II/include/deal.II/dofs/number_cache.h
deal.II/include/deal.II/fe/fe_values.h
deal.II/include/deal.II/hp/dof_handler.h
deal.II/include/deal.II/lac/constraint_matrix.h
deal.II/include/deal.II/multigrid/mg_dof_handler.h
deal.II/include/deal.II/numerics/vector_tools.templates.h
deal.II/source/base/index_set.cc
deal.II/source/base/utilities.cc
deal.II/source/dofs/dof_handler.cc
deal.II/source/dofs/dof_renumbering.cc
deal.II/source/dofs/dof_tools.cc
deal.II/source/dofs/dof_tools.inst.in
deal.II/source/fe/fe_tools.cc
deal.II/source/fe/fe_values.cc
deal.II/source/grid/grid_tools.cc
deal.II/source/hp/dof_handler.cc
deal.II/source/numerics/matrix_tools.cc

Simple merge
Simple merge
Simple merge
index ae4b831bc8064b5e724a0ab3c64ba1bf6ae12400,cefac603cc35b427951e6209697664cd3866093f..8643bd7694a2b063202ec638297bb0b2214783a1
@@@ -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
                                               // 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
index 8fabf6a013e541d6d5d4967519b47cc9b587e331,a79738f20e351ab81372272b1df5c9b5a941e334..af4f3ad83dd03b50f4f512456843450b83b272bd
@@@ -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);
  
                                       /**
index c7c30bc3647e97fdc3b619d586a136453f22da2e,ae40ad9d46a8d3cca0c0aa2a8f6f95e14a93aac9..10882b2e97c069e966f0099c42e862004c9d3919
@@@ -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.;
+           }
      }
    }
  
    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);
    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);
    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);
    }
Simple merge
Simple merge
index 0a68adaa5cca92b7532404823b687b584be7020e,afa53fe64e312a959f22b44ceabf0c19005c6662..10ea886f826bc8bf0cb05d4940b68a2b63ae8395
@@@ -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"
index e969f3aabd30cefe8cbadaebcfecfbd6e994927b,05f978abeb4e9258040a88c6f283e69698d7b7a6..f3edb6c3abd7918bf33110c8592b23f09f116ee3
@@@ -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;
                             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);
                             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);
Simple merge
index 9d6f34d4f8908bc4a10310c0ac8ca2dd1060972c,b0f889f3f6e41730c68491d420de0b17e7031c26..ba94bf15cdad9535bced0870f5de6be8bfc4f5d1
@@@ -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>&);
  
 -DoFTools::extract_level_dofs<1,3>
+ template
+ void
++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 &,
Simple merge
Simple merge
index d67815fbe81e6946ed603fca7b025c300c588358,321023958a63007f63adcb6c486e8d76e1b00f60..d5a2255cc8a435127f5b93c366db89cb7361662e
@@@ -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>
      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));
                                           // 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));
                                     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;
                                 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,
      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() &&
  
      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;
                                 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);
    }
  
  
index 41abf8114ccd344149ad09f5984f9049cbdfe456,1e162dcfb941618aebc1f5d8ce2ca5f42b61d9bb..8847f900c129dbd3ef72708fc2375a51c88058f7
@@@ -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
                                       // 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);
    }
  
  
    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
                                       // 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);
    }
  
  
    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>());
    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>());
    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>());
    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>());

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.