From: kronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Fri, 4 Jan 2013 10:10:48 +0000 (+0000)
Subject: Delete MGDoFHandler from matrix free class. These features are now present in DoFHandler.
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=513c5e724279ed3ae29a75017d38ac509a3e380b;p=dealii-svn.git

Delete MGDoFHandler from matrix free class. These features are now present in DoFHandler.

git-svn-id: https://svn.dealii.org/trunk@27919 0785d39b-7218-0410-832d-ea1e28bc413d
---

diff --git a/deal.II/examples/step-37/step-37.cc b/deal.II/examples/step-37/step-37.cc
index 71655eae6f..fb8834e003 100644
--- a/deal.II/examples/step-37/step-37.cc
+++ b/deal.II/examples/step-37/step-37.cc
@@ -31,7 +31,6 @@
 #include <deal.II/grid/grid_generator.h>
 
 #include <deal.II/multigrid/multigrid.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
 #include <deal.II/multigrid/mg_dof_accessor.h>
 #include <deal.II/multigrid/mg_transfer.h>
 #include <deal.II/multigrid/mg_tools.h>
@@ -273,9 +272,9 @@ namespace Step37
 
     void clear();
 
-    void reinit (const MGDoFHandler<dim> &dof_handler,
+    void reinit (const DoFHandler<dim>  &dof_handler,
                  const ConstraintMatrix &constraints,
-                 const unsigned int       level = numbers::invalid_unsigned_int);
+                 const unsigned int      level = numbers::invalid_unsigned_int);
 
     unsigned int m () const;
     unsigned int n () const;
@@ -398,7 +397,7 @@ namespace Step37
   // FEEvaluation class (and its template arguments) will be explained below.
   template <int dim, int fe_degree, typename number>
   void
-  LaplaceOperator<dim,fe_degree,number>::reinit (const MGDoFHandler<dim> &dof_handler,
+  LaplaceOperator<dim,fe_degree,number>::reinit (const DoFHandler<dim>  &dof_handler,
                                                  const ConstraintMatrix &constraints,
                                                  const unsigned int      level)
   {
@@ -734,7 +733,7 @@ namespace Step37
 
     Triangulation<dim>               triangulation;
     FE_Q<dim>                        fe;
-    MGDoFHandler<dim>                mg_dof_handler;
+    DoFHandler<dim>                  dof_handler;
     ConstraintMatrix                 constraints;
 
     SystemMatrixType                 system_matrix;
@@ -760,7 +759,7 @@ namespace Step37
   LaplaceProblem<dim>::LaplaceProblem ()
     :
     fe (degree_finite_element),
-    mg_dof_handler (triangulation),
+    dof_handler (triangulation),
     time_details (std::cout, false)
   {}
 
@@ -789,14 +788,15 @@ namespace Step37
     mg_matrices.clear();
     mg_constraints.clear();
 
-    mg_dof_handler.distribute_dofs (fe);
+    dof_handler.distribute_dofs (fe);
+    dof_handler.distribute_mg_dofs (fe);
 
     std::cout << "Number of degrees of freedom: "
-              << mg_dof_handler.n_dofs()
+              << dof_handler.n_dofs()
               << std::endl;
 
     constraints.clear();
-    VectorTools::interpolate_boundary_values (mg_dof_handler,
+    VectorTools::interpolate_boundary_values (dof_handler,
                                               0,
                                               ZeroFunction<dim>(),
                                               constraints);
@@ -806,15 +806,15 @@ namespace Step37
                  << time() << "s/" << time.wall_time() << "s" << std::endl;
     time.restart();
 
-    system_matrix.reinit (mg_dof_handler, constraints);
+    system_matrix.reinit (dof_handler, constraints);
     std::cout.precision(4);
     std::cout << "System matrix memory consumption:     "
               << system_matrix.memory_consumption()*1e-6
               << " MB."
               << std::endl;
 
-    solution.reinit (mg_dof_handler.n_dofs());
-    system_rhs.reinit (mg_dof_handler.n_dofs());
+    solution.reinit (dof_handler.n_dofs());
+    system_rhs.reinit (dof_handler.n_dofs());
 
     setup_time += time.wall_time();
     time_details << "Setup matrix-free system   (CPU/wall) "
@@ -837,7 +837,7 @@ namespace Step37
     ZeroFunction<dim>               homogeneous_dirichlet_bc (1);
     dirichlet_boundary[0] = &homogeneous_dirichlet_bc;
     std::vector<std::set<unsigned int> > boundary_indices(triangulation.n_levels());
-    MGTools::make_boundary_list (mg_dof_handler,
+    MGTools::make_boundary_list (dof_handler,
                                  dirichlet_boundary,
                                  boundary_indices);
     for (unsigned int level=0; level<nlevels; ++level)
@@ -847,12 +847,12 @@ namespace Step37
           mg_constraints[level].add_line(*bc_it);
 
         mg_constraints[level].close();
-        mg_matrices[level].reinit(mg_dof_handler,
+        mg_matrices[level].reinit(dof_handler,
                                   mg_constraints[level],
                                   level);
       }
-    coarse_matrix.reinit (mg_dof_handler.n_dofs(0),
-                          mg_dof_handler.n_dofs(0));
+    coarse_matrix.reinit (dof_handler.n_dofs(0),
+                          dof_handler.n_dofs(0));
     setup_time += time.wall_time();
     time_details << "Setup matrix-free levels   (CPU/wall) "
                  << time() << "s/" << time.wall_time() << "s" << std::endl;
@@ -881,8 +881,8 @@ namespace Step37
     const Coefficient<dim> coefficient;
     std::vector<double>    coefficient_values (n_q_points);
 
-    typename DoFHandler<dim>::active_cell_iterator cell = mg_dof_handler.begin_active(),
-                                                   endc = mg_dof_handler.end();
+    typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+                                                   endc = dof_handler.end();
     for (; cell!=endc; ++cell)
       {
         cell->get_dof_indices (local_dof_indices);
@@ -932,11 +932,11 @@ namespace Step37
     const unsigned int n_levels = triangulation.n_levels();
     std::vector<Vector<float> > diagonals (n_levels);
     for (unsigned int level=0; level<n_levels; ++level)
-      diagonals[level].reinit (mg_dof_handler.n_dofs(level));
+      diagonals[level].reinit (dof_handler.n_dofs(level));
 
     std::vector<unsigned int> cell_no(triangulation.n_levels());
-    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                              endc = mg_dof_handler.end();
+    typename DoFHandler<dim>::cell_iterator cell = dof_handler.begin(),
+                                            endc = dof_handler.end();
     for (; cell!=endc; ++cell)
       {
         const unsigned int level = cell->level();
@@ -1004,7 +1004,7 @@ namespace Step37
     GrowingVectorMemory<>   vector_memory;
 
     MGTransferPrebuilt<Vector<double> > mg_transfer;
-    mg_transfer.build_matrices(mg_dof_handler);
+    mg_transfer.build_matrices(dof_handler);
     setup_time += time.wall_time();
     time_details << "MG build transfer time     (CPU/wall) " << time()
                  << "s/" << time.wall_time() << "s\n";
@@ -1041,7 +1041,7 @@ namespace Step37
     MGMatrix<LevelMatrixType, Vector<double> >
     mg_matrix(&mg_matrices);
 
-    Multigrid<Vector<double> > mg(mg_dof_handler,
+    Multigrid<Vector<double> > mg(dof_handler,
                                   mg_matrix,
                                   mg_coarse,
                                   mg_transfer,
@@ -1049,7 +1049,7 @@ namespace Step37
                                   mg_smoother);
     PreconditionMG<dim, Vector<double>,
                    MGTransferPrebuilt<Vector<double> > >
-                   preconditioner(mg_dof_handler, mg, mg_transfer);
+                   preconditioner(dof_handler, mg, mg_transfer);
 
     // Finally, write out the memory consumption of the Multigrid object (or
     // rather, of its most significant components, since there is no built-in
@@ -1102,7 +1102,7 @@ namespace Step37
   {
     DataOut<dim> data_out;
 
-    data_out.attach_dof_handler (mg_dof_handler);
+    data_out.attach_dof_handler (dof_handler);
     data_out.add_data_vector (solution, "solution");
     data_out.build_patches ();
 
diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.h b/deal.II/include/deal.II/matrix_free/matrix_free.h
index 6a9ad8cafd..525a33b3df 100644
--- a/deal.II/include/deal.II/matrix_free/matrix_free.h
+++ b/deal.II/include/deal.II/matrix_free/matrix_free.h
@@ -224,10 +224,12 @@ public:
     UpdateFlags         mapping_update_flags;
 
     /**
-     * If working on a MGDoFHandler, this option can be used to define whether
-     * we work on a certain level of the mesh, and not the active cells. If
-     * set to invalid_unsigned_int (which is the default value), the active
-     * cells are gone through, otherwise the level given by this parameter.
+     * This option can be used to define whether we work on a certain level of
+     * the mesh, and not the active cells. If set to invalid_unsigned_int
+     * (which is the default value), the active cells are gone through,
+     * otherwise the level given by this parameter. Note that if you specify
+     * to work on a level, its dofs must be distributed by using
+     * <code>dof_handler.distribute_mg_dofs(fe);</code>.
      */
     unsigned int        level_mg_handler;
 
@@ -558,7 +560,7 @@ public:
   get_constrained_dofs (const unsigned int fe_component = 0) const;
 
   /**
-   * Calls renumber_dofs function in dof info which renumbers the the degrees
+   * Calls renumber_dofs function in dof info which renumbers the degrees
    * of freedom according to the ordering for parallelization.
    */
   void renumber_dofs (std::vector<unsigned int> &renumbering,
@@ -599,27 +601,14 @@ public:
 
   /**
    * In case this structure was built based on a DoFHandler, this returns the
-   * DoFHandler. Note that this function returns an exception in case the
-   * structure was based on MGDoFHandler and a level has been specified in
-   * InitializationOption.
+   * DoFHandler.
    */
   const DoFHandler<dim> &
   get_dof_handler (const unsigned int fe_component = 0) const;
 
-  /**
-   * In case this structure was built based on a DoFHandler, this returns the
-   * DoFHandler. Note that this function returns an exception in case the
-   * structure was based on MGDoFHandler and a level has been specified in
-   * InitializationOption.
-   */
-  const MGDoFHandler<dim> &
-  get_mg_dof_handler (const unsigned int fe_component = 0) const;
-
   /**
    * This returns the cell iterator in deal.II speak to a given cell in the
-   * renumbering of this structure. This function returns an exception in case
-   * the structure was constructed based on an MGDoFHandler with level
-   * specified, as these cells are in general not active.
+   * renumbering of this structure.
    *
    * Note that the cell iterators in deal.II go through cells differently to
    * what the cell loop of this class does. This is because several cells are
@@ -632,22 +621,6 @@ public:
                      const unsigned int vector_number,
                      const unsigned int fe_component = 0) const;
 
-  /**
-   * This returns the cell iterator in deal.II speak to a given cell in the
-   * renumbering of this structure. This function returns an exception in case
-   * the structure was constructed based on a DoFHandler and not MGDoFHandler.
-   *
-   * Note that the cell iterators in deal.II go through cells differently to
-   * what the cell loop of this class does. This is because several cells are
-   * worked on together (vectorization), and since cells with neighbors on
-   * different MPI processors need to be accessed at a certain time when
-   * accessing remote data and overlapping communication with computation.
-   */
-  typename MGDoFHandler<dim>::cell_iterator
-  get_mg_cell_iterator (const unsigned int macro_cell_number,
-                        const unsigned int vector_number,
-                        const unsigned int fe_component = 0) const;
-
   /**
    * This returns the cell iterator in deal.II speak to a given cell in the
    * renumbering of this structure. This function returns an exception in case
@@ -829,11 +802,10 @@ private:
 
   /**
    * This is the actual reinit function that sets up the indices for the
-   * DoFHandler and MGDoFHandler case.
+   * DoFHandler case.
    */
-  template <typename DoFHandler>
   void internal_reinit (const Mapping<dim>                &mapping,
-                        const std::vector<const DoFHandler *> &dof_handler,
+                        const std::vector<const DoFHandler<dim> *> &dof_handler,
                         const std::vector<const ConstraintMatrix *> &constraint,
                         const std::vector<IndexSet>       &locally_owned_set,
                         const std::vector<hp::QCollection<1> > &quad,
@@ -866,13 +838,7 @@ private:
                                 const unsigned int                         level);
 
   /**
-   * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
-   */
-  void initialize_dof_handlers (const std::vector<const MGDoFHandler<dim>*> &dof_handlers,
-                                const unsigned int                           level);
-
-  /**
-   * Initializes the DoFHandlers based on a DoFHandler<dim> argument.
+   * Initializes the DoFHandlers based on a hp::DoFHandler<dim> argument.
    */
   void initialize_dof_handlers (const std::vector<const hp::DoFHandler<dim>*> &dof_handlers,
                                 const unsigned int                             level);
@@ -886,9 +852,8 @@ private:
   {
     DoFHandlers () : n_dof_handlers (0), level (numbers::invalid_unsigned_int) {};
     std::vector<SmartPointer<const DoFHandler<dim> > >   dof_handler;
-    std::vector<SmartPointer<const MGDoFHandler<dim> > > mg_dof_handler;
     std::vector<SmartPointer<const hp::DoFHandler<dim> > > hp_dof_handler;
-    enum ActiveDoFHandler { usual, multigrid, hp } active_dof_handler;
+    enum ActiveDoFHandler { usual, hp } active_dof_handler;
     unsigned int n_dof_handlers;
     unsigned int level;
   };
@@ -1213,15 +1178,6 @@ MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
                        dof_handlers.n_dof_handlers);
       return *dof_handlers.dof_handler[dof_index];
     }
-  else if (dof_handlers.active_dof_handler == DoFHandlers::multigrid)
-    {
-      Assert (dof_handlers.level == numbers::invalid_unsigned_int,
-              ExcMessage("Cannot return DoFHandler when MGDoFHandler and "
-                         "level are specified."));
-      AssertDimension (dof_handlers.mg_dof_handler.size(),
-                       dof_handlers.n_dof_handlers);
-      return *dof_handlers.mg_dof_handler[dof_index];
-    }
   else
     {
       Assert (false, ExcNotImplemented());
@@ -1233,19 +1189,6 @@ MatrixFree<dim,Number>::get_dof_handler (const unsigned int dof_index) const
 
 
 
-template <int dim, typename Number>
-inline
-const MGDoFHandler<dim> &
-MatrixFree<dim,Number>::get_mg_dof_handler (const unsigned int dof_index) const
-{
-  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
-  Assert (dof_handlers.active_dof_handler == DoFHandlers::multigrid,
-          ExcNotImplemented());
-  return *dof_handlers.mg_dof_handler[dof_index];
-}
-
-
-
 template <int dim, typename Number>
 inline
 typename DoFHandler<dim>::active_cell_iterator
@@ -1271,15 +1214,6 @@ MatrixFree<dim,Number>::get_cell_iterator(const unsigned int macro_cell_number,
                        dof_handlers.n_dof_handlers);
       dofh = dof_handlers.dof_handler[dof_index];
     }
-  else if (dof_handlers.active_dof_handler == DoFHandlers::multigrid)
-    {
-      Assert (dof_handlers.level == numbers::invalid_unsigned_int,
-              ExcMessage("Cannot return DoFHandler when MGDoFHandler and "
-                         "level are specified."));
-      AssertDimension (dof_handlers.mg_dof_handler.size(),
-                       dof_handlers.n_dof_handlers);
-      dofh = dof_handlers.mg_dof_handler[dof_index];
-    }
   else
     {
       Assert (false, ExcMessage ("Cannot return DoFHandler<dim>::cell_iterator "
@@ -1294,36 +1228,6 @@ MatrixFree<dim,Number>::get_cell_iterator(const unsigned int macro_cell_number,
 
 
 
-template <int dim, typename Number>
-inline
-typename MGDoFHandler<dim>::cell_iterator
-MatrixFree<dim,Number>::get_mg_cell_iterator(const unsigned int macro_cell_number,
-                                             const unsigned int vector_number,
-                                             const unsigned int dof_index) const
-{
-  const unsigned int vectorization_length=VectorizedArray<Number>::n_array_elements;
-#ifdef DEBUG
-  AssertIndexRange (dof_index, dof_handlers.n_dof_handlers);
-  AssertIndexRange (macro_cell_number, size_info.n_macro_cells);
-  AssertIndexRange (vector_number, vectorization_length);
-  const unsigned int irreg_filled =
-    std_cxx1x::get<2>(dof_info[dof_index].row_starts[macro_cell_number]);
-  if (irreg_filled > 0)
-    AssertIndexRange (vector_number, irreg_filled);
-#endif
-
-  Assert (dof_handlers.active_dof_handler == DoFHandlers::multigrid,
-          ExcNotImplemented());
-  const MGDoFHandler<dim> *dofh = dof_handlers.mg_dof_handler[dof_index];
-
-  std::pair<unsigned int,unsigned int> index =
-    cell_level_index[macro_cell_number*vectorization_length+vector_number];
-  return typename DoFHandler<dim>::cell_iterator
-         (&dofh->get_tria(), index.first, index.second, dofh);
-}
-
-
-
 template <int dim, typename Number>
 inline
 typename hp::DoFHandler<dim>::active_cell_iterator
@@ -1520,37 +1424,28 @@ MatrixFree<dim,Number>::mapping_initialized () const
 
 namespace internal
 {
-  template <typename DH>
-  std::vector<IndexSet>
-  extract_locally_owned_index_sets (const std::vector<const DH *> &dofh,
-                                    const unsigned int)
-  {
-    std::vector<IndexSet> locally_owned_set;
-    locally_owned_set.reserve (dofh.size());
-    for (unsigned int j=0; j<dofh.size(); j++)
-      locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
-    return locally_owned_set;
-  }
-
-
-
-  template <int dim>
-  std::vector<IndexSet>
-  extract_locally_owned_index_sets (const std::vector<const dealii::MGDoFHandler<dim>*> &dofh,
-                                    const unsigned int level)
+  namespace MatrixFree
   {
-    std::vector<IndexSet> locally_owned_set;
-    locally_owned_set.reserve (dofh.size());
-    for (unsigned int j=0; j<dofh.size(); j++)
-      if (level == numbers::invalid_unsigned_int)
-        locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
-      else
-        {
-          IndexSet new_set (dofh[j]->n_dofs(level));
-          new_set.add_range (0, dofh[j]->n_dofs(level));
-          locally_owned_set.push_back(new_set);
-        }
-    return locally_owned_set;
+    template <typename DH>
+    inline
+    std::vector<IndexSet>
+    extract_locally_owned_index_sets (const std::vector<const DH *> &dofh,
+                                      const unsigned int level)
+    {
+      std::vector<IndexSet> locally_owned_set;
+      locally_owned_set.reserve (dofh.size());
+      for (unsigned int j=0; j<dofh.size(); j++)
+        if (level == numbers::invalid_unsigned_int)
+          locally_owned_set.push_back(dofh[j]->locally_owned_dofs());
+        else
+          {
+            // TODO: not distributed yet
+            IndexSet new_set (dofh[j]->n_dofs(level));
+            new_set.add_range (0, dofh[j]->n_dofs(level));
+            locally_owned_set.push_back(new_set);
+          }
+      return locally_owned_set;
+    }
   }
 }
 
@@ -1574,8 +1469,8 @@ reinit(const DH               &dof_handler,
   quads.push_back (quad);
 
   std::vector<IndexSet> locally_owned_sets =
-    internal::extract_locally_owned_index_sets (dof_handlers,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handlers, additional_data.level_mg_handler);
   reinit(mapping, dof_handlers,constraints, locally_owned_sets, quads,
          additional_data);
 }
@@ -1600,8 +1495,8 @@ reinit(const Mapping<dim>     &mapping,
   quads.push_back (quad);
 
   std::vector<IndexSet> locally_owned_sets =
-    internal::extract_locally_owned_index_sets (dof_handlers,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handlers, additional_data.level_mg_handler);
   reinit(mapping, dof_handlers,constraints,locally_owned_sets, quads,
          additional_data);
 }
@@ -1618,8 +1513,8 @@ reinit(const std::vector<const DH *>               &dof_handler,
 {
   MappingQ1<dim> mapping;
   std::vector<IndexSet> locally_owned_set =
-    internal::extract_locally_owned_index_sets (dof_handler,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handler, additional_data.level_mg_handler);
   reinit(mapping, dof_handler,constraint,locally_owned_set,
          static_cast<const std::vector<Quadrature<1> >&>(quad),
          additional_data);
@@ -1639,8 +1534,8 @@ reinit(const std::vector<const DH *>               &dof_handler,
   std::vector<Quad> quads;
   quads.push_back(quad);
   std::vector<IndexSet> locally_owned_set =
-    internal::extract_locally_owned_index_sets (dof_handler,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handler, additional_data.level_mg_handler);
   reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
          additional_data);
 }
@@ -1659,8 +1554,8 @@ reinit(const Mapping<dim>                         &mapping,
   std::vector<Quad> quads;
   quads.push_back(quad);
   std::vector<IndexSet> locally_owned_set =
-    internal::extract_locally_owned_index_sets (dof_handler,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handler, additional_data.level_mg_handler);
   reinit(mapping, dof_handler,constraint,locally_owned_set, quads,
          additional_data);
 }
@@ -1677,14 +1572,46 @@ reinit(const Mapping<dim>                         &mapping,
        const MatrixFree<dim,Number>::AdditionalData additional_data)
 {
   std::vector<IndexSet> locally_owned_set =
-    internal::extract_locally_owned_index_sets (dof_handler,
-                                                additional_data.level_mg_handler);
+    internal::MatrixFree::extract_locally_owned_index_sets
+    (dof_handler, additional_data.level_mg_handler);
   reinit(mapping, dof_handler,constraint,locally_owned_set,
          quad, additional_data);
 }
 
 
 
+namespace internal
+{
+  namespace MatrixFree
+  {
+    // resolve DoFHandler types
+
+    // MGDoFHandler is deprecated in deal.II but might still be present in
+    // user code, so we need to resolve its type (fortunately, it is derived
+    // from DoFHandler, so we can static_cast it to a DoFHandler<dim>)
+    template <typename DH>
+    inline
+    std::vector<const dealii::DoFHandler<DH::dimension> *>
+    resolve_dof_handler (const std::vector<const DH *> &dof_handler)
+    {
+      std::vector<const dealii::DoFHandler<DH::dimension> *> conversion(dof_handler.size());
+      for (unsigned int i=0; i<dof_handler.size(); ++i)
+        conversion[i] = static_cast<const dealii::DoFHandler<DH::dimension> *>(dof_handler[i]);
+      return conversion;
+    }
+
+    template <int dim>
+    inline
+    std::vector<const dealii::hp::DoFHandler<dim> *>
+    resolve_dof_handler (const std::vector<const dealii::hp::DoFHandler<dim> *> &dof_handler)
+    {
+      return dof_handler;
+    }
+  }
+}
+
+
+
 template <int dim, typename Number>
 template <typename DH, typename Quad>
 void MatrixFree<dim,Number>::
@@ -1699,8 +1626,9 @@ reinit(const Mapping<dim>                         &mapping,
   std::vector<hp::QCollection<1> > quad_hp;
   for (unsigned int q=0; q<quad.size(); ++q)
     quad_hp.push_back (hp::QCollection<1>(quad[q]));
-  internal_reinit (mapping, dof_handler, constraint, locally_owned_set,
-                   quad_hp, additional_data);
+  internal_reinit (mapping,
+                   internal::MatrixFree::resolve_dof_handler(dof_handler),
+                   constraint, locally_owned_set, quad_hp, additional_data);
 }
 
 
diff --git a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h
index c5a7efc502..385d0f546b 100644
--- a/deal.II/include/deal.II/matrix_free/matrix_free.templates.h
+++ b/deal.II/include/deal.II/matrix_free/matrix_free.templates.h
@@ -99,13 +99,12 @@ namespace internal
 
 
 template <int dim, typename Number>
-template <typename DH>
 void MatrixFree<dim,Number>::
-internal_reinit(const Mapping<dim>                         &mapping,
-                const std::vector<const DH *>               &dof_handler,
+internal_reinit(const Mapping<dim>                          &mapping,
+                const std::vector<const DoFHandler<dim> *>  &dof_handler,
                 const std::vector<const ConstraintMatrix *> &constraint,
-                const std::vector<IndexSet>                &locally_owned_set,
-                const std::vector<hp::QCollection<1> >     &quad,
+                const std::vector<IndexSet>                 &locally_owned_set,
+                const std::vector<hp::QCollection<1> >      &quad,
                 const MatrixFree<dim,Number>::AdditionalData additional_data)
 {
   if (additional_data.initialize_indices == true)
@@ -323,9 +322,10 @@ namespace internal
 template <int dim, typename Number>
 void MatrixFree<dim,Number>::
 initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handler,
-                         const unsigned int)
+                         const unsigned int level)
 {
   dof_handlers.active_dof_handler = DoFHandlers::usual;
+  dof_handlers.level = level;
   dof_handlers.n_dof_handlers = dof_handler.size();
   dof_handlers.dof_handler.resize (dof_handlers.n_dof_handlers);
   for (unsigned int no=0; no<dof_handlers.n_dof_handlers; ++no)
@@ -340,41 +340,6 @@ initialize_dof_handlers (const std::vector<const DoFHandler<dim>*> &dof_handler,
   const unsigned int my_pid = size_info.my_pid;
 
   const Triangulation<dim> &tria = dof_handlers.dof_handler[0]->get_tria();
-  {
-    if (n_mpi_procs == 1)
-      cell_level_index.reserve (tria.n_active_cells());
-    typename Triangulation<dim>::cell_iterator cell = tria.begin(0),
-                                               end_cell = tria.end(0);
-    for ( ; cell != end_cell; ++cell)
-      internal::resolve_cell (cell, cell_level_index, my_pid);
-  }
-}
-
-
-template <int dim, typename Number>
-void MatrixFree<dim,Number>::
-initialize_dof_handlers (const std::vector<const MGDoFHandler<dim>*> &dof_handler,
-                         const unsigned int                           level)
-{
-  dof_handlers.active_dof_handler = DoFHandlers::multigrid;
-  dof_handlers.level = level;
-  dof_handlers.n_dof_handlers = dof_handler.size();
-  dof_handlers.mg_dof_handler.resize (dof_handlers.n_dof_handlers);
-  for (unsigned int no=0; no<dof_handlers.n_dof_handlers; ++no)
-    dof_handlers.mg_dof_handler[no] = dof_handler[no];
-
-  dof_info.resize (dof_handlers.n_dof_handlers);
-
-  // go through cells on zeroth level and then successively step down into
-  // children. This gives a z-ordering of the cells, which is beneficial when
-  // setting up neighboring relations between cells for thread parallelization
-  const unsigned int n_mpi_procs = size_info.n_procs;
-  const unsigned int my_pid = size_info.my_pid;
-
-  // if we have no level given, use the same as for the standard DoFHandler,
-  // otherwise we must loop through the respective level
-  const Triangulation<dim> &tria = dof_handlers.mg_dof_handler[0]->get_tria();
-
   if (level == numbers::invalid_unsigned_int)
     {
       if (n_mpi_procs == 1)
@@ -476,9 +441,7 @@ void MatrixFree<dim,Number>::initialize_indices
         }
       else
         {
-          const DoFHandler<dim> *dofh =
-            dof_handlers.active_dof_handler == DoFHandlers::usual ?
-            &*dof_handlers.dof_handler[no] : &*dof_handlers.mg_dof_handler[no];
+          const DoFHandler<dim> *dofh =&*dof_handlers.dof_handler[no];
           fes.push_back (&dofh->get_fe());
           dof_info[no].max_fe_index = 1;
           dof_info[no].fe_index_conversion.resize (1);
@@ -577,15 +540,11 @@ void MatrixFree<dim,Number>::initialize_indices
       bool cell_at_boundary = false;
       for (unsigned int no=0; no<n_fe; ++no)
         {
-          // OK, read indices from standard DoFHandler or active indices in
-          // MGDoFHandler. That is the usual stuff
-          if (dof_handlers.active_dof_handler == DoFHandlers::usual ||
-              (dof_handlers.active_dof_handler == DoFHandlers::multigrid &&
-               dof_handlers.level == numbers::invalid_unsigned_int))
+          // OK, read indices from standard DoFHandler in the usual way
+          if (dof_handlers.active_dof_handler == DoFHandlers::usual &&
+              dof_handlers.level == numbers::invalid_unsigned_int)
             {
-              const DoFHandler<dim> *dofh =
-                dof_handlers.active_dof_handler == DoFHandlers::usual ?
-                &*dof_handlers.dof_handler[no] : &*dof_handlers.mg_dof_handler[no];
+              const DoFHandler<dim> *dofh = &*dof_handlers.dof_handler[no];
               typename DoFHandler<dim>::active_cell_iterator
               cell_it (&dofh->get_tria(),
                        cell_level_index[counter].first,
@@ -600,13 +559,12 @@ void MatrixFree<dim,Number>::initialize_indices
                                              cell_at_boundary);
             }
           // ok, now we are requested to use a level in a MGDoFHandler
-          else if (dof_handlers.active_dof_handler == DoFHandlers::multigrid &&
+          else if (dof_handlers.active_dof_handler == DoFHandlers::usual &&
                    dof_handlers.level != numbers::invalid_unsigned_int)
             {
-              const MGDoFHandler<dim> *dofh =
-                dof_handlers.mg_dof_handler[no];
+              const DoFHandler<dim> *dofh = dof_handlers.dof_handler[no];
               AssertIndexRange (dof_handlers.level, dofh->get_tria().n_levels());
-              typename MGDoFHandler<dim>::cell_iterator
+              typename DoFHandler<dim>::cell_iterator
               cell_it (&dofh->get_tria(),
                        cell_level_index[counter].first,
                        cell_level_index[counter].second,
diff --git a/deal.II/source/numerics/matrix_free.inst.in b/deal.II/source/numerics/matrix_free.inst.in
index a9965f513b..3c2dee08e8 100644
--- a/deal.II/source/numerics/matrix_free.inst.in
+++ b/deal.II/source/numerics/matrix_free.inst.in
@@ -16,39 +16,6 @@ for (deal_II_dimension : DIMENSIONS)
   template class MatrixFree<deal_II_dimension,double>;
   template class MatrixFree<deal_II_dimension,float>;
 
-  // reinit for DoFHandler
-  template void MatrixFree<deal_II_dimension,double>::internal_reinit
-  (const Mapping<deal_II_dimension>                           &,
-   const std::vector<const DoFHandler<deal_II_dimension>*>    &,
-   const std::vector<const ConstraintMatrix*>                 &,
-   const std::vector<IndexSet>                                &,
-   const std::vector<hp::QCollection<1> >                     &,
-   const MatrixFree<deal_II_dimension,double>::AdditionalData);
-  template void MatrixFree<deal_II_dimension,float>::internal_reinit
-  (const Mapping<deal_II_dimension>                           &,
-   const std::vector<const DoFHandler<deal_II_dimension>*>    &,
-   const std::vector<const ConstraintMatrix*>                 &,
-   const std::vector<IndexSet>                                &,
-   const std::vector<hp::QCollection<1> >                     &,
-   const MatrixFree<deal_II_dimension,float>::AdditionalData);
-
-
-  // reinit for MGDoFHandler
-  template void MatrixFree<deal_II_dimension,double>::internal_reinit
-  (const Mapping<deal_II_dimension>                           &,
-   const std::vector<const MGDoFHandler<deal_II_dimension>*>  &,
-   const std::vector<const ConstraintMatrix*>                 &,
-   const std::vector<IndexSet>                                &,
-   const std::vector<hp::QCollection<1> >                     &,
-   const MatrixFree<deal_II_dimension,double>::AdditionalData);
-  template void MatrixFree<deal_II_dimension,float>::internal_reinit
-  (const Mapping<deal_II_dimension>                           &,
-   const std::vector<const MGDoFHandler<deal_II_dimension>*>  &,
-   const std::vector<const ConstraintMatrix*>                 &,
-   const std::vector<IndexSet>                                &,
-   const std::vector<hp::QCollection<1> >                     &,
-   const MatrixFree<deal_II_dimension,float>::AdditionalData);
-
   template void MatrixFree<deal_II_dimension,double>::
   print_memory_consumption<std::ostream> (std::ostream &) const;
   template void MatrixFree<deal_II_dimension,double>::
diff --git a/tests/matrix_free/matrix_vector_mg.cc b/tests/matrix_free/matrix_vector_mg.cc
index 3544ad0728..e4dbf08f17 100644
--- a/tests/matrix_free/matrix_vector_mg.cc
+++ b/tests/matrix_free/matrix_vector_mg.cc
@@ -15,7 +15,6 @@
 std::ofstream logfile("matrix_vector_mg/output");
 
 #include "matrix_vector_common.h"
-#include <deal.II/multigrid/mg_dof_handler.h>
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_matrix.h>
 #include <deal.II/base/function.h>
@@ -34,8 +33,9 @@ void test ()
   FE_Q<dim> fe (fe_degree);
 
                                 // setup DoFs
-  MGDoFHandler<dim> dof(tria);
+  DoFHandler<dim> dof(tria);
   dof.distribute_dofs(fe);
+  dof.distribute_mg_dofs(fe);
   ConstraintMatrix constraints;
   VectorTools::interpolate_boundary_values (dof,
                                             0,
@@ -132,7 +132,7 @@ void test ()
       }
 
     // now to the MG assembly
-    typename MGDoFHandler<dim>::cell_iterator
+    typename DoFHandler<dim>::cell_iterator
       cellm = dof.begin(),
       endcm = dof.end();
     for (; cellm!=endcm; ++cellm)
diff --git a/tests/matrix_free/step-37.cc b/tests/matrix_free/step-37.cc
index 005c227d17..03454d07c4 100644
--- a/tests/matrix_free/step-37.cc
+++ b/tests/matrix_free/step-37.cc
@@ -33,7 +33,6 @@
 #include <deal.II/grid/grid_generator.h>
 
 #include <deal.II/multigrid/multigrid.h>
-#include <deal.II/multigrid/mg_dof_handler.h>
 #include <deal.II/multigrid/mg_dof_accessor.h>
 #include <deal.II/multigrid/mg_transfer.h>
 #include <deal.II/multigrid/mg_tools.h>
@@ -124,7 +123,7 @@ namespace Step37
 
     void clear();
 
-    void reinit (const MGDoFHandler<dim> &dof_handler,
+    void reinit (const DoFHandler<dim>  &dof_handler,
                  const ConstraintMatrix &constraints,
                  const unsigned int       level = numbers::invalid_unsigned_int);
 
@@ -202,7 +201,7 @@ namespace Step37
 
   template <int dim, int fe_degree, typename number>
   void
-  LaplaceOperator<dim,fe_degree,number>::reinit (const MGDoFHandler<dim> &dof_handler,
+  LaplaceOperator<dim,fe_degree,number>::reinit (const DoFHandler<dim>  &dof_handler,
                                                  const ConstraintMatrix &constraints,
                                                  const unsigned int      level)
   {
@@ -376,7 +375,7 @@ namespace Step37
 
     Triangulation<dim>               triangulation;
     FE_Q<dim>                        fe;
-    MGDoFHandler<dim>                mg_dof_handler;
+    DoFHandler<dim>                  mg_dof_handler;
     ConstraintMatrix                 constraints;
 
     SystemMatrixType                 system_matrix;
@@ -408,6 +407,7 @@ namespace Step37
     mg_constraints.clear();
 
     mg_dof_handler.distribute_dofs (fe);
+    mg_dof_handler.distribute_mg_dofs (fe);
 
     deallog << "Number of degrees of freedom: "
             << mg_dof_handler.n_dofs()
@@ -511,8 +511,8 @@ namespace Step37
       diagonals[level].reinit (mg_dof_handler.n_dofs(level));
 
     std::vector<unsigned int> cell_no(triangulation.n_levels());
-    typename MGDoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
-                                              endc = mg_dof_handler.end();
+    typename DoFHandler<dim>::cell_iterator cell = mg_dof_handler.begin(),
+                                            endc = mg_dof_handler.end();
     for (; cell!=endc; ++cell)
       {
         const unsigned int level = cell->level();