From: heister <heister@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Mon, 1 Oct 2012 15:23:44 +0000 (+0000)
Subject: merge from mainline
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=3f33529ca6da777b29784f3b4d2fd294b745bee8;p=dealii-svn.git

merge from mainline

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

3f33529ca6da777b29784f3b4d2fd294b745bee8
diff --cc deal.II/include/deal.II/dofs/dof_tools.h
index 1c568de5b5,3565e81af0..5c71d01791
--- a/deal.II/include/deal.II/dofs/dof_tools.h
+++ b/deal.II/include/deal.II/dofs/dof_tools.h
@@@ -1234,22 -1250,61 +1250,61 @@@ namespace DoFTool
    template <int dim, int spacedim>
    void
    extract_dofs (const hp::DoFHandler<dim,spacedim>   &dof_handler,
-                 const std::vector<bool> &select,
-                 std::vector<bool>       &selected_dofs,
-                 const bool               blocks = false);
+                 const ComponentMask     &component_mask,
+                 std::vector<bool>       &selected_dofs);
  
-                                    /**
-                                     * Do the same thing as
-                                     * extract_dofs() for one level
-                                     * of a multi-grid DoF numbering.
-                                     */
+   /**
+    * This function is the equivalent to the DoFTools::extract_dofs() functions above
+    * except that the selection of which degrees of freedom to extract is not done
+    * based on components (see @ref GlossComponent) but instead based on whether they
+    * are part of a particular block (see @ref GlossBlock). Consequently, the second
+    * argument is not a ComponentMask but a BlockMask object.
+    *
+    * @param dof_handler The DoFHandler object from which to extract degrees of freedom
+    * @param block_mask The block mask that describes which blocks to consider (see
+    *       @ref GlossBlockMask)
+    * @param selected_dofs A vector of length DoFHandler::n_dofs() in which those
+    *       entries are true that correspond to the selected blocks.
+    */
+   template <int dim, int spacedim>
+   void
+   extract_dofs (const DoFHandler<dim,spacedim>   &dof_handler,
+                 const BlockMask &block_mask,
+                 std::vector<bool>       &selected_dofs);
+ 
+   /**
+    * The same function as above,
+    * but for a hp::DoFHandler.
+    */
+   template <int dim, int spacedim>
+   void
+   extract_dofs (const hp::DoFHandler<dim,spacedim>   &dof_handler,
+                 const BlockMask        &block_mask,
+                 std::vector<bool>       &selected_dofs);
+ 
+   /**
+    * Do the same thing as the corresponding
+    * extract_dofs() functino for one level
+    * of a multi-grid DoF numbering.
+    */
 -  template <int dim, int spacedim>
++  template <class DH>
+   void
+   extract_level_dofs (const unsigned int       level,
 -                      const MGDoFHandler<dim,spacedim> &dof,
++                      const DH &dof,
+                       const ComponentMask     &component_mask,
+                       std::vector<bool>       &selected_dofs);
+ 
+   /**
+    * Do the same thing as the corresponding
+    * extract_dofs() functino for one level
+    * of a multi-grid DoF numbering.
+    */
 -  template <int dim, int spacedim>
 +  template <class DH>
    void
    extract_level_dofs (const unsigned int       level,
 -                      const MGDoFHandler<dim,spacedim> &dof,
 +                      const DH &dof,
-                       const std::vector<bool> &select,
-                       std::vector<bool>       &selected_dofs,
-                       const bool               blocks = false);
+                       const BlockMask     &component_mask,
+                       std::vector<bool>       &selected_dofs);
  
    /**
     * Extract all degrees of freedom
diff --cc deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
index 4288c6d92c,e011893a8e..5b90025982
--- a/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
+++ b/deal.II/include/deal.II/multigrid/mg_constrained_dofs.h
@@@ -59,14 -58,9 +59,14 @@@ class MGConstrainedDoFs : public Subscr
                                        */
      template <int dim, int spacedim>
      void initialize(const MGDoFHandler<dim,spacedim>& dof,
- 		            const typename FunctionMap<dim>::type& function_map,
- 		            const std::vector<bool>& component_mask = std::vector<bool>());
+                     const typename FunctionMap<dim>::type& function_map,
+                     const ComponentMask &component_mask = ComponentMask());
  
 +    template <int dim, int spacedim>
 +    void initialize(const DoFHandler<dim,spacedim>& dof,
 +		            const typename FunctionMap<dim>::type& function_map,
 +		            const std::vector<bool>& component_mask = std::vector<bool>());
 +
                                       /**
                                        * Reset the data structures.
                                        */
diff --cc deal.II/source/dofs/dof_handler.cc
index 2bbf7f7c5b,1d4cacf434..ca50bf37d3
--- a/deal.II/source/dofs/dof_handler.cc
+++ b/deal.II/source/dofs/dof_handler.cc
@@@ -1137,15 -714,15 +1044,15 @@@ DoFHandler<dim,spacedim>::last_active_f
    switch (dim)
      {
        case 1:
 -            Assert (false, ExcImpossibleInDim(1));
 -            return raw_face_iterator();
 +	        Assert (false, ExcImpossibleInDim(1));
 +	        return raw_face_iterator();
        case 2:
 -            return last_active_line ();
 +	       return last_active_line ();
        case 3:
 -            return last_active_quad ();
 +	       return last_active_quad ();
        default:
- 	        Assert (false, ExcNotImplemented());
- 	        return raw_face_iterator ();
+             Assert (false, ExcNotImplemented());
+             return face_iterator ();
      }
  }
  
@@@ -1155,30 -741,34 +1071,32 @@@ DoFHandler<dim,spacedim>::begin_raw_lin
  
  
  template <int dim, int spacedim>
- typename DoFHandler<dim, spacedim>::raw_line_iterator
- DoFHandler<dim, spacedim>::begin_raw_line (const unsigned int level) const
+ typename DoFHandler<dim,spacedim>::line_iterator
+ DoFHandler<dim,spacedim>::begin_line (const unsigned int level) const
  {
-   Assert(tria != 0, ExcNotInitialized());
-   switch (dim)
-     {
-       case 1:
- 	        Assert (level<tria->n_levels(), ExcInvalidLevel(level));
+   typename Triangulation<dim,spacedim>::line_iterator ti = tria->begin_line(level);
+   typename DoFHandler<dim, spacedim>::line_iterator ri(*ti, this);
+   return ri;
+ }
  
- 	        if (tria->n_raw_lines(level) == 0)
- 	          return end_line ();
 -
 -
+ template <int dim, int spacedim>
+ typename DoFHandler<dim, spacedim>::active_line_iterator
+ DoFHandler<dim, spacedim>::begin_active_line (const unsigned int level) const
+ {
+   typename Triangulation<dim,spacedim>::active_line_iterator ti = tria->begin_active_line(level);
+   typename DoFHandler<dim, spacedim>::active_line_iterator ri(*ti, this);
+   return ri;
+ }
  
- 	        return raw_line_iterator (tria,
- 				                      level,
- 				                      0,
- 				                      this);
  
-       default:
- 	        Assert (level == 0, ExcFacesHaveNoLevel());
- 	        return raw_line_iterator (tria,
- 				                      0,
- 				                      0,
- 				                      this);
-     }
+ 
+ template <int dim, int spacedim>
+ typename DoFHandler<dim, spacedim>::raw_line_iterator
+ DoFHandler<dim, spacedim>::end_raw_line (const unsigned int level) const
+ {
+   typename Triangulation<dim,spacedim>::raw_line_iterator ti = tria->end_raw_line(level);
+   typename DoFHandler<dim, spacedim>::raw_line_iterator ri(ti.access_any(), this);
+   return ri;
  }
  
  
diff --cc deal.II/source/dofs/dof_tools.cc
index a10b59bcf2,a955d6d26b..20805db545
--- a/deal.II/source/dofs/dof_tools.cc
+++ b/deal.II/source/dofs/dof_tools.cc
@@@ -3361,16 -3344,19 +3344,19 @@@ namespace DoFTool
      face_1->get_dof_indices(dofs_1, face_1_index);
      face_2->get_dof_indices(dofs_2, face_2_index);
  
 -                                   // .. and constrain them (respecting
 -                                   // component_mask):
 +                                     // .. and constrain them (respecting
 +                                     // component_mask):
      for (unsigned int i = 0; i < dofs_per_face; ++i) {
-       if (component_mask.size() == 0 ||
-           component_mask[fe.face_system_to_component_index(i).first]) {
-         if(!constraint_matrix.is_constrained(dofs_1[i])) {
-           constraint_matrix.add_line(dofs_1[i]);
-           constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+       if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components())
+           ||
+           (component_mask[fe.face_system_to_component_index(i).first] == true))
+         {
+           if(!constraint_matrix.is_constrained(dofs_1[i]))
+             {
+             constraint_matrix.add_line(dofs_1[i]);
+             constraint_matrix.add_entry(dofs_1[i], dofs_2[i], 1.0);
+             }
          }
-       }
      }
    }
  
@@@ -3399,9 -3385,9 +3385,9 @@@
                                  const types::boundary_id     boundary_component,
                                  const int                      direction,
                                  dealii::Tensor<1,DH::space_dimension>
 -                                                               &offset,
 +                                &offset,
                                  dealii::ConstraintMatrix       &constraint_matrix,
-                                 const std::vector<bool>        &component_mask)
+                                 const ComponentMask            &component_mask)
    {
      static const int space_dim = DH::space_dimension;
      Assert (0<=direction && direction<space_dim,
@@@ -3799,28 -3842,44 +3842,44 @@@
  
  
  
+   template <int dim, int spacedim>
+   void
+   extract_dofs (const DoFHandler<dim,spacedim>   &dof,
+                 const BlockMask     &block_mask,
+                 std::vector<bool>       &selected_dofs)
+   {
+     // simply forward to the function that works based on a component mask
+     extract_dofs (dof, dof.get_fe().component_mask (block_mask),
+                   selected_dofs);
+   }
+ 
+ 
+ 
+   template <int dim, int spacedim>
+   void
+   extract_dofs (const hp::DoFHandler<dim,spacedim>   &dof,
+                 const BlockMask     &block_mask,
+                 std::vector<bool>       &selected_dofs)
+   {
+     // simply forward to the function that works based on a component mask
+     extract_dofs (dof, dof.get_fe().component_mask (block_mask),
+                   selected_dofs);
+   }
+ 
+ 
+ 
 -  template<int dim, int spacedim>
 +  template<class DH>
    void
-   extract_level_dofs(
-     const unsigned int       level,
-     const DH &dof,
-     const std::vector<bool> &component_select,
-     std::vector<bool>       &selected_dofs,
-     const bool               count_by_blocks)
+   extract_level_dofs(const unsigned int       level,
 -                     const MGDoFHandler<dim,spacedim> &dof,
++                     const DH &dof,
+                      const ComponentMask &component_mask,
+                      std::vector<bool>       &selected_dofs)
    {
 -    const FiniteElement<dim,spacedim>& fe = dof.get_fe();
 +    const FiniteElement<DH::dimension,DH::space_dimension>& fe = dof.get_fe();
  
-     if (count_by_blocks == true)
-       {
-         Assert(component_select.size() == fe.n_blocks(),
-                ExcDimensionMismatch(component_select.size(), fe.n_blocks()));
-       }
-     else
-       {
-         Assert(component_select.size() == fe.n_components(),
-                ExcDimensionMismatch(component_select.size(), fe.n_components()));
-       }
- 
+     Assert(component_mask.represents_n_components(n_components(dof)),
+            ExcMessage ("The given component mask is not sized correctly to represent the "
+                        "components of the given finite element."));
      Assert(selected_dofs.size() == dof.n_dofs(level),
             ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level)));
  
@@@ -3902,8 -3914,8 +3914,8 @@@
                                       // then loop over all cells and do
                                       // work
      std::vector<unsigned int> indices(fe.dofs_per_cell);
 -    typename MGDoFHandler<dim,spacedim>::cell_iterator c;
 -    for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
 +    typename DH::cell_iterator c;
-     for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)
++    for (c = dof.begin(level) ; c != dof.end(level) ; ++ c)	
        {
          c->get_mg_dof_indices(indices);
          for (unsigned int i=0; i<fe.dofs_per_cell; ++i)
@@@ -3913,6 -3925,20 +3925,20 @@@
  
  
  
 -  template<int dim, int spacedim>
++  template<class DH>
+   void
+   extract_level_dofs(const unsigned int       level,
 -                     const MGDoFHandler<dim,spacedim> &dof,
++                     const DH &dof,
+                      const BlockMask         &block_mask,
+                      std::vector<bool>       &selected_dofs)
+   {
+     // simply defer to the other extract_level_dofs() function
+     extract_level_dofs (level, dof, dof.get_fe().component_mask(block_mask),
+                         selected_dofs);
+   }
+ 
+ 
+ 
    template <class DH>
    void
    extract_boundary_dofs (const DH                      &dof_handler,
@@@ -3928,10 -3954,10 +3954,10 @@@
                          "See the documentation for more information."));
  
      IndexSet indices;
-     extract_boundary_dofs (dof_handler, component_select,
+     extract_boundary_dofs (dof_handler, component_mask,
                             indices, boundary_indicators);
  
 -    // clear and reset array by default values
 +                                     // clear and reset array by default values
      selected_dofs.clear ();
      selected_dofs.resize (dof_handler.n_dofs(), false);
  
@@@ -4008,44 -4037,44 +4037,44 @@@
                    if (!check_vector_component)
                      selected_dofs.add_index (face_dof_indices[i]);
                    else
 -                    // check for
 -                    // component is
 -                    // required. somewhat
 -                    // tricky as usual
 -                    // for the case that
 -                    // the shape function
 -                    // is non-primitive,
 -                    // but use usual
 -                    // convention (see
 -                    // docs)
 +                                                     // check for
 +                                                     // component is
 +                                                     // required. somewhat
 +                                                     // tricky as usual
 +                                                     // for the case that
 +                                                     // the shape function
 +                                                     // is non-primitive,
 +                                                     // but use usual
 +                                                     // convention (see
 +                                                     // docs)
                      {
 -                      // first get at the
 -                      // cell-global
 -                      // number of a face
 -                      // dof, to ask the
 -                      // fe certain
 -                      // questions
 +                                                       // first get at the
 +                                                       // cell-global
 +                                                       // number of a face
 +                                                       // dof, to ask the
 +                                                       // fe certain
 +                                                       // questions
                        const unsigned int cell_index
 -                      = (dim == 1 ?
 -                          i
 -                          :
 -                          (dim == 2 ?
 -                              (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
 +                        = (dim == 1 ?
 +                           i
 +                           :
 +                           (dim == 2 ?
 +                            (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
 +                            :
 +                            (dim == 3 ?
 +                             (i<4*fe.dofs_per_vertex ?
 +                              i
                                :
 -                              (dim == 3 ?
 -                                  (i<4*fe.dofs_per_vertex ?
 -                                      i
 -                                      :
 -                                      (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
 -                                          i+4*fe.dofs_per_vertex
 -                                          :
 -                                          i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
 -                                          :
 -                                          numbers::invalid_unsigned_int)));
 +                              (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
 +                               i+4*fe.dofs_per_vertex
 +                               :
 +                               i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
 +                             :
 +                             numbers::invalid_unsigned_int)));
                        if (fe.is_primitive (cell_index))
                          {
-                           if (component_select[fe.face_system_to_component_index(i).first]
+                           if (component_mask[fe.face_system_to_component_index(i).first]
 -                                               == true)
 +                              == true)
                              selected_dofs.add_index (face_dof_indices[i]);
                          }
                        else // not primitive
diff --cc deal.II/source/dofs/dof_tools.inst.in
index 206641e784,b0f889f3f6..68536b3d24
--- a/deal.II/source/dofs/dof_tools.inst.in
+++ b/deal.II/source/dofs/dof_tools.inst.in
@@@ -425,15 -431,25 +431,37 @@@ template void DoFTools::extract_dofs<de
  
  template void DoFTools::extract_dofs<deal_II_dimension>
  (const hp::DoFHandler<deal_II_dimension>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+  const ComponentMask &,
+  std::vector<bool>&);
+ 
+ template void DoFTools::extract_dofs<deal_II_dimension>
+ (const hp::DoFHandler<deal_II_dimension>&,
+  const BlockMask &,
+  std::vector<bool>&);
+ 
 -template void DoFTools::extract_level_dofs<deal_II_dimension>
++template void DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension> >
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension>&,
++ const ComponentMask &,
++ std::vector<bool>&);
 +
 +template void DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension> >
- (const unsigned int, const DoFHandler<deal_II_dimension>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
++(const unsigned int level,
++ const DoFHandler<deal_II_dimension>&,
++ const BlockMask &,
++ std::vector<bool>&);
 +
 +template void DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension> >
- (const unsigned int, const MGDoFHandler<deal_II_dimension>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+  const MGDoFHandler<deal_II_dimension>&,
+  const ComponentMask &,
+  std::vector<bool>&);
+ 
 -template void DoFTools::extract_level_dofs<deal_II_dimension>
++template void DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension> >
+ (const unsigned int level,
+  const MGDoFHandler<deal_II_dimension>&,
+  const BlockMask &,
+  std::vector<bool>&);
  
  template
  void
@@@ -734,15 -750,19 +762,19 @@@ DoFTools::count_dofs_per_component<hp::
  
  template
  void
 -DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
 +DoFTools::extract_level_dofs<DoFHandler<deal_II_dimension, deal_II_dimension+1> >
- (const unsigned int, const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
 - const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
++ const DoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+  const ComponentMask&,
+  std::vector<bool>&);
  
  template
  void
 -DoFTools::extract_level_dofs<deal_II_dimension, deal_II_dimension+1>
 +DoFTools::extract_level_dofs<MGDoFHandler<deal_II_dimension, deal_II_dimension+1> >
- (const unsigned int, const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+  const MGDoFHandler<deal_II_dimension, deal_II_dimension+1>&,
+  const BlockMask&,
+  std::vector<bool>&);
  
  #endif
  
@@@ -762,15 -782,19 +794,19 @@@ DoFTools::count_dofs_per_component<hp::
  
  template
  void
 -DoFTools::extract_level_dofs<1,3>
 +DoFTools::extract_level_dofs<DoFHandler<1, 3> >
- (const unsigned int, const DoFHandler<1, 3>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
 - const MGDoFHandler<1,3>&,
++ const DoFHandler<1,3>&,
+  const ComponentMask &,
+  std::vector<bool>&);
  
  template
  void
 -DoFTools::extract_level_dofs<1,3>
 +DoFTools::extract_level_dofs<MGDoFHandler<1, 3> >
- (const unsigned int, const MGDoFHandler<1, 3>&,
-  const std::vector<bool>&, std::vector<bool>&, bool);
+ (const unsigned int level,
+  const MGDoFHandler<1,3>&,
+  const BlockMask &,
+  std::vector<bool>&);
  
  #endif
  
diff --cc deal.II/source/multigrid/mg_tools.cc
index c81d632032,15d96c3c86..7d2e640eba
--- a/deal.II/source/multigrid/mg_tools.cc
+++ b/deal.II/source/multigrid/mg_tools.cc
@@@ -30,9 -30,9 +30,10 @@@
  #include <deal.II/multigrid/mg_tools.h>
  #include <deal.II/multigrid/mg_base.h>
  #include <deal.II/base/mg_level_object.h>
 +#include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/fe/fe.h>
+ #include <deal.II/fe/component_mask.h>
  
  #include <vector>
  #include <algorithm>
@@@ -977,60 -978,61 +978,61 @@@ namespace MGTool
                                             // otherwise determine the number
                                             // of dofs in each component
                                             // separately. do so in parallel
- 	      std::vector<std::vector<bool> >
- 	        dofs_in_component (n_components,
- 			                   std::vector<bool>(dof_handler.n_dofs(l),
- 						                         false));
- 	      std::vector<std::vector<bool> >
- 	        component_select (n_components,
- 			                  std::vector<bool>(n_components, false));
- 	      Threads::TaskGroup<> tasks;
- 	      for (unsigned int i=0; i<n_components; ++i)
- 	        {
- 	          void (*fun_ptr) (const unsigned int       level,
- 			                   const MGDoFHandler<dim,spacedim>    &,
- 			                   const std::vector<bool>    &,
- 			                   std::vector<bool>          &,
- 			                   bool)
- 		        = &DoFTools::template extract_level_dofs<MGDoFHandler<dim, spacedim> >;
- 	          component_select[i][i] = true;
- 	          tasks += Threads::new_task (fun_ptr,
- 					                      l, dof_handler,
- 					                      component_select[i],
- 					                      dofs_in_component[i], false);
- 	        }
- 	      tasks.join_all();
+           std::vector<std::vector<bool> >
+             dofs_in_component (n_components,
+                                std::vector<bool>(dof_handler.n_dofs(l),
+                                                  false));
+           std::vector<ComponentMask> component_select (n_components);
+           Threads::TaskGroup<> tasks;
+           for (unsigned int i=0; i<n_components; ++i)
+             {
+               void (*fun_ptr) (const unsigned int       level,
+                                const MGDoFHandler<dim,spacedim>    &,
+                                const ComponentMask    &,
+                                std::vector<bool>          &)
 -                = &DoFTools::template extract_level_dofs<dim>;
++                = &DoFTools::template extract_level_dofs<MGDoFHandler<dim,spacedim> >;
+ 
+ 	      std::vector<bool> tmp(n_components, false);
+ 	      tmp[i] = true;
+ 	      component_select[i] = ComponentMask(tmp);
+ 
+               tasks += Threads::new_task (fun_ptr,
+                                           l, dof_handler,
+                                           component_select[i],
+                                           dofs_in_component[i]);
+             }
+           tasks.join_all();
  
                                             // next count what we got
 -          unsigned int component = 0;
 -          for (unsigned int b=0;b<fe.n_base_elements();++b)
 -            {
 -              const FiniteElement<dim>& base = fe.base_element(b);
 +	      unsigned int component = 0;
 +	      for (unsigned int b=0;b<fe.n_base_elements();++b)
 +	        {
 +	          const FiniteElement<dim>& base = fe.base_element(b);
                                                 // Dimension of base element
 -              unsigned int d = base.n_components();
 -
 -              for (unsigned int m=0;m<fe.element_multiplicity(b);++m)
 -                {
 -                  for (unsigned int dd=0;dd<d;++dd)
 -                    {
 -                      if (base.is_primitive() || (!only_once || dd==0))
 -                        result[l][target_component[component]]
 -                          += std::count(dofs_in_component[component].begin(),
 -                                        dofs_in_component[component].end(),
 -                                        true);
 -                      ++component;
 -                    }
 -                }
 -            }
 +	          unsigned int d = base.n_components();
 +
 +	          for (unsigned int m=0;m<fe.element_multiplicity(b);++m)
 +	  	        {
 +		          for (unsigned int dd=0;dd<d;++dd)
 +		            {
 +		              if (base.is_primitive() || (!only_once || dd==0))
 +			            result[l][target_component[component]]
 +			              += std::count(dofs_in_component[component].begin(),
 +					                    dofs_in_component[component].end(),
 +					                    true);
 +		              ++component;
 +		            }
 +		        }
 +	        }
                                             // finally sanity check
 -          Assert (!dof_handler.get_fe().is_primitive()
 -                  ||
 -                  std::accumulate (result[l].begin(),
 -                                   result[l].end(), 0U)
 -                  ==
 -                  dof_handler.n_dofs(l),
 -                  ExcInternalError());
 -        }
 +	      Assert (!dof_handler.get_fe().is_primitive()
 +		          ||
 +		          std::accumulate (result[l].begin(),
 +				                   result[l].end(), 0U)
 +		          ==
 +		          dof_handler.n_dofs(l),
 +		          ExcInternalError());
 +	    }
        }
    }
  
@@@ -1099,23 -1101,25 +1101,25 @@@
      for (unsigned int l=0;l<n_levels;++l)
        {
          std::vector<std::vector<bool> >
 -          dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
 +	      dofs_in_block (n_blocks, std::vector<bool>(dof_handler.n_dofs(l), false));
-         std::vector<std::vector<bool> >
- 	      block_select (n_blocks, std::vector<bool>(n_blocks, false));
+         std::vector<BlockMask> block_select (n_blocks);
          Threads::TaskGroup<> tasks;
          for (unsigned int i=0; i<n_blocks; ++i)
- 	      {
- 	        void (*fun_ptr) (const unsigned int level,
- 			                 const DH&,
- 			                 const std::vector<bool>&,
- 			                 std::vector<bool>&,
- 			                 bool)
- 	          = &DoFTools::template extract_level_dofs<DH>;
- 	        block_select[i][i] = true;
- 	        tasks += Threads::new_task (fun_ptr,
- 				                        l, dof_handler, block_select[i],
- 				                        dofs_in_block[i], true);
- 	      };
+           {
+             void (*fun_ptr) (const unsigned int level,
 -                             const MGDoFHandler<dim,spacedim>&,
++                             const DH&,
+                              const BlockMask &,
+                              std::vector<bool>&)
 -              = &DoFTools::template extract_level_dofs<dim>;
++              = &DoFTools::template extract_level_dofs<DH>;
+ 
+             std::vector<bool> tmp(n_blocks, false);
+ 	    tmp[i] = true;
+ 	    block_select[i] = tmp;
+ 
+             tasks += Threads::new_task (fun_ptr,
+                                         l, dof_handler, block_select[i],
+                                         dofs_in_block[i]);
+           }
          tasks.join_all ();
  
                                           // next count what we got
@@@ -1179,58 -1184,59 +1183,59 @@@
  
      std::vector<unsigned int> local_dofs;
      local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
-     std::fill (local_dofs.begin (), local_dofs.end (),
+     std::fill (local_dofs.begin (),
+     		   local_dofs.end (),
 -               DoFHandler<dim,spacedim>::invalid_dof_index);
 +	           DoFHandler<dim,spacedim>::invalid_dof_index);
  
                                       // First, deal with the simpler
                                       // case when we have to identify
                                       // all boundary dofs
-     if (component_mask.size() == 0)
+     if (component_mask.n_selected_components(n_components) == n_components)
        {
          typename MGDoFHandler<dim,spacedim>::cell_iterator
 -          cell = dof.begin(),
 -          endc = dof.end();
 +  	      cell = dof.begin(),
 +	      endc = dof.end();
          for (; cell!=endc; ++cell)
 -          {
 -            const FiniteElement<dim> &fe = cell->get_fe();
 -            const unsigned int level = cell->level();
 -            local_dofs.resize(fe.dofs_per_face);
 -
 -            for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 -                 ++face_no)
 -              if (cell->at_boundary(face_no) == true)
 -                {
 -                  const typename MGDoFHandler<dim,spacedim>::face_iterator
 -                    face = cell->face(face_no);
 -                  const types::boundary_id bi = face->boundary_indicator();
 +	      {
 +	        const FiniteElement<dim> &fe = cell->get_fe();
 +	        const unsigned int level = cell->level();
 +	        local_dofs.resize(fe.dofs_per_face);
 +
 +	        for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 +	             ++face_no)
 +	          if (cell->at_boundary(face_no) == true)
 +	            {
 +		          const typename MGDoFHandler<dim,spacedim>::face_iterator
 +		            face = cell->face(face_no);
 +		          const types::boundary_id bi = face->boundary_indicator();
                                                     // Face is listed in
                                                     // boundary map
 -                  if (function_map.find(bi) != function_map.end())
 -                    {
 -                      face->get_mg_dof_indices(level, local_dofs);
 -                      for (unsigned int i=0;i<fe.dofs_per_face;++i)
 -                        boundary_indices[level].insert(local_dofs[i]);
 -                    }
 -                }
 -          }
 +		          if (function_map.find(bi) != function_map.end())
 +		            {
 +		              face->get_mg_dof_indices(level, local_dofs);
 +		              for (unsigned int i=0;i<fe.dofs_per_face;++i)
 +		                boundary_indices[level].insert(local_dofs[i]);
 +		            }
 +	            }
 +	      }
        }
      else
        {
-         Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
+         Assert (component_mask.n_selected_components(n_components) > 0,
 -                ExcMessage("It's probably worthwhile to select at least one component."));
 +	            ExcMessage("It's probably worthwhile to select at least one component."));
  
          typename MGDoFHandler<dim,spacedim>::cell_iterator
 -          cell = dof.begin(),
 -          endc = dof.end();
 +	      cell = dof.begin(),
 +	      endc = dof.end();
          for (; cell!=endc; ++cell)
 -          for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 -               ++face_no)
 -            {
 -              if (cell->at_boundary(face_no) == false)
 -                continue;
 +	      for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 +	           ++face_no)
 +	        {
 +	          if (cell->at_boundary(face_no) == false)
 +	            continue;
  
 -              const FiniteElement<dim> &fe = cell->get_fe();
 -              const unsigned int level = cell->level();
 +	          const FiniteElement<dim> &fe = cell->get_fe();
 +	          const unsigned int level = cell->level();
  
                                                 // we can presently deal only with
                                                 // primitive elements for boundary
@@@ -1242,25 -1248,25 +1247,25 @@@
                                                 // that are non-zero for the
                                                 // components we are interested in,
                                                 // are in fact primitive
- 	          for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
- 	            {
- 		          const std::vector<bool> &nonzero_component_array
- 		            = cell->get_fe().get_nonzero_components (i);
- 		          for (unsigned int c=0; c<n_components; ++c)
- 		            if ((nonzero_component_array[c] == true)
- 		                &&
- 		                (component_mask[c] == true))
- 		              Assert (cell->get_fe().is_primitive (i),
- 			                  ExcMessage ("This function can only deal with requested boundary "
- 					                      "values that correspond to primitive (scalar) base "
- 					                      "elements"));
- 	            }
+               for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+                 {
+                   const ComponentMask &nonzero_component_array
+                     = cell->get_fe().get_nonzero_components (i);
+                   for (unsigned int c=0; c<n_components; ++c)
+                     if ((nonzero_component_array[c] == true)
+                         &&
+                         (component_mask[c] == true))
+                       Assert (cell->get_fe().is_primitive (i),
+                               ExcMessage ("This function can only deal with requested boundary "
+                                           "values that correspond to primitive (scalar) base "
+                                           "elements"));
+                 }
  
 -              typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
 -              const types::boundary_id boundary_component = face->boundary_indicator();
 -              if (function_map.find(boundary_component) != function_map.end())
 +	          typename MGDoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
 +	          const types::boundary_id boundary_component = face->boundary_indicator();
 +	          if (function_map.find(boundary_component) != function_map.end())
                                                   // face is of the right component
 -                {
 +	            {
                                                     // get indices, physical location and
                                                     // boundary values of dofs on this
                                                     // face
@@@ -1347,22 -1354,18 +1352,18 @@@
  // components. if shape function is non-primitive, then we will ignore
  // the result in the following anyway, otherwise there's only one
  // non-zero component which we will use
- 			                  component = (std::find (fe.get_nonzero_components(cell_i).begin(),
- 						                              fe.get_nonzero_components(cell_i).end(),
- 						                              true)
- 					                       -
- 					                       fe.get_nonzero_components(cell_i).begin());
- 			                }
+                               component = fe.get_nonzero_components(cell_i).first_selected_component();
+                             }
  
 -                          if (component_mask[component] == true)
 -                            boundary_indices[level].insert(local_dofs[i]);
 -                        }
 -                    }
 -                  else
 -                    for (unsigned int i=0; i<local_dofs.size(); ++i)
 -                      boundary_indices[level].insert(local_dofs[i]);
 -                }
 -            }
 +		  	              if (component_mask[component] == true)
 +			                boundary_indices[level].insert(local_dofs[i]);
 +		                }
 +		            }
 +		          else
 +		            for (unsigned int i=0; i<local_dofs.size(); ++i)
 +		              boundary_indices[level].insert(local_dofs[i]);
 +	            }
 +	        }
        }
    }
  
@@@ -1370,229 -1373,15 +1371,225 @@@
  
    template <int dim, int spacedim>
    void
 +  make_boundary_list(
 +    const DoFHandler<dim,spacedim>& dof,
 +    const typename FunctionMap<dim>::type& function_map,
 +    std::vector<std::set<unsigned int> >& boundary_indices,
 +    const std::vector<bool>& component_mask)
 +  {
 +                                   // if for whatever reason we were
 +                                   // passed an empty map, return
 +                                   // immediately
 +    if (function_map.size() == 0)
 +      return;
 +
 +    const unsigned int n_levels = dof.get_tria().n_levels();
 +
 +
  
 +    const unsigned int n_components = DoFTools::n_components(dof);
 +    const bool          fe_is_system = (n_components != 1);
 +
 +    AssertDimension (boundary_indices.size(), n_levels);
 +
 +    std::vector<unsigned int> local_dofs;
 +    local_dofs.reserve (DoFTools::max_dofs_per_face(dof));
 +    std::fill (local_dofs.begin (), local_dofs.end (),
 +	       DoFHandler<dim,spacedim>::invalid_dof_index);
 +
 +				   // First, deal with the simpler
 +				   // case when we have to identify
 +				   // all boundary dofs
 +    if (component_mask.size() == 0)
 +      {
 +        typename DoFHandler<dim,spacedim>::cell_iterator
 +	  cell = dof.begin(),
 +	  endc = dof.end();
 +        for (; cell!=endc; ++cell)
 +	  {
 +	    const FiniteElement<dim> &fe = cell->get_fe();
 +	    const unsigned int level = cell->level();
 +	    local_dofs.resize(fe.dofs_per_face);
 +
 +	    for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 +	         ++face_no)
 +	      if (cell->at_boundary(face_no))
 +	        {
 +		  const typename DoFHandler<dim,spacedim>::face_iterator
 +		    face = cell->face(face_no);
 +		  const unsigned char bi = face->boundary_indicator();
 +						 // Face is listed in
 +						 // boundary map
 +		  if (function_map.find(bi) != function_map.end())
 +		    {
 +		      face->get_mg_dof_indices(level, local_dofs);
 +		      for (unsigned int i=0;i<fe.dofs_per_face;++i)
 +		        boundary_indices[level].insert(local_dofs[i]);
 +		    }
 +	        }
 +	  }
 +      }
 +    else
 +      {
 +        Assert (std::count(component_mask.begin(), component_mask.end(), true) > 0,
 +	        ExcMessage("It's probably worthwhile to select at least one component."));
 +
 +        typename DoFHandler<dim,spacedim>::cell_iterator
 +	  cell = dof.begin(),
 +	  endc = dof.end();
 +        for (; cell!=endc; ++cell)
 +	  for (unsigned int face_no = 0; face_no < GeometryInfo<dim>::faces_per_cell;
 +	       ++face_no)
 +	    {
 +	      if (!(cell->at_boundary(face_no)))
 +	        continue;
 +
 +	      const FiniteElement<dim> &fe = cell->get_fe();
 +	      const unsigned int level = cell->level();
 +
 +					     // we can presently deal only with
 +					     // primitive elements for boundary
 +					     // values. this does not preclude
 +					     // us using non-primitive elements
 +					     // in components that we aren't
 +					     // interested in, however. make
 +					     // sure that all shape functions
 +					     // that are non-zero for the
 +					     // components we are interested in,
 +					     // are in fact primitive
 +	      for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
 +	        {
- 		  const std::vector<bool> &nonzero_component_array
- 		    = cell->get_fe().get_nonzero_components (i);
++                  const ComponentMask &nonzero_component_array
++                    = cell->get_fe().get_nonzero_components (i);
 +		  for (unsigned int c=0; c<n_components; ++c)
 +		    if ((nonzero_component_array[c] == true)
 +		        &&
 +		        (component_mask[c] == true))
 +		      Assert (cell->get_fe().is_primitive (i),
 +			      ExcMessage ("This function can only deal with requested boundary "
 +				 	  "values that correspond to primitive (scalar) base "
 +					  "elements"));
 +	        }
 +
 +	      typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_no);
 +	      const unsigned char boundary_component = face->boundary_indicator();
 +	      if (function_map.find(boundary_component) != function_map.end())
 +					       // face is of the right component
 +	        {
 +						 // get indices, physical location and
 +						 // boundary values of dofs on this
 +						 // face
 +		  local_dofs.resize (fe.dofs_per_face);
 +		  face->get_mg_dof_indices (level, local_dofs);
 +		  if (fe_is_system)
 +		    {
 +						     // enter those dofs
 +						     // into the list that
 +						     // match the
 +						     // component
 +						     // signature. avoid
 +						     // the usual
 +						     // complication that
 +						     // we can't just use
 +						     // *_system_to_component_index
 +						     // for non-primitive
 +						     // FEs
 +		      for (unsigned int i=0; i<local_dofs.size(); ++i)
 +		        {
 +			  unsigned int component;
 +			  if (fe.is_primitive())
 +			    component = fe.face_system_to_component_index(i).first;
 +			  else
 +			    {
 +							     // non-primitive
 +							     // case. make
 +							     // sure that
 +							     // this
 +							     // particular
 +							     // shape
 +							     // function
 +							     // _is_
 +							     // primitive,
 +							     // and get at
 +							     // it's
 +							     // component. use
 +							     // usual
 +							     // trick to
 +							     // transfer
 +							     // face dof
 +							     // index to
 +							     // cell dof
 +
 +							     // index
 +			      const unsigned int cell_i
 +			        = (dim == 1 ?
 +				   i
 +				   :
 +				   (dim == 2 ?
 +				    (i<2*fe.dofs_per_vertex ? i : i+2*fe.dofs_per_vertex)
 +				    :
 +				    (dim == 3 ?
 +				     (i<4*fe.dofs_per_vertex ?
 +				      i
 +				      :
 +				      (i<4*fe.dofs_per_vertex+4*fe.dofs_per_line ?
 +				       i+4*fe.dofs_per_vertex
 +				       :
 +				       i+4*fe.dofs_per_vertex+8*fe.dofs_per_line))
 +				     :
 +				     numbers::invalid_unsigned_int)));
 +			      Assert (cell_i < fe.dofs_per_cell, ExcInternalError());
 +
 +							     // make sure
 +							     // that if
 +							     // this is
 +							     // not a
 +							     // primitive
 +							     // shape function,
 +							     // then all
 +							     // the
 +							     // corresponding
 +							     // components
 +							     // in the
 +							     // mask are
 +							     // not set
 +//                         if (!fe.is_primitive(cell_i))
 +//                           for (unsigned int c=0; c<n_components; ++c)
 +//                             if (fe.get_nonzero_components(cell_i)[c])
 +//                               Assert (component_mask[c] == false,
 +//                                       ExcFENotPrimitive());
 +
 +// let's pick the first of possibly more than one non-zero
 +// components. if shape function is non-primitive, then we will ignore
 +// the result in the following anyway, otherwise there's only one
 +// non-zero component which we will use
- 			      component = (std::find (fe.get_nonzero_components(cell_i).begin(),
- 						      fe.get_nonzero_components(cell_i).end(),
- 						      true)
- 					   -
- 					   fe.get_nonzero_components(cell_i).begin());
++                              component = fe.get_nonzero_components(cell_i).first_selected_component();
 +			    }
 +
 +			  if (component_mask[component] == true)
 +			    boundary_indices[level].insert(local_dofs[i]);
 +		        }
 +		    }
 +		  else
 +		    for (unsigned int i=0; i<local_dofs.size(); ++i)
 +		      boundary_indices[level].insert(local_dofs[i]);
 +	        }
 +	    }
 +      }
 +  }
 +
 +
 +
 +  template <int dim, int spacedim>
 +  void
    make_boundary_list(const MGDoFHandler<dim,spacedim>& dof,
- 		             const typename FunctionMap<dim>::type& function_map,
- 		             std::vector<IndexSet>& boundary_indices,
- 		             const std::vector<bool>& component_mask)
+                      const typename FunctionMap<dim>::type& function_map,
+                      std::vector<IndexSet>& boundary_indices,
+                      const ComponentMask & component_mask)
    {
      Assert (boundary_indices.size() == dof.get_tria().n_levels(),
 -            ExcDimensionMismatch (boundary_indices.size(),
 -                                  dof.get_tria().n_levels()));
 +	        ExcDimensionMismatch (boundary_indices.size(),
 +				                  dof.get_tria().n_levels()));
  
      std::vector<std::set<unsigned int> >
        my_boundary_indices (dof.get_tria().n_levels());
diff --cc deal.II/source/multigrid/mg_tools.inst.in
index 4f657ada47,9ec731c1e7..24932ea63e
--- a/deal.II/source/multigrid/mg_tools.inst.in
+++ b/deal.II/source/multigrid/mg_tools.inst.in
@@@ -101,16 -92,10 +101,16 @@@ for (deal_II_dimension : DIMENSIONS
  	const MGDoFHandler<deal_II_dimension>&,
  	const FunctionMap<deal_II_dimension>::type&,
  	std::vector<std::set<unsigned int> >&,
- 	const std::vector<bool>&);
+ 	const ComponentMask &);
 +
 +      template void make_boundary_list(
 +	const DoFHandler<deal_II_dimension>&,
 +	const FunctionMap<deal_II_dimension>::type&,
 +	std::vector<std::set<unsigned int> >&,
 +	const std::vector<bool>&);
  #endif
  
-       template void make_boundary_list(
+     template void make_boundary_list(
  	const MGDoFHandler<deal_II_dimension>&,
  	const FunctionMap<deal_II_dimension>::type&,
  	std::vector<IndexSet>&,
@@@ -121,12 -106,7 +121,12 @@@
  	extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
  				      std::vector<std::vector<bool> >  &interface_dofs,
  				      std::vector<std::vector<bool> >  &boundary_interface_dofs);
-       template
+     template
 +	void
 +	extract_inner_interface_dofs (const DoFHandler<deal_II_dimension> &mg_dof_handler,
 +				      std::vector<std::vector<bool> >  &interface_dofs,
 +				      std::vector<std::vector<bool> >  &boundary_interface_dofs);
 +      template
  	void
  	extract_inner_interface_dofs (const MGDoFHandler<deal_II_dimension> &mg_dof_handler,
  				      std::vector<std::vector<bool> >  &interface_dofs);