From: Karl Ljungkvist <k.ljungkvist@gmail.com>
Date: Thu, 17 Nov 2016 16:26:57 +0000 (+0100)
Subject: Extract some common MG transfer setup code
X-Git-Tag: v8.5.0-rc1~373^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f49d41e35dc0ed2771f8608e48ae75f3b63007cc;p=dealii.git

Extract some common MG transfer setup code
---

diff --git a/include/deal.II/multigrid/mg_transfer_internal.h b/include/deal.II/multigrid/mg_transfer_internal.h
index 27532a3ed6..dec3ea3b11 100644
--- a/include/deal.II/multigrid/mg_transfer_internal.h
+++ b/include/deal.II/multigrid/mg_transfer_internal.h
@@ -17,10 +17,12 @@
 #ifndef dealii__mg_transfer_internal_h
 #define dealii__mg_transfer_internal_h
 
+#include <deal.II/base/mg_level_object.h>
 #include <deal.II/dofs/dof_handler.h>
+#include <deal.II/lac/la_parallel_vector.h>
+#include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/multigrid/mg_constrained_dofs.h>
 
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace internal
@@ -38,6 +40,76 @@ namespace internal
                            std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > &copy_indices,
                            std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > &copy_indices_global_mine,
                            std::vector<std::vector<std::pair<types::global_dof_index, types::global_dof_index> > > &copy_indices_level_mine);
+
+
+
+    /**
+     * Given the collection of child cells in lexicographic ordering as seen
+     * from the parent, this function computes the first index of the given
+     * child
+     */
+    template <int dim>
+    unsigned int
+    compute_shift_within_children(const unsigned int child,
+                                  const unsigned int fe_shift_1d,
+                                  const unsigned int fe_degree);
+
+    /**
+     * Stores data related to the finite element contained in the
+     * DoFHandler. Used only for the initialization using
+     * <tt>setup_transfer</tt>.
+     */
+    template <typename Number>
+    struct ElementInfo
+    {
+
+      /**
+       * Stores the degree of the finite element. The selection of the
+       * computational kernel is based on this number.
+       */
+      unsigned int fe_degree;
+
+      /**
+       * Stores whether the element is continuous and there is a joint degree of
+       * freedom in the center of the 1D line.
+       */
+      bool element_is_continuous;
+
+      /**
+       * Stores the number of components in the finite element.
+       */
+      unsigned int n_components;
+
+      /**
+       * Stores the number of degrees of freedom on all child cells. It is
+       * <tt>2<sup>dim</sup>*fe.dofs_per_cell</tt> for DG elements and somewhat
+       * less for continuous elements.
+       */
+      unsigned int n_child_cell_dofs;
+
+      /**
+       * Holds the one-dimensional embedding (prolongation) matrix from mother
+       * element to the children.
+       */
+      internal::MatrixFreeFunctions::ShapeInfo<Number> shape_info;
+
+    };
+
+    /**
+     * Sets up most of the internal data structures of MGTransferMatrixFree
+     */
+    template <int dim, typename Number>
+    void setup_transfer(const dealii::DoFHandler<dim>                                       &mg_dof,
+                        const MGConstrainedDoFs                                             *mg_constrained_dofs,
+                        ElementInfo<Number>                                                 &elem_info,
+                        std::vector<std::vector<unsigned int> >                             &level_dof_indices,
+                        std::vector<std::vector<std::pair<unsigned int,unsigned int> > >    &parent_child_connect,
+                        std::vector<unsigned int>                                           &n_owned_level_cells,
+                        std::vector<std::vector<std::vector<unsigned short> > >             &dirichlet_indices,
+                        std::vector<AlignedVector<VectorizedArray<Number> > >               &weights_on_refined,
+                        std::vector<std::vector<std::pair<unsigned int, unsigned int> > >   &copy_indices_global_mine,
+                        MGLevelObject<LinearAlgebra::distributed::Vector<Number> >          &ghosted_level_vector);
+
   }
 }
 
diff --git a/source/multigrid/mg_transfer_internal.cc b/source/multigrid/mg_transfer_internal.cc
index 1c3696e078..6561d033a0 100644
--- a/source/multigrid/mg_transfer_internal.cc
+++ b/source/multigrid/mg_transfer_internal.cc
@@ -16,6 +16,7 @@
 
 #include <deal.II/distributed/tria.h>
 #include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_tools.h>
 #include <deal.II/multigrid/mg_transfer_internal.h>
 
 DEAL_II_NAMESPACE_OPEN
@@ -24,11 +25,10 @@ namespace internal
 {
   namespace MGTransfer
   {
-    /**
-     * Internal data structure that is used in the MPI communication in
-     * fill_copy_indices().  It represents an entry in the copy_indices* map,
-     * that associates a level dof index with a global dof index.
-     */
+
+    // Internal data structure that is used in the MPI communication in
+    // fill_copy_indices().  It represents an entry in the copy_indices* map,
+    // that associates a level dof index with a global dof index.
     struct DoFPair
     {
       unsigned int level;
@@ -52,6 +52,8 @@ namespace internal
       {}
     };
 
+    // Internal function for filling the copy indices from global to level
+    // indices
     template <int dim, int spacedim>
     void fill_copy_indices(const dealii::DoFHandler<dim,spacedim>                                                  &mg_dof,
                            const MGConstrainedDoFs                                                                 *mg_constrained_dofs,
@@ -272,6 +274,467 @@ namespace internal
       for (unsigned int level=0; level<copy_indices_global_mine.size(); ++level)
         std::sort(copy_indices_global_mine[level].begin(), copy_indices_global_mine[level].end(), compare);
     }
+
+
+
+    // initialize the vectors needed for the transfer (and merge with the
+    // content in copy_indices_global_mine)
+    template <typename Number>
+    void
+    reinit_ghosted_vector(const IndexSet &locally_owned,
+                          std::vector<types::global_dof_index> &ghosted_level_dofs,
+                          const MPI_Comm &communicator,
+                          LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
+                          std::vector<std::pair<unsigned int,unsigned int> > &copy_indices_global_mine)
+    {
+      std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
+      IndexSet ghosted_dofs(locally_owned.size());
+      ghosted_dofs.add_indices(ghosted_level_dofs.begin(),
+                               std::unique(ghosted_level_dofs.begin(),
+                                           ghosted_level_dofs.end()));
+      ghosted_dofs.compress();
+
+      // Add possible ghosts from the previous content in the vector
+      if (ghosted_level_vector.size() == locally_owned.size())
+        {
+          // shift the local number of the copy indices according to the new
+          // partitioner that we are going to use for the vector
+          const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> part
+            = ghosted_level_vector.get_partitioner();
+          ghosted_dofs.add_indices(part->ghost_indices());
+          for (unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
+            copy_indices_global_mine[i].second =
+              locally_owned.n_elements() +
+              ghosted_dofs.index_within_set(part->local_to_global(copy_indices_global_mine[i].second));
+        }
+      ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator);
+    }
+
+    // Transform the ghost indices to local index space for the vector
+    inline void
+    copy_indices_to_mpi_local_numbers(const Utilities::MPI::Partitioner &part,
+                                      const std::vector<types::global_dof_index> &mine,
+                                      const std::vector<types::global_dof_index> &remote,
+                                      std::vector<unsigned int> &localized_indices)
+    {
+      localized_indices.resize(mine.size()+remote.size(),
+                               numbers::invalid_unsigned_int);
+      for (unsigned int i=0; i<mine.size(); ++i)
+        if (mine[i] != numbers::invalid_dof_index)
+          localized_indices[i] = part.global_to_local(mine[i]);
+
+      for (unsigned int i=0; i<remote.size(); ++i)
+        if (remote[i] != numbers::invalid_dof_index)
+          localized_indices[i+mine.size()] = part.global_to_local(remote[i]);
+    }
+
+    // given the collection of child cells in lexicographic ordering as seen
+    // from the parent, compute the first index of the given child
+    template <int dim>
+    unsigned int
+    compute_shift_within_children(const unsigned int child,
+                                  const unsigned int fe_shift_1d,
+                                  const unsigned int fe_degree)
+    {
+      // we put the degrees of freedom of all child cells in lexicographic
+      // ordering
+      unsigned int c_tensor_index[dim];
+      unsigned int tmp = child;
+      for (unsigned int d=0; d<dim; ++d)
+        {
+          c_tensor_index[d] = tmp % 2;
+          tmp /= 2;
+        }
+      const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
+      unsigned int factor = 1;
+      unsigned int shift = fe_shift_1d * c_tensor_index[0];
+      for (unsigned int d=1; d<dim; ++d)
+        {
+          factor *= n_child_dofs_1d;
+          shift = shift + factor * fe_shift_1d * c_tensor_index[d];
+        }
+      return shift;
+    }
+
+    // puts the indices on the given child cell in lexicographic ordering with
+    // respect to the collection of all child cells as seen from the parent
+    template <int dim>
+    void add_child_indices(const unsigned int child,
+                           const unsigned int fe_shift_1d,
+                           const unsigned int fe_degree,
+                           const std::vector<unsigned int> &lexicographic_numbering,
+                           const std::vector<types::global_dof_index> &local_dof_indices,
+                           types::global_dof_index *target_indices)
+    {
+      const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
+      const unsigned int shift =
+        compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
+      const unsigned int n_components =
+        local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
+      types::global_dof_index *indices = target_indices + shift;
+      const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
+      for (unsigned int c=0, m=0; c<n_components; ++c)
+        for (unsigned int k=0; k<(dim>2 ? (fe_degree+1) : 1); ++k)
+          for (unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j)
+            for (unsigned int i=0; i<(fe_degree+1); ++i, ++m)
+              {
+                const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d*
+                                           n_child_dofs_1d+j*n_child_dofs_1d+i;
+                Assert(indices[index] == numbers::invalid_dof_index ||
+                       indices[index] == local_dof_indices[lexicographic_numbering[m]],
+                       ExcInternalError());
+                indices[index] = local_dof_indices[lexicographic_numbering[m]];
+              }
+    }
+
+    template <int dim, typename Number>
+    void setup_element_info(ElementInfo<Number> &elem_info,const FiniteElement<1> &fe,
+                            const dealii::DoFHandler<dim> &mg_dof)
+    {
+      // currently, we have only FE_Q and FE_DGQ type elements implemented
+      elem_info.n_components = mg_dof.get_fe().element_multiplicity(0);
+      AssertDimension(Utilities::fixed_power<dim>(fe.dofs_per_cell)*elem_info.n_components,
+                      mg_dof.get_fe().dofs_per_cell);
+      AssertDimension(fe.degree, mg_dof.get_fe().degree);
+      elem_info.fe_degree = fe.degree;
+      elem_info.element_is_continuous = fe.dofs_per_vertex > 0;
+      Assert(fe.dofs_per_vertex < 2, ExcNotImplemented());
+
+      // step 1.2: get renumbering of 1D basis functions to lexicographic
+      // numbers. The distinction according to fe.dofs_per_vertex is to support
+      // both continuous and discontinuous bases.
+      std::vector<unsigned int> renumbering(fe.dofs_per_cell);
+      {
+        AssertIndexRange(fe.dofs_per_vertex, 2);
+        renumbering[0] = 0;
+        for (unsigned int i=0; i<fe.dofs_per_line; ++i)
+          renumbering[i+fe.dofs_per_vertex] =
+            GeometryInfo<1>::vertices_per_cell*fe.dofs_per_vertex + i;
+        if (fe.dofs_per_vertex > 0)
+          renumbering[fe.dofs_per_cell-fe.dofs_per_vertex] = fe.dofs_per_vertex;
+      }
+
+      // step 1.3: create a 1D quadrature formula from the finite element that
+      // collects the support points of the basis functions on the two children.
+      std::vector<Point<1> > basic_support_points = fe.get_unit_support_points();
+      Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
+             ExcNotImplemented());
+      std::vector<Point<1> > points_refined(fe.dofs_per_vertex > 0 ?
+                                            (2 * fe.dofs_per_cell - 1) :
+                                            (2 * fe.dofs_per_cell));
+      const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
+      for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
+        for (unsigned int j=0; j<basic_support_points.size(); ++j)
+          points_refined[shift*c+j][0] =
+            c*0.5 + 0.5 * basic_support_points[renumbering[j]][0];
+
+      unsigned int n_child_dofs_1d = points_refined.size();
+      elem_info.n_child_cell_dofs = elem_info.n_components*Utilities::fixed_power<dim>(n_child_dofs_1d);
+
+      // step 1.4: evaluate the polynomials and store the data in ShapeInfo
+      const Quadrature<1> quadrature(points_refined);
+      elem_info.shape_info.reinit(quadrature, mg_dof.get_fe(), 0);
+
+      for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
+        for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
+          for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
+            Assert(std::abs(elem_info.shape_info.shape_values[i*n_child_dofs_1d+j+c*shift][0] -
+                            fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]))
+                   < std::max(2.*(double)std::numeric_limits<Number>::epsilon(),1e-12),
+                   ExcInternalError());
+    }
+
+
+    // Sets up most of the internal data structures of the MGTransferMatrixFree
+    // class
+    template <int dim, typename Number>
+    void setup_transfer(const dealii::DoFHandler<dim>                                       &mg_dof,
+                        const MGConstrainedDoFs                                             *mg_constrained_dofs,
+                        ElementInfo<Number>                                                 &elem_info,
+                        std::vector<std::vector<unsigned int> >                             &level_dof_indices,
+                        std::vector<std::vector<std::pair<unsigned int,unsigned int> > >    &parent_child_connect,
+                        std::vector<unsigned int>                                           &n_owned_level_cells,
+                        std::vector<std::vector<std::vector<unsigned short> > >             &dirichlet_indices,
+                        std::vector<AlignedVector<VectorizedArray<Number> > >               &weights_on_refined,
+                        std::vector<std::vector<std::pair<unsigned int, unsigned int> > >   &copy_indices_global_mine,
+                        MGLevelObject<LinearAlgebra::distributed::Vector<Number> >          &ghosted_level_vector)
+    {
+
+      // we collect all child DoFs of a mother cell together. For faster
+      // tensorized operations, we align the degrees of freedom
+      // lexicographically. We distinguish FE_Q elements and FE_DGQ elements
+
+      const dealii::Triangulation<dim> &tria = mg_dof.get_triangulation();
+
+      // ---------------------------- 1. Extract 1D info about the finite element
+      // step 1.1: create a 1D copy of the finite element from FETools where we
+      // substitute the template argument
+      AssertDimension(mg_dof.get_fe().n_base_elements(), 1);
+      std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
+      {
+        const std::size_t template_starts = fe_name.find_first_of('<');
+        Assert (fe_name[template_starts+1] == (dim==1?'1':(dim==2?'2':'3')),
+                ExcInternalError());
+        fe_name[template_starts+1] = '1';
+      }
+      std_cxx11::shared_ptr<FiniteElement<1> > fe_1d
+      (FETools::get_fe_by_name<1,1>(fe_name));
+      const FiniteElement<1> &fe = *fe_1d;
+
+      setup_element_info(elem_info,fe,mg_dof);
+
+      unsigned int n_child_dofs_1d = elem_info.shape_info.n_q_points;
+
+      // -------------- 2. Extract and match dof indices between child and parent
+      const unsigned int n_levels = tria.n_global_levels();
+      level_dof_indices.resize(n_levels);
+      parent_child_connect.resize(n_levels-1);
+      n_owned_level_cells.resize(n_levels-1);
+      std::vector<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
+      for (unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
+        coarse_level_indices[level].resize(tria.n_raw_cells(level),
+                                           numbers::invalid_unsigned_int);
+      std::vector<types::global_dof_index> local_dof_indices(mg_dof.get_fe().dofs_per_cell);
+      dirichlet_indices.resize(n_levels-1);
+
+      // We use the vectors stored ghosted_level_vector in the base class for
+      // keeping ghosted transfer indices. To avoid keeping two very similar
+      // vectors, we merge them here.
+      if (ghosted_level_vector.max_level() != n_levels-1)
+        ghosted_level_vector.resize(0, n_levels-1);
+
+      for (unsigned int level=n_levels-1; level > 0; --level)
+        {
+          unsigned int counter = 0;
+          std::vector<types::global_dof_index> global_level_dof_indices;
+          std::vector<types::global_dof_index> global_level_dof_indices_remote;
+          std::vector<types::global_dof_index> ghosted_level_dofs;
+          std::vector<types::global_dof_index> global_level_dof_indices_l0;
+          std::vector<types::global_dof_index> ghosted_level_dofs_l0;
+
+          // step 2.1: loop over the cells on the coarse side
+          typename dealii::DoFHandler<dim>::cell_iterator cell, endc = mg_dof.end(level-1);
+          for (cell = mg_dof.begin(level-1); cell != endc; ++cell)
+            {
+              // need to look into a cell if it has children and it is locally
+              // owned
+              if (!cell->has_children())
+                continue;
+
+              bool consider_cell = false;
+              if (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id
+                  || cell->level_subdomain_id()==tria.locally_owned_subdomain()
+                 )
+                consider_cell = true;
+
+              // due to the particular way we store DoF indices (via children),
+              // we also need to add the DoF indices for coarse cells where we
+              // own at least one child
+              bool cell_is_remote = !consider_cell;
+              for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+                if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain())
+                  {
+                    consider_cell = true;
+                    break;
+                  }
+
+              if (!consider_cell)
+                continue;
+
+              // step 2.2: loop through children and append the dof indices to
+              // the appropriate list. We need separate lists for the owned
+              // coarse cell case (which will be part of
+              // restriction/prolongation between level-1 and level) and the
+              // remote case (which needs to store DoF indices for the
+              // operations between level and level+1).
+              AssertDimension(cell->n_children(),
+                              GeometryInfo<dim>::max_children_per_cell);
+              std::vector<types::global_dof_index> &next_indices =
+                cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices;
+              const std::size_t start_index = next_indices.size();
+              next_indices.resize(start_index + elem_info.n_child_cell_dofs,
+                                  numbers::invalid_dof_index);
+              for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
+                {
+                  if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
+                      tria.locally_owned_subdomain())
+                    continue;
+                  cell->child(c)->get_mg_dof_indices(local_dof_indices);
+
+                  const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
+                  for (unsigned int i=0; i<local_dof_indices.size(); ++i)
+                    if (!owned_level_dofs.is_element(local_dof_indices[i]))
+                      ghosted_level_dofs.push_back(local_dof_indices[i]);
+
+                  add_child_indices<dim>(c, fe.dofs_per_cell - fe.dofs_per_vertex,
+                                         fe.degree, elem_info.shape_info.lexicographic_numbering,
+                                         local_dof_indices,
+                                         &next_indices[start_index]);
+
+                  // step 2.3 store the connectivity to the parent
+                  if (cell->child(c)->has_children() &&
+                      (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id
+                       || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()
+                      ))
+                    {
+                      const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()];
+                      AssertIndexRange(child_index, parent_child_connect[level].size());
+                      unsigned int parent_index = counter;
+                      // remote cells, i.e., cells where we work on a further
+                      // level but are not treated on the current level, need to
+                      // be placed at the end of the list; however, we do not
+                      // yet know the exact position in the array, so shift
+                      // their parent index by the number of cells so we can set
+                      // the correct number after the end of this loop
+                      if (cell_is_remote)
+                        parent_index = start_index/elem_info.n_child_cell_dofs + tria.n_cells(level);
+                      parent_child_connect[level][child_index] =
+                        std::make_pair(parent_index, c);
+                      AssertIndexRange(mg_dof.get_fe().dofs_per_cell,
+                                       static_cast<unsigned short>(-1));
+
+                      // set Dirichlet boundary conditions (as a list of
+                      // constrained DoFs) for the child
+                      if (mg_constrained_dofs != 0)
+                        for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
+                          if (mg_constrained_dofs->is_boundary_index(level,
+                                                                     local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]]))
+                            dirichlet_indices[level][child_index].push_back(i);
+                    }
+                }
+              if (!cell_is_remote)
+                {
+                  AssertIndexRange(static_cast<unsigned int>(cell->index()),
+                                   coarse_level_indices[level-1].size());
+                  coarse_level_indices[level-1][cell->index()] = counter++;
+                }
+
+              // step 2.4: include indices for the coarsest cells. we still
+              // insert the indices as if they were from a child in order to use
+              // the same code (the coarsest level does not matter much in terms
+              // of memory, so we gain in code simplicity)
+              if (level == 1 && !cell_is_remote)
+                {
+                  cell->get_mg_dof_indices(local_dof_indices);
+
+                  const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
+                  for (unsigned int i=0; i<local_dof_indices.size(); ++i)
+                    if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
+                      ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
+
+                  const std::size_t start_index = global_level_dof_indices_l0.size();
+                  global_level_dof_indices_l0.resize(start_index+elem_info.n_child_cell_dofs,
+                                                     numbers::invalid_dof_index);
+                  add_child_indices<dim>(0, fe.dofs_per_cell - fe.dofs_per_vertex,
+                                         fe.degree, elem_info.shape_info.lexicographic_numbering,
+                                         local_dof_indices,
+                                         &global_level_dof_indices_l0[start_index]);
+
+                  dirichlet_indices[0].push_back(std::vector<unsigned short>());
+                  if (mg_constrained_dofs != 0)
+                    for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
+                      if (mg_constrained_dofs->is_boundary_index(0, local_dof_indices[elem_info.shape_info.lexicographic_numbering[i]]))
+                        dirichlet_indices[0].back().push_back(i);
+                }
+            }
+
+          // step 2.5: store information about the current level and prepare the
+          // Dirichlet indices and parent-child relationship for the next
+          // coarser level
+          AssertDimension(counter*elem_info.n_child_cell_dofs, global_level_dof_indices.size());
+          n_owned_level_cells[level-1] = counter;
+          dirichlet_indices[level-1].resize(counter);
+          parent_child_connect[level-1].
+          resize(counter, std::make_pair(numbers::invalid_unsigned_int,
+                                         numbers::invalid_unsigned_int));
+
+          // step 2.6: put the cells with remotely owned parent to the end of
+          // the list (these are needed for the transfer from level to level+1
+          // but not for the transfer from level-1 to level).
+          if (level < n_levels-1)
+            for (std::vector<std::pair<unsigned int,unsigned int> >::iterator
+                 i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i)
+              if (i->first >= tria.n_cells(level))
+                {
+                  i->first -= tria.n_cells(level);
+                  i->first += counter;
+                }
+
+          // step 2.7: Initialize the ghosted vector
+          const parallel::Triangulation<dim,dim> *ptria =
+            (dynamic_cast<const parallel::Triangulation<dim,dim>*> (&tria));
+          const MPI_Comm communicator =
+            ptria != 0 ? ptria->get_communicator() : MPI_COMM_SELF;
+
+          reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(level),
+                                 ghosted_level_dofs, communicator,
+                                 ghosted_level_vector[level],
+                                 copy_indices_global_mine[level]);
+
+          copy_indices_to_mpi_local_numbers(*ghosted_level_vector[level].get_partitioner(),
+                                            global_level_dof_indices,
+                                            global_level_dof_indices_remote,
+                                            level_dof_indices[level]);
+          // step 2.8: Initialize the ghosted vector for level 0
+          if (level == 1)
+            {
+              for (unsigned int i = 0; i<parent_child_connect[0].size(); ++i)
+                parent_child_connect[0][i] = std::make_pair(i, 0U);
+
+              reinit_ghosted_vector (mg_dof.locally_owned_mg_dofs(0),
+                                     ghosted_level_dofs_l0, communicator,
+                                     ghosted_level_vector[0],
+                                     copy_indices_global_mine[0]);
+
+              copy_indices_to_mpi_local_numbers(*ghosted_level_vector[0].get_partitioner(),
+                                                global_level_dof_indices_l0,
+                                                std::vector<types::global_dof_index>(),
+                                                level_dof_indices[0]);
+            }
+        }
+
+      // ---------------------- 3. compute weights to make restriction additive
+
+      // get the valence of the individual components and compute the weights as
+      // the inverse of the valence
+      weights_on_refined.resize(n_levels-1);
+      for (unsigned int level = 1; level<n_levels; ++level)
+        {
+          ghosted_level_vector[level] = 0;
+          for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
+            for (unsigned int j=0; j<elem_info.n_child_cell_dofs; ++j)
+              ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+j]) += Number(1.);
+          ghosted_level_vector[level].compress(VectorOperation::add);
+          ghosted_level_vector[level].update_ghost_values();
+
+          const unsigned int vec_size = VectorizedArray<Number>::n_array_elements;
+          std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
+          degree_to_3[0] = 0;
+          for (unsigned int i=1; i<n_child_dofs_1d-1; ++i)
+            degree_to_3[i] = 1;
+          degree_to_3.back() = 2;
+
+          // we only store 3^dim weights because all dofs on a line have the
+          // same valence, and all dofs on a quad have the same valence.
+          weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*Utilities::fixed_power<dim>(3));
+          for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
+            {
+              const unsigned int comp = c/vec_size;
+              const unsigned int v = c%vec_size;
+
+              for (unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
+                for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
+                  {
+                    unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
+                    for (unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
+                      weights_on_refined[level-1][comp*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]][v] = Number(1.)/
+                          ghosted_level_vector[level].local_element(level_dof_indices[level][elem_info.n_child_cell_dofs*c+m]);
+                  }
+            }
+        }
+
+    }
+
   }
 }
 
diff --git a/source/multigrid/mg_transfer_internal.inst.in b/source/multigrid/mg_transfer_internal.inst.in
index bab949b153..96c231acbd 100644
--- a/source/multigrid/mg_transfer_internal.inst.in
+++ b/source/multigrid/mg_transfer_internal.inst.in
@@ -35,3 +35,55 @@ for (deal_II_dimension : DIMENSIONS;
     \}
     \}
 }
+
+
+for (deal_II_dimension : DIMENSIONS)
+{
+    namespace internal
+    \{
+    namespace MGTransfer
+    \{
+
+    template
+    unsigned int
+    compute_shift_within_children<deal_II_dimension>(const unsigned int,
+            const unsigned int,
+            const unsigned int);
+    \}
+    \}
+}
+
+for (S : REAL_SCALARS)
+{
+    namespace internal
+    \{
+    namespace MGTransfer
+    \{
+
+    template
+    struct ElementInfo<S>;
+    \}
+    \}
+}
+
+for (deal_II_dimension : DIMENSIONS; S : REAL_SCALARS)
+{
+    namespace internal
+    \{
+    namespace MGTransfer
+    \{
+
+    template
+    void setup_transfer<deal_II_dimension>(const dealii::DoFHandler<deal_II_dimension>&,
+                                           const MGConstrainedDoFs*,
+                                           ElementInfo<S>&,
+                                           std::vector<std::vector<unsigned int> >&,
+                                           std::vector<std::vector<std::pair<unsigned int,unsigned int> > >&,
+                                           std::vector<unsigned int>&,
+                                           std::vector<std::vector<std::vector<unsigned short> > >&,
+                                           std::vector<AlignedVector<VectorizedArray<S> > >&,
+                                           std::vector<std::vector<std::pair<unsigned int, unsigned int> > >&,
+                                           MGLevelObject<LinearAlgebra::distributed::Vector<S> >&);
+    \}
+    \}
+}
diff --git a/source/multigrid/mg_transfer_matrix_free.cc b/source/multigrid/mg_transfer_matrix_free.cc
index f8117424b3..e23881f1b3 100644
--- a/source/multigrid/mg_transfer_matrix_free.cc
+++ b/source/multigrid/mg_transfer_matrix_free.cc
@@ -22,10 +22,10 @@
 #include <deal.II/grid/tria_iterator.h>
 #include <deal.II/dofs/dof_tools.h>
 #include <deal.II/fe/fe.h>
-#include <deal.II/fe/fe_tools.h>
 #include <deal.II/dofs/dof_accessor.h>
 #include <deal.II/multigrid/mg_tools.h>
 #include <deal.II/multigrid/mg_transfer_matrix_free.h>
+#include <deal.II/multigrid/mg_transfer_internal.h>
 
 #include <deal.II/matrix_free/shape_info.h>
 #include <deal.II/matrix_free/fe_evaluation.h>
@@ -91,453 +91,29 @@ void MGTransferMatrixFree<dim,Number>::clear ()
 }
 
 
-
-namespace
-{
-  // given the collection of child cells in lexicographic ordering as seen
-  // from the parent, compute the first index of the given child
-  template <int dim>
-  unsigned int
-  compute_shift_within_children(const unsigned int child,
-                                const unsigned int fe_shift_1d,
-                                const unsigned int fe_degree)
-  {
-    // we put the degrees of freedom of all child cells in
-    // lexicographic ordering
-    unsigned int c_tensor_index[dim];
-    unsigned int tmp = child;
-    for (unsigned int d=0; d<dim; ++d)
-      {
-        c_tensor_index[d] = tmp % 2;
-        tmp /= 2;
-      }
-    const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
-    unsigned int factor = 1;
-    unsigned int shift = fe_shift_1d * c_tensor_index[0];
-    for (unsigned int d=1; d<dim; ++d)
-      {
-        factor *= n_child_dofs_1d;
-        shift = shift + factor * fe_shift_1d * c_tensor_index[d];
-      }
-    return shift;
-  }
-
-
-
-  // puts the indices on the given child cell in lexicographic ordering with
-  // respect to the collection of all child cells as seen from the parent
-  template <int dim>
-  void add_child_indices(const unsigned int child,
-                         const unsigned int fe_shift_1d,
-                         const unsigned int fe_degree,
-                         const std::vector<unsigned int> &lexicographic_numbering,
-                         const std::vector<types::global_dof_index> &local_dof_indices,
-                         types::global_dof_index *target_indices)
-  {
-    const unsigned int n_child_dofs_1d = fe_degree + 1 + fe_shift_1d;
-    const unsigned int shift =
-      compute_shift_within_children<dim>(child, fe_shift_1d, fe_degree);
-    const unsigned int n_components =
-      local_dof_indices.size()/Utilities::fixed_power<dim>(fe_degree+1);
-    types::global_dof_index *indices = target_indices + shift;
-    const unsigned int n_scalar_cell_dofs = Utilities::fixed_power<dim>(n_child_dofs_1d);
-    for (unsigned int c=0, m=0; c<n_components; ++c)
-      for (unsigned int k=0; k<(dim>2 ? (fe_degree+1) : 1); ++k)
-        for (unsigned int j=0; j<(dim>1 ? (fe_degree+1) : 1); ++j)
-          for (unsigned int i=0; i<(fe_degree+1); ++i, ++m)
-            {
-              const unsigned int index = c*n_scalar_cell_dofs+k*n_child_dofs_1d*
-                                         n_child_dofs_1d+j*n_child_dofs_1d+i;
-              Assert(indices[index] == numbers::invalid_dof_index ||
-                     indices[index] == local_dof_indices[lexicographic_numbering[m]],
-                     ExcInternalError());
-              indices[index] = local_dof_indices[lexicographic_numbering[m]];
-            }
-  }
-
-
-
-  // initialize the vectors needed for the transfer (and merge with the
-  // content in copy_indices_global_mine)
-  template <typename Number>
-  void
-  reinit_ghosted_vector(const IndexSet &locally_owned,
-                        std::vector<types::global_dof_index> &ghosted_level_dofs,
-                        const MPI_Comm &communicator,
-                        LinearAlgebra::distributed::Vector<Number> &ghosted_level_vector,
-                        std::vector<std::pair<unsigned int,unsigned int> > &copy_indices_global_mine)
-  {
-    std::sort(ghosted_level_dofs.begin(), ghosted_level_dofs.end());
-    IndexSet ghosted_dofs(locally_owned.size());
-    ghosted_dofs.add_indices(ghosted_level_dofs.begin(),
-                             std::unique(ghosted_level_dofs.begin(),
-                                         ghosted_level_dofs.end()));
-    ghosted_dofs.compress();
-
-    // Add possible ghosts from the previous content in the vector
-    if (ghosted_level_vector.size() == locally_owned.size())
-      {
-        // shift the local number of the copy indices according to the new
-        // partitioner that we are going to use for the vector
-        const std_cxx11::shared_ptr<const Utilities::MPI::Partitioner> part
-          = ghosted_level_vector.get_partitioner();
-        ghosted_dofs.add_indices(part->ghost_indices());
-        for (unsigned int i=0; i<copy_indices_global_mine.size(); ++i)
-          copy_indices_global_mine[i].second =
-            locally_owned.n_elements() +
-            ghosted_dofs.index_within_set(part->local_to_global(copy_indices_global_mine[i].second));
-      }
-    ghosted_level_vector.reinit(locally_owned, ghosted_dofs, communicator);
-  }
-
-  // Transform the ghost indices to local index space for the vector
-  void
-  copy_indices_to_mpi_local_numbers(const Utilities::MPI::Partitioner &part,
-                                    const std::vector<types::global_dof_index> &mine,
-                                    const std::vector<types::global_dof_index> &remote,
-                                    std::vector<unsigned int> &localized_indices)
-  {
-    localized_indices.resize(mine.size()+remote.size(),
-                             numbers::invalid_unsigned_int);
-    for (unsigned int i=0; i<mine.size(); ++i)
-      if (mine[i] != numbers::invalid_dof_index)
-        localized_indices[i] = part.global_to_local(mine[i]);
-
-    for (unsigned int i=0; i<remote.size(); ++i)
-      if (remote[i] != numbers::invalid_dof_index)
-        localized_indices[i+mine.size()] = part.global_to_local(remote[i]);
-  }
-}
-
-
-
 template <int dim, typename Number>
 void MGTransferMatrixFree<dim,Number>::build
 (const DoFHandler<dim,dim>  &mg_dof)
 {
   this->fill_and_communicate_copy_indices(mg_dof);
 
-  // we collect all child DoFs of a mother cell together. For faster
-  // tensorized operations, we align the degrees of freedom
-  // lexicographically. We distinguish FE_Q elements and FE_DGQ elements
-
-  const Triangulation<dim> &tria = mg_dof.get_triangulation();
-
-  // ---------------------------- 1. Extract 1D info about the finite element
-  // step 1.1: create a 1D copy of the finite element from FETools where we
-  // substitute the template argument
-  AssertDimension(mg_dof.get_fe().n_base_elements(), 1);
-  std::string fe_name = mg_dof.get_fe().base_element(0).get_name();
-  {
-    const std::size_t template_starts = fe_name.find_first_of('<');
-    Assert (fe_name[template_starts+1] == (dim==1?'1':(dim==2?'2':'3')),
-            ExcInternalError());
-    fe_name[template_starts+1] = '1';
-  }
-  std_cxx11::shared_ptr<FiniteElement<1> > fe_1d
-  (FETools::get_fe_by_name<1,1>(fe_name));
-  const FiniteElement<1> &fe = *fe_1d;
-  unsigned int n_child_dofs_1d = numbers::invalid_unsigned_int;
-
-  {
-    // currently, we have only FE_Q and FE_DGQ type elements implemented
-    n_components = mg_dof.get_fe().element_multiplicity(0);
-    AssertDimension(Utilities::fixed_power<dim>(fe.dofs_per_cell)*n_components,
-                    mg_dof.get_fe().dofs_per_cell);
-    AssertDimension(fe.degree, mg_dof.get_fe().degree);
-    fe_degree = fe.degree;
-    element_is_continuous = fe.dofs_per_vertex > 0;
-    Assert(fe.dofs_per_vertex < 2, ExcNotImplemented());
-
-    // step 1.2: get renumbering of 1D basis functions to lexicographic
-    // numbers. The distinction according to fe.dofs_per_vertex is to support
-    // both continuous and discontinuous bases.
-    std::vector<unsigned int> renumbering(fe.dofs_per_cell);
-    {
-      AssertIndexRange(fe.dofs_per_vertex, 2);
-      renumbering[0] = 0;
-      for (unsigned int i=0; i<fe.dofs_per_line; ++i)
-        renumbering[i+fe.dofs_per_vertex] =
-          GeometryInfo<1>::vertices_per_cell*fe.dofs_per_vertex + i;
-      if (fe.dofs_per_vertex > 0)
-        renumbering[fe.dofs_per_cell-fe.dofs_per_vertex] = fe.dofs_per_vertex;
-    }
-
-    // step 1.3: create a 1D quadrature formula from the finite element that
-    // collects the support points of the basis functions on the two children.
-    std::vector<Point<1> > basic_support_points = fe.get_unit_support_points();
-    Assert(fe.dofs_per_vertex == 0 || fe.dofs_per_vertex == 1,
-           ExcNotImplemented());
-    std::vector<Point<1> > points_refined(fe.dofs_per_vertex > 0 ?
-                                          (2 * fe.dofs_per_cell - 1) :
-                                          (2 * fe.dofs_per_cell));
-    const unsigned int shift = fe.dofs_per_cell - fe.dofs_per_vertex;
-    for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
-      for (unsigned int j=0; j<basic_support_points.size(); ++j)
-        points_refined[shift*c+j][0] =
-          c*0.5 + 0.5 * basic_support_points[renumbering[j]][0];
-
-    n_child_dofs_1d = points_refined.size();
-    n_child_cell_dofs = n_components*Utilities::fixed_power<dim>(n_child_dofs_1d);
-
-    // step 1.4: evaluate the polynomials and store the data in ShapeInfo
-    const Quadrature<1> quadrature(points_refined);
-    shape_info.reinit(quadrature, mg_dof.get_fe(), 0);
-
-    for (unsigned int c=0; c<GeometryInfo<1>::max_children_per_cell; ++c)
-      for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
-        for (unsigned int j=0; j<fe.dofs_per_cell; ++j)
-          Assert(std::abs(shape_info.shape_values[i*n_child_dofs_1d+j+c*shift][0] -
-                          fe.get_prolongation_matrix(c)(renumbering[j],renumbering[i]))
-                 < std::max(2.*(double)std::numeric_limits<Number>::epsilon(),1e-12),
-                 ExcInternalError());
-  }
-
-  // -------------- 2. Extract and match dof indices between child and parent
-  const unsigned int n_levels = tria.n_global_levels();
-  level_dof_indices.resize(n_levels);
-  parent_child_connect.resize(n_levels-1);
-  n_owned_level_cells.resize(n_levels-1);
-  std::vector<std::vector<unsigned int> > coarse_level_indices(n_levels-1);
-  for (unsigned int level=0; level<std::min(tria.n_levels(),n_levels-1); ++level)
-    coarse_level_indices[level].resize(tria.n_raw_cells(level),
-                                       numbers::invalid_unsigned_int);
-  std::vector<types::global_dof_index> local_dof_indices(mg_dof.get_fe().dofs_per_cell);
-  dirichlet_indices.resize(n_levels-1);
-
-  // We use the vectors stored ghosted_level_vector in the base class for
-  // keeping ghosted transfer indices. To avoid keeping two very similar
-  // vectors, we merge them here.
-  if (this->ghosted_level_vector.max_level() != n_levels-1)
-    this->ghosted_level_vector.resize(0, n_levels-1);
-
-  for (unsigned int level=n_levels-1; level > 0; --level)
-    {
-      unsigned int counter = 0;
-      std::vector<types::global_dof_index> global_level_dof_indices;
-      std::vector<types::global_dof_index> global_level_dof_indices_remote;
-      std::vector<types::global_dof_index> ghosted_level_dofs;
-      std::vector<types::global_dof_index> global_level_dof_indices_l0;
-      std::vector<types::global_dof_index> ghosted_level_dofs_l0;
-
-      // step 2.1: loop over the cells on the coarse side
-      typename DoFHandler<dim>::cell_iterator cell, endc = mg_dof.end(level-1);
-      for (cell = mg_dof.begin(level-1); cell != endc; ++cell)
-        {
-          // need to look into a cell if it has children and it is locally owned
-          if (!cell->has_children())
-            continue;
-
-          bool consider_cell = false;
-          if (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id
-              || cell->level_subdomain_id()==tria.locally_owned_subdomain()
-             )
-            consider_cell = true;
-
-          // due to the particular way we store DoF indices (via children), we
-          // also need to add the DoF indices for coarse cells where we own at
-          // least one child
-          bool cell_is_remote = !consider_cell;
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            if (cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain())
-              {
-                consider_cell = true;
-                break;
-              }
-
-          if (!consider_cell)
-            continue;
-
-          // step 2.2: loop through children and append the dof indices to the
-          // appropriate list. We need separate lists for the owned coarse
-          // cell case (which will be part of restriction/prolongation between
-          // level-1 and level) and the remote case (which needs to store DoF
-          // indices for the operations between level and level+1).
-          AssertDimension(cell->n_children(),
-                          GeometryInfo<dim>::max_children_per_cell);
-          std::vector<types::global_dof_index> &next_indices =
-            cell_is_remote ? global_level_dof_indices_remote : global_level_dof_indices;
-          const std::size_t start_index = next_indices.size();
-          next_indices.resize(start_index + n_child_cell_dofs,
-                              numbers::invalid_dof_index);
-          for (unsigned int c=0; c<GeometryInfo<dim>::max_children_per_cell; ++c)
-            {
-              if (cell_is_remote && cell->child(c)->level_subdomain_id() !=
-                  tria.locally_owned_subdomain())
-                continue;
-              cell->child(c)->get_mg_dof_indices(local_dof_indices);
-
-              const IndexSet &owned_level_dofs = mg_dof.locally_owned_mg_dofs(level);
-              for (unsigned int i=0; i<local_dof_indices.size(); ++i)
-                if (!owned_level_dofs.is_element(local_dof_indices[i]))
-                  ghosted_level_dofs.push_back(local_dof_indices[i]);
-
-              add_child_indices<dim>(c, fe.dofs_per_cell - fe.dofs_per_vertex,
-                                     fe.degree, shape_info.lexicographic_numbering,
-                                     local_dof_indices,
-                                     &next_indices[start_index]);
-
-              // step 2.3 store the connectivity to the parent
-              if (cell->child(c)->has_children() &&
-                  (tria.locally_owned_subdomain()==numbers::invalid_subdomain_id
-                   || cell->child(c)->level_subdomain_id()==tria.locally_owned_subdomain()
-                  ))
-                {
-                  const unsigned int child_index = coarse_level_indices[level][cell->child(c)->index()];
-                  AssertIndexRange(child_index, parent_child_connect[level].size());
-                  unsigned int parent_index = counter;
-                  // remote cells, i.e., cells where we work on a further
-                  // level but are not treated on the current level, need to
-                  // be placed at the end of the list; however, we do not yet
-                  // know the exact position in the array, so shift their
-                  // parent index by the number of cells so we can set the
-                  // correct number after the end of this loop
-                  if (cell_is_remote)
-                    parent_index = start_index/n_child_cell_dofs + tria.n_cells(level);
-                  parent_child_connect[level][child_index] =
-                    std::make_pair(parent_index, c);
-                  AssertIndexRange(mg_dof.get_fe().dofs_per_cell,
-                                   static_cast<unsigned short>(-1));
-
-                  // set Dirichlet boundary conditions (as a list of
-                  // constrained DoFs) for the child
-                  if (this->mg_constrained_dofs != 0)
-                    for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
-                      if (this->mg_constrained_dofs->is_boundary_index(level, local_dof_indices[shape_info.lexicographic_numbering[i]]))
-                        dirichlet_indices[level][child_index].push_back(i);
-                }
-            }
-          if (!cell_is_remote)
-            {
-              AssertIndexRange(static_cast<unsigned int>(cell->index()),
-                               coarse_level_indices[level-1].size());
-              coarse_level_indices[level-1][cell->index()] = counter++;
-            }
-
-          // step 2.4: include indices for the coarsest cells. we still insert
-          // the indices as if they were from a child in order to use the same
-          // code (the coarsest level does not matter much in terms of memory,
-          // so we gain in code simplicity)
-          if (level == 1 && !cell_is_remote)
-            {
-              cell->get_mg_dof_indices(local_dof_indices);
-
-              const IndexSet &owned_level_dofs_l0 = mg_dof.locally_owned_mg_dofs(0);
-              for (unsigned int i=0; i<local_dof_indices.size(); ++i)
-                if (!owned_level_dofs_l0.is_element(local_dof_indices[i]))
-                  ghosted_level_dofs_l0.push_back(local_dof_indices[i]);
-
-              const std::size_t start_index = global_level_dof_indices_l0.size();
-              global_level_dof_indices_l0.resize(start_index+n_child_cell_dofs,
-                                                 numbers::invalid_dof_index);
-              add_child_indices<dim>(0, fe.dofs_per_cell - fe.dofs_per_vertex,
-                                     fe.degree, shape_info.lexicographic_numbering,
-                                     local_dof_indices,
-                                     &global_level_dof_indices_l0[start_index]);
-
-              dirichlet_indices[0].push_back(std::vector<unsigned short>());
-              if (this->mg_constrained_dofs != 0)
-                for (unsigned int i=0; i<mg_dof.get_fe().dofs_per_cell; ++i)
-                  if (this->mg_constrained_dofs->is_boundary_index(0, local_dof_indices[shape_info.lexicographic_numbering[i]]))
-                    dirichlet_indices[0].back().push_back(i);
-            }
-        }
-
-      // step 2.5: store information about the current level and prepare the
-      // Dirichlet indices and parent-child relationship for the next coarser
-      // level
-      AssertDimension(counter*n_child_cell_dofs, global_level_dof_indices.size());
-      n_owned_level_cells[level-1] = counter;
-      dirichlet_indices[level-1].resize(counter);
-      parent_child_connect[level-1].
-      resize(counter, std::make_pair(numbers::invalid_unsigned_int,
-                                     numbers::invalid_unsigned_int));
-
-      // step 2.6: put the cells with remotely owned parent to the end of the
-      // list (these are needed for the transfer from level to level+1 but not
-      // for the transfer from level-1 to level).
-      if (level < n_levels-1)
-        for (std::vector<std::pair<unsigned int,unsigned int> >::iterator
-             i=parent_child_connect[level].begin(); i!=parent_child_connect[level].end(); ++i)
-          if (i->first >= tria.n_cells(level))
-            {
-              i->first -= tria.n_cells(level);
-              i->first += counter;
-            }
-
-      // step 2.7: Initialize the ghosted vector
-      const parallel::Triangulation<dim,dim> *ptria =
-        (dynamic_cast<const parallel::Triangulation<dim,dim>*> (&tria));
-      const MPI_Comm communicator =
-        ptria != 0 ? ptria->get_communicator() : MPI_COMM_SELF;
-
-      reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(level),
-                            ghosted_level_dofs, communicator,
-                            this->ghosted_level_vector[level],
-                            this->copy_indices_global_mine[level]);
-
-      copy_indices_to_mpi_local_numbers(*this->ghosted_level_vector[level].get_partitioner(),
-                                        global_level_dof_indices,
-                                        global_level_dof_indices_remote,
-                                        level_dof_indices[level]);
-
-      // step 2.8: Initialize the ghosted vector for level 0
-      if (level == 1)
-        {
-          for (unsigned int i = 0; i<parent_child_connect[0].size(); ++i)
-            parent_child_connect[0][i] = std::make_pair(i, 0U);
-
-          reinit_ghosted_vector(mg_dof.locally_owned_mg_dofs(0),
-                                ghosted_level_dofs_l0, communicator,
-                                this->ghosted_level_vector[0],
-                                this->copy_indices_global_mine[0]);
-
-          copy_indices_to_mpi_local_numbers(*this->ghosted_level_vector[0].get_partitioner(),
-                                            global_level_dof_indices_l0,
-                                            std::vector<types::global_dof_index>(),
-                                            level_dof_indices[0]);
-        }
-    }
-
-  // ------------------------ 3. compute weights to make restriction additive
-  //
-  // get the valence of the individual components and compute the weights as
-  // the inverse of the valence
-  weights_on_refined.resize(n_levels-1);
-  for (unsigned int level = 1; level<n_levels; ++level)
-    {
-      this->ghosted_level_vector[level] = 0;
-      for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
-        for (unsigned int j=0; j<n_child_cell_dofs; ++j)
-          this->ghosted_level_vector[level].local_element(level_dof_indices[level][n_child_cell_dofs*c+j]) += Number(1.);
-      this->ghosted_level_vector[level].compress(VectorOperation::add);
-      this->ghosted_level_vector[level].update_ghost_values();
-
-      const unsigned int vec_size = VectorizedArray<Number>::n_array_elements;
-      std::vector<unsigned int> degree_to_3 (n_child_dofs_1d);
-      degree_to_3[0] = 0;
-      for (unsigned int i=1; i<n_child_dofs_1d-1; ++i)
-        degree_to_3[i] = 1;
-      degree_to_3.back() = 2;
-
-      // we only store 3^dim weights because all dofs on a line have the same
-      // valence, and all dofs on a quad have the same valence.
-      weights_on_refined[level-1].resize(((n_owned_level_cells[level-1]+vec_size-1)/vec_size)*Utilities::fixed_power<dim>(3));
-      for (unsigned int c=0; c<n_owned_level_cells[level-1]; ++c)
-        {
-          const unsigned int comp = c/vec_size;
-          const unsigned int v = c%vec_size;
-
-          for (unsigned int k=0, m=0; k<(dim>2 ? n_child_dofs_1d : 1); ++k)
-            for (unsigned int j=0; j<(dim>1 ? n_child_dofs_1d : 1); ++j)
-              {
-                unsigned int shift = 9*degree_to_3[k] + 3*degree_to_3[j];
-                for (unsigned int i=0; i<n_child_dofs_1d; ++i, ++m)
-                  weights_on_refined[level-1][comp*Utilities::fixed_power<dim>(3)+shift+degree_to_3[i]][v] = Number(1.)/
-                      this->ghosted_level_vector[level].local_element(level_dof_indices[level][n_child_cell_dofs*c+m]);
-              }
-        }
-    }
+  internal::MGTransfer::ElementInfo<Number> elem_info;
+  internal::MGTransfer::setup_transfer<dim,Number>(mg_dof,
+                                                   (const MGConstrainedDoFs *)this->mg_constrained_dofs,
+                                                   elem_info,
+                                                   level_dof_indices,
+                                                   parent_child_connect,
+                                                   n_owned_level_cells,
+                                                   dirichlet_indices,
+                                                   weights_on_refined,
+                                                   this->copy_indices_global_mine,
+                                                   this->ghosted_level_vector);
+  // unpack element info data
+  fe_degree                = elem_info.fe_degree;
+  element_is_continuous    = elem_info.element_is_continuous;
+  n_components             = elem_info.n_components;
+  n_child_cell_dofs        = elem_info.n_child_cell_dofs;
+  shape_info               = elem_info.shape_info;
 
   evaluation_data.resize(3*n_child_cell_dofs);
 }
@@ -764,7 +340,7 @@ void MGTransferMatrixFree<dim,Number>
       // read from source vector
       for (unsigned int v=0; v<n_chunks; ++v)
         {
-          const unsigned int shift = compute_shift_within_children<dim>
+          const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
                                      (parent_child_connect[to_level-1][cell+v].second,
                                       degree+1-element_is_continuous, degree);
           const unsigned int *indices = &level_dof_indices[to_level-1][parent_child_connect[to_level-1][cell+v].first*n_child_cell_dofs+shift];
@@ -897,7 +473,7 @@ void MGTransferMatrixFree<dim,Number>
       // write into dst vector
       for (unsigned int v=0; v<n_chunks; ++v)
         {
-          const unsigned int shift = compute_shift_within_children<dim>
+          const unsigned int shift = internal::MGTransfer::compute_shift_within_children<dim>
                                      (parent_child_connect[from_level-1][cell+v].second,
                                       degree+1-element_is_continuous, degree);
           AssertIndexRange(parent_child_connect[from_level-1][cell+v].first*