]> https://gitweb.dealii.org/ - dealii.git/commitdiff
selection with block grouping
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 8 Aug 2003 19:41:54 +0000 (19:41 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Fri, 8 Aug 2003 19:41:54 +0000 (19:41 +0000)
git-svn-id: https://svn.dealii.org/trunk@7907 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/mg_dof_tools.h
deal.II/deal.II/include/multigrid/mg_transfer.h
deal.II/deal.II/include/multigrid/mg_transfer.templates.h
deal.II/deal.II/source/multigrid/mg_dof_tools.cc
deal.II/deal.II/source/multigrid/mg_transfer_block.cc
deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
tests/multigrid/transfer.cc

index 2be97591b50cd7f084a7868e2b07355c1f71c5fd..bc0859e92fbac33b2012b03998388ac2223380a5 100644 (file)
@@ -110,6 +110,30 @@ class MGCoarseGridBase : public Subscriptor
  * context. This class is an abstract one and has no implementations of
  * possible algorithms for these operations.
  *
+ * There are several derived classes, reflecting the fact that vector
+ * types and numbering of the fine-grid discretization and of the
+ * multi-level implementation are independent.
+ *
+ * If you use multigrid for your complete system of equations, you
+ * will use @ref{MGTransferPrebuilt} together with Multigrid for
+ * Vector objects. The fine-grid vectors may still be of type
+ * BlockVector.
+ *
+ * For mixed systems, it may be required to do multigrid only for a
+ * single component or for some components. The classes
+ * @ref{MGTransferSelect} and @ref{MGTransferBlock} handle these cases.
+ *
+ * MGTransferSelect is used if you use mutligrid (on Vector oblects)
+ * for a single component, possibly grouped using
+ * @p{mg_target_component}.
+ *
+ * The class MGTransferBlock handles the case where your multigrid
+ * method operates on BlockVector objects. These can contain all or a
+ * consecutive set of the blocks of the complete system. Since most
+ * smoothers cannot operate on block structures, it is not clear
+ * whether this case is really useful. Therefore, a tested
+ * implementation of this case will be supplied when needed.
+ *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2002
  */
  template <class VECTOR>
index 99ce9dcb268036bb3eaeb15032d4a64f06fceb18..3976262247c1fd0675fa632628ef169772199ca2 100644 (file)
@@ -180,6 +180,38 @@ class MGTools
                     MGLevelObject<BlockVector<number> >& v,
                     const std::vector<bool>& selected,
                     const std::vector<unsigned int>& target_component);
+                                    /**
+                                     * Adjust vectors on all levels
+                                     * to correct size.  Count the
+                                     * numbers of degrees of freedom
+                                     * on each level component-wise
+                                     * in a single component. Then,
+                                     * assign @p{vector} the
+                                     * corresponding size.
+                                     *
+                                     * The boolean field @p{selected}
+                                     * may be nonzero in a single
+                                     * component, indicating the
+                                     * block of a block vector the
+                                     * argument @p{v} corresponds to.
+                                     *
+                                     * Degrees of freedom must be
+                                     * sorted by component in order
+                                     * to obtain reasonable results
+                                     * from this function.
+                                     *
+                                     * The argument
+                                     * @p{target_component} allows to
+                                     * re-sort and groupt components
+                                     * as in
+                                     * @p{DoFRenumbering::component_wise}.
+                                     */
+    template <int dim, typename number>
+      static void
+      reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                    MGLevelObject<Vector<number> >& v,
+                    const std::vector<bool>& selected,
+                    const std::vector<unsigned int>& target_component);
 };
 
 
index 3d405e64c36bc8ec6adf3c1122d8140631887597..8a5388648c4396190a036aebe47a9c7a663688bc 100644 (file)
@@ -41,6 +41,9 @@ template <int dim> class MGDoFHandler;
  * once by looping over all cells and storing the result in a matrix for
  * each level, but requires additional memory.
  *
+ * See @ref{MGTransferBase} to find out which of the transfer classes
+ * is best for your needs.
+ *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999-2003
  */
 template <typename number>
@@ -163,11 +166,6 @@ class MGTransferPrebuilt : public MGTransferBase<Vector<number> >
 
   private:
 
-                                  /**
-                                   * Selected component.
-                                   */
-    unsigned int selected;
-
                                   /**
                                    * Sizes of the multi-level vectors.
                                    */
@@ -253,7 +251,7 @@ class MGTransferPrebuilt : public MGTransferBase<Vector<number> >
 
 /**
  * Implementation of matrix generation for @ref{MGTransferBlock} and
- * @p{MGTransferSelected}.
+ * @p{MGTransferSelect}.
  *
  * @author Guido Kanschat, 2001-2003
  */
@@ -291,6 +289,13 @@ class MGTransferBlockBase
                                      */
     std::vector<unsigned int> target_component;
     
+                                    /**
+                                     * Target component if
+                                     * renumbering of level vectors
+                                     * is required.
+                                     */
+    std::vector<unsigned int> mg_target_component;
+    
                                   /**
                                    * Sizes of the multi-level vectors.
                                    */
@@ -343,11 +348,18 @@ class MGTransferBlockBase
 #endif
 };
 
-
+//TODO:[GK] Update this class
 
 /**
  * Implementation of the @p{MGTransferBase} interface for block
- * matrices and block vectors.  In addition to the functionality of
+ * matrices and block vectors.
+ *
+ * Warning! Due to additional requirements on MGTransferSelect, the
+ * implementation in the base class has changed. Therefore, this class
+ * is left in an untested state. If you use it and you encounter
+ * problems, please contact Guido Kanschat.
+ *
+ * In addition to the functionality of
  * @ref{MGTransferPrebuilt}, the operation may be restricted to
  * certain blocks of the vector.
  *
@@ -355,6 +367,9 @@ class MGTransferBlockBase
  * transfer routines may only have as many components as there are
  * @p{true}s in the selected-field.
  *
+ * See @ref{MGTransferBase} to find out which of the transfer classes
+ * is best for your needs.
+ *
  * @author Guido Kanschat, 2001, 2002
  */
 template <typename number>
@@ -391,6 +406,8 @@ class MGTransferBlock : public MGTransferBase<BlockVector<number> >,
     void build_matrices (const MGDoFHandler<dim> &mg_dof,
                         std::vector<bool> selected = std::vector<bool>(),
                         const std::vector<unsigned int>& target_component
+                        = std::vector<unsigned int>(),
+                        const std::vector<unsigned int>& mg_target_component
                         =std::vector<unsigned int>());
 
                                     /**
@@ -534,13 +551,17 @@ class MGTransferBlock : public MGTransferBase<BlockVector<number> >,
 };
 
 
+//TODO:[GK] Update documentation for copy_* functions
 
 /**
  * Implementation of the @p{MGTransferBase} interface for block
- * matrices and simple vectors. This class uses @ref{MGTransferBlock}
+ * matrices and simple vectors. This class uses @ref{MGTransferBlockBase}
  * selecting a single component or grouping several components into a
  * single block. The transfer operators themselves are implemented for
- * simple vectors again.
+ * Vector and BlockVector objects.
+ *
+ * See @ref{MGTransferBase} to find out which of the transfer classes
+ * is best for your needs.
  *
  * @author Guido Kanschat, 2001, 2002, 2003
  */
@@ -558,8 +579,13 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      * Actually build the prolongation
                                      * matrices for grouped components.
                                      *
-                                     * @p{selected} is the number of
-                                     * the component for which the
+                                     * @p{selected} tells the copy
+                                     * functions operating on single
+                                     * vectors, which component this
+                                     * vector belongs to.
+                                     *
+                                     * @p{mg_selected} is the number
+                                     * of the component for which the
                                      * transfer matrices should be
                                      * built.
                                      *
@@ -580,9 +606,18 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
     template <int dim>
     void build_matrices (const MGDoFHandler<dim> &mg_dof,
                         unsigned int selected,
+                        unsigned int mg_selected,
                         const std::vector<unsigned int>& target_component
+                        = std::vector<unsigned int>(),
+                        const std::vector<unsigned int>& mg_target_component
                         = std::vector<unsigned int>());
 
+                                    /**
+                                     * Change selected
+                                     * component. Handle with care!
+                                     */
+    void select (const unsigned int component);
+    
                                     /**
                                      * Prolongate a vector from level
                                      * @p{to_level-1} to level
@@ -648,11 +683,11 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      * on each of the levels
                                      * separately, i.a. an @p{MGVector}.
                                      */
-    template <int dim, class InVector>
+    template <int dim, typename number2>
     void
     copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
                MGLevelObject<Vector<number> > &dst,
-               const InVector                 &src) const;
+               const Vector<number2>          &src) const;
 
                                     /**
                                      * Transfer from multi-level vector to
@@ -666,10 +701,10 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      * constrained degrees of freedom
                                      * are set to zero.
                                      */
-    template <int dim, class OutVector>
+    template <int dim, typename number2>
     void
     copy_from_mg (const MGDoFHandler<dim>              &mg_dof,
-                 OutVector                            &dst,
+                 Vector<number2>                      &dst,
                  const MGLevelObject<Vector<number> > &src) const;
 
                                     /**
@@ -680,18 +715,90 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      * function, but probably not for
                                      * continuous elements.
                                      */
-    template <int dim, class OutVector>
+    template <int dim, typename number2>
     void
     copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof,
-                     OutVector                            &dst,
+                     Vector<number2>                      &dst,
+                     const MGLevelObject<Vector<number> > &src) const;
+
+                                    /**
+                                     * Transfer from a vector on the
+                                     * global grid to vectors defined
+                                     * on each of the levels
+                                     * separately, i.a. an @p{MGVector}.
+                                     */
+    template <int dim, typename number2>
+    void
+    copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
+               MGLevelObject<Vector<number> > &dst,
+               const BlockVector<number2>     &src) const;
+
+                                    /**
+                                     * Transfer from multi-level vector to
+                                     * normal vector.
+                                     *
+                                     * Copies data from active
+                                     * portions of an MGVector into
+                                     * the respective positions of a
+                                     * @p{Vector<number>}. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom
+                                     * are set to zero.
+                                     */
+    template <int dim, typename number2>
+    void
+    copy_from_mg (const MGDoFHandler<dim>              &mg_dof,
+                 BlockVector<number2>                 &dst,
+                 const MGLevelObject<Vector<number> > &src) const;
+
+                                    /**
+                                     * Add a multi-level vector to a
+                                     * normal vector.
+                                     *
+                                     * Works as the previous
+                                     * function, but probably not for
+                                     * continuous elements.
+                                     */
+    template <int dim, typename number2>
+    void
+    copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof,
+                     BlockVector<number2>                 &dst,
                      const MGLevelObject<Vector<number> > &src) const;
 
   private:
-                                     /**
-                                      * Selected component.
-                                      */
-    unsigned int selected;
+                                    /**
+                                     * Transfer from multi-level vector to
+                                     * normal vector.
+                                     *
+                                     * Copies data from active
+                                     * portions of an MGVector into
+                                     * the respective positions of a
+                                     * @p{Vector<number>}. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom
+                                     * are set to zero.
+                                     */
+    template <int dim, class OutVector>
+    void
+    do_copy_from_mg (const MGDoFHandler<dim>              &mg_dof,
+                    OutVector                            &dst,
+                    const MGLevelObject<Vector<number> > &src,
+                    const unsigned int offset) const;
 
+                                    /**
+                                     * Add a multi-level vector to a
+                                     * normal vector.
+                                     *
+                                     * Works as the previous
+                                     * function, but probably not for
+                                     * continuous elements.
+                                     */
+    template <int dim, class OutVector>
+    void
+    do_copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof,
+                        OutVector                            &dst,
+                        const MGLevelObject<Vector<number> > &src,
+                        const unsigned int offset) const;
 
                                     /**
                                      * Implementation of the
@@ -712,10 +819,11 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      */
     template <int dim, class InVector>
     void
-    copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
-               MGLevelObject<Vector<number> > &dst,
-               const InVector                 &src,
-               const is_1d<true>              &) const;
+    do_copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
+                  MGLevelObject<Vector<number> > &dst,
+                  const InVector                 &src,
+                  const unsigned int              offset,
+                  const is_1d<true>              &) const;
 
                                     /**
                                      * Same for all other space
@@ -723,12 +831,25 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
                                      */
     template <int dim, class InVector>
     void
-    copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
-               MGLevelObject<Vector<number> > &dst,
-               const InVector                 &src,
-               const is_1d<false>             &) const;
+    do_copy_to_mg (const MGDoFHandler<dim>        &mg_dof,
+                  MGLevelObject<Vector<number> > &dst,
+                  const InVector                 &src,
+                  const unsigned int              offset,
+                  const is_1d<false>             &) const;
+
+                                     /**
+                                      * Selected component.
+                                      */
+    unsigned int selected_component;
 };
 
+//----------------------------------------------------------------------//
+template <typename number>
+inline void
+MGTransferSelect<number>::select(const unsigned int component)
+{
+  selected_component = component;
+}
 
 
 #endif
index 173ad810273225d7b9c244eb95465cc510937925..53c3b513b5ac2f32af62443d7ece54190bc3a902 100644 (file)
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
-                                       MGLevelObject<Vector<number> > &dst,
-                                       const InVector                 &src) const
+MGTransferPrebuilt<number>::copy_to_mg (
+  const MGDoFHandler<dim>        &mg_dof_handler,
+  MGLevelObject<Vector<number> > &dst,
+  const InVector                 &src) const
 {
                                   // forward to the correct
                                   // specialization
@@ -47,10 +48,11 @@ MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &mg_dof_h
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &,
-                                       MGLevelObject<Vector<number> > &,
-                                       const InVector                 &,
-                                       const is_1d<true>              &) const
+MGTransferPrebuilt<number>::copy_to_mg (
+  const MGDoFHandler<dim>        &,
+  MGLevelObject<Vector<number> > &,
+  const InVector                 &,
+  const is_1d<true>              &) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -59,10 +61,11 @@ MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
-                                       MGLevelObject<Vector<number> > &dst,
-                                       const InVector                 &src,
-                                       const is_1d<false>             &) const
+MGTransferPrebuilt<number>::copy_to_mg (
+  const MGDoFHandler<dim>        &mg_dof_handler,
+  MGLevelObject<Vector<number> > &dst,
+  const InVector                 &src,
+  const is_1d<false>             &) const
 {
                                   // Make src a real finite element function
 //  InVector src = osrc;
@@ -159,10 +162,10 @@ MGTransferPrebuilt<number>::copy_to_mg (const MGDoFHandler<dim>        &mg_dof_h
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferPrebuilt<number>::
-copy_from_mg(const MGDoFHandler<dim>              &mg_dof_handler,
-             OutVector                            &dst,
-             const MGLevelObject<Vector<number> > &src) const
+MGTransferPrebuilt<number>::copy_from_mg(
+  const MGDoFHandler<dim>              &mg_dof_handler,
+  OutVector                            &dst,
+  const MGLevelObject<Vector<number> > &src) const
 {
   const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
 
@@ -207,10 +210,10 @@ copy_from_mg(const MGDoFHandler<dim>              &mg_dof_handler,
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferPrebuilt<number>::
-copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
-                  OutVector                            &dst,
-                  const MGLevelObject<Vector<number> > &src) const
+MGTransferPrebuilt<number>::copy_from_mg_add (
+  const MGDoFHandler<dim>              &mg_dof_handler,
+  OutVector                            &dst,
+  const MGLevelObject<Vector<number> > &src) const
 {
   const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
 
@@ -259,16 +262,33 @@ copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
 
 
 template <typename number>
-template <int dim, class InVector>
+template <int dim, typename number2>
 void
-MGTransferSelect<number>::
-copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
-            MGLevelObject<Vector<number> > &dst,
-            const InVector                 &src) const
+MGTransferSelect<number>::copy_to_mg (
+  const MGDoFHandler<dim>        &mg_dof_handler,
+  MGLevelObject<Vector<number> > &dst,
+  const BlockVector<number2>     &src) const
 {
                                   // forward to the correct
                                   // specialization
-  copy_to_mg (mg_dof_handler, dst, src, is_1d<(dim==1)>());
+  do_copy_to_mg (mg_dof_handler, dst, src, 0, is_1d<(dim==1)>());
+}
+
+
+
+template <typename number>
+template <int dim, typename number2>
+void
+MGTransferSelect<number>::copy_to_mg (
+  const MGDoFHandler<dim>        &mg_dof_handler,
+  MGLevelObject<Vector<number> > &dst,
+  const Vector<number2>          &src) const
+{
+                                  // forward to the correct
+                                  // specialization
+  do_copy_to_mg (mg_dof_handler, dst, src,
+                component_start[selected_component],
+                is_1d<(dim==1)>());
 }
 
 
@@ -276,11 +296,12 @@ copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferSelect<number>::
-copy_to_mg (const MGDoFHandler<dim>        &,
-            MGLevelObject<Vector<number> > &,
-            const InVector                 &,
-           const is_1d<true>              &) const
+MGTransferSelect<number>::do_copy_to_mg (
+  const MGDoFHandler<dim>&,
+  MGLevelObject<Vector<number> >&,
+  const InVector&,
+  const unsigned int,
+  const is_1d<true>&) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -290,14 +311,15 @@ copy_to_mg (const MGDoFHandler<dim>        &,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferSelect<number>::
-copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
-            MGLevelObject<Vector<number> > &dst,
-            const InVector                 &osrc,
-           const is_1d<false>             &) const
+MGTransferSelect<number>::do_copy_to_mg (
+  const MGDoFHandler<dim>&        mg_dof_handler,
+  MGLevelObject<Vector<number> >& dst,
+  const InVector&                 src,
+  const unsigned int              offset,
+  const is_1d<false>&) const
 {
                                   // Make src a real finite element function
-  InVector src = osrc;
+//  InVector src = osrc;
 //  constraints->distribute(src);
 
   const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
@@ -313,17 +335,13 @@ copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
   
   Assert(sizes.size()==mg_dof_handler.get_tria().n_levels(),
         ExcMatricesNotBuilt());
-  for (unsigned int l=minlevel;l<=maxlevel;++l)
-    dst[l].reinit(sizes[l][selected]);
+
+  MGTools::reinit_vector(mg_dof_handler, dst, selected, mg_target_component);
   
   std::vector<unsigned int> global_dof_indices (dofs_per_cell);
   std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
   std::vector<unsigned int> level_face_indices (dofs_per_face);
 
-                                  // Start of input vector inside a
-                                  // block vector.
-  const unsigned int start = component_start[selected];
-
                                   // traverse the grid top-down
                                   // (i.e. starting with the most
                                   // refined grid). this way, we can
@@ -336,10 +354,6 @@ copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
                                   // already have built.
   for (int level=maxlevel; level>=static_cast<signed int>(minlevel); --level)
     {
-                                      // Start of treated component
-                                      // for this level
-      const unsigned int level_start = mg_component_start[level][selected];
-
       typename MGDoFHandler<dim>::active_cell_iterator
        level_cell = mg_dof_handler.begin_active(level);
       const typename MGDoFHandler<dim>::active_cell_iterator
@@ -366,9 +380,16 @@ copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
                                           // into the level-wise one
          for (unsigned int i=0; i<dofs_per_cell; ++i)
            {
-             if (target_component[fe.system_to_component_index(i).first] == selected)
-               dst[level](level_dof_indices[i] - level_start)
-                 = src(global_dof_indices[i] - start);
+             const unsigned int component
+               = mg_target_component[fe.system_to_component_index(i).first];
+             if (selected[component])
+               {
+                 const unsigned int level_start
+                   = mg_component_start[level][component];
+                 
+                 dst[level](level_dof_indices[i] - level_start)
+                   = src(global_dof_indices[i] - offset);
+               }
            }
          
          for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
@@ -402,13 +423,68 @@ copy_to_mg (const MGDoFHandler<dim>        &mg_dof_handler,
 
 
 
+template <typename number>
+template <int dim, typename number2>
+void
+MGTransferSelect<number>::copy_from_mg (
+  const MGDoFHandler<dim>&              mg_dof_handler,
+  BlockVector<number2>&                 dst,
+  const MGLevelObject<Vector<number> >& src) const
+{
+  do_copy_from_mg (mg_dof_handler, dst, src, 0);
+}
+
+
+
+template <typename number>
+template <int dim, typename number2>
+void
+MGTransferSelect<number>::copy_from_mg (
+  const MGDoFHandler<dim>&              mg_dof_handler,
+  Vector<number2>&                      dst,
+  const MGLevelObject<Vector<number> >& src) const
+{
+  do_copy_from_mg (mg_dof_handler, dst, src,
+                  component_start[selected_component]);
+}
+
+
+
+template <typename number>
+template <int dim, typename number2>
+void
+MGTransferSelect<number>::copy_from_mg_add (
+  const MGDoFHandler<dim>&              mg_dof_handler,
+  BlockVector<number2>&                 dst,
+  const MGLevelObject<Vector<number> >& src) const
+{
+  do_copy_from_mg_add (mg_dof_handler, dst, src, 0);
+}
+
+
+
+template <typename number>
+template <int dim, typename number2>
+void
+MGTransferSelect<number>::copy_from_mg_add (
+  const MGDoFHandler<dim>&              mg_dof_handler,
+  Vector<number2>&                      dst,
+  const MGLevelObject<Vector<number> >& src) const
+{
+  do_copy_from_mg_add (mg_dof_handler, dst, src,
+                      component_start[selected_component]);
+}
+
+
+
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferSelect<number>::
-copy_from_mg (const MGDoFHandler<dim>              &mg_dof_handler,
-              OutVector                            &dst,
-              const MGLevelObject<Vector<number> > &src) const
+MGTransferSelect<number>::do_copy_from_mg (
+  const MGDoFHandler<dim>              &mg_dof_handler,
+  OutVector                            &dst,
+  const MGLevelObject<Vector<number> > &src,
+  const unsigned int offset) const
 {
 
   const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
@@ -422,10 +498,6 @@ copy_from_mg (const MGDoFHandler<dim>              &mg_dof_handler,
   const typename MGDoFHandler<dim>::active_cell_iterator
     endc = mg_dof_handler.end();
 
-                                  // Start of input vector inside a
-                                  // block vector.
-  const unsigned int start = component_start[selected];
-
                                   // traverse all cells and copy the
                                   // data appropriately to the output
                                   // vector
@@ -437,10 +509,6 @@ copy_from_mg (const MGDoFHandler<dim>              &mg_dof_handler,
        DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
       const unsigned int level = level_cell->level();
 
-                                      // Start of treated component
-                                      // for this level
-      const unsigned int level_start = mg_component_start[level][selected];
-      
                                       // get the dof numbers of
                                       // this cell for the global
                                       // and the level-wise
@@ -451,9 +519,17 @@ copy_from_mg (const MGDoFHandler<dim>              &mg_dof_handler,
                                       // copy level-wise data to
                                       // global vector
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-       if (target_component[fe.system_to_component_index(i).first] == selected)
-         dst(global_dof_indices[i]-start)
-           = src[level](level_dof_indices[i]-level_start);
+       {
+         const unsigned int component
+           = mg_target_component[fe.system_to_component_index(i).first];
+       if (selected[component])
+         {
+           const unsigned int level_start
+             = mg_component_start[level][component];
+           dst(global_dof_indices[i] - offset)
+             = src[level](level_dof_indices[i]-level_start);
+         }
+       }
     };
 
                                   // clear constrained nodes
@@ -465,10 +541,11 @@ copy_from_mg (const MGDoFHandler<dim>              &mg_dof_handler,
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferSelect<number>::
-copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
-                  OutVector                            &dst,
-                  const MGLevelObject<Vector<number> > &src) const
+MGTransferSelect<number>::do_copy_from_mg_add (
+  const MGDoFHandler<dim>              &mg_dof_handler,
+  OutVector                            &dst,
+  const MGLevelObject<Vector<number> > &src,
+  const unsigned int offset) const
 {
 
   const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
@@ -477,10 +554,6 @@ copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
   std::vector<unsigned int> global_dof_indices (dofs_per_cell);
   std::vector<unsigned int> level_dof_indices (dofs_per_cell);
 
-                                  // Start of input vector inside a
-                                  // block vector.
-  const unsigned int start = component_start[selected];
-
   typename MGDoFHandler<dim>::active_cell_iterator
     level_cell = mg_dof_handler.begin_active();
   const typename MGDoFHandler<dim>::active_cell_iterator
@@ -497,10 +570,6 @@ copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
        DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
       const unsigned int level = level_cell->level();
 
-                                      // Start of treated component
-                                      // for this level
-      const unsigned int level_start = mg_component_start[level][selected];
-      
                                       // get the dof numbers of
                                       // this cell for the global
                                       // and the level-wise
@@ -508,18 +577,25 @@ copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
       global_cell.get_dof_indices (global_dof_indices);
       level_cell->get_mg_dof_indices(level_dof_indices);
 
-//TODO:[GK] Probably wrong for continuous elements
+//TODO:[GK+WB] Probably wrong for continuous elements
 
                                       // copy level-wise data to
                                       // global vector
       for (unsigned int i=0; i<dofs_per_cell; ++i)
-       if (target_component[fe.system_to_component_index(i).first] == selected)
-         dst(global_dof_indices[i]-start)
-           += src[level](level_dof_indices[i]-level_start);
-    };
-
+       {
+         const unsigned int component
+           = mg_target_component[fe.system_to_component_index(i).first];
+         if (selected[component])
+           {
+             const unsigned int level_start
+               = mg_component_start[level][component];
+             dst(global_dof_indices[i] - offset)
+               += src[level](level_dof_indices[i] - level_start);
+           }
+       }
+    }
                                   // clear constrained nodes
-//TODO:[GK]  constraints->set_zero(dst);
+//TODO:[GK+WB]  constraints->set_zero(dst);
 }
 
 
@@ -531,10 +607,10 @@ copy_from_mg_add (const MGDoFHandler<dim>              &mg_dof_handler,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferBlock<number>::
-copy_to_mg (const MGDoFHandler<dim>             &mg_dof_handler,
-            MGLevelObject<BlockVector<number> > &dst,
-            const InVector                      &src) const
+MGTransferBlock<number>::copy_to_mg (
+  const MGDoFHandler<dim>             &mg_dof_handler,
+  MGLevelObject<BlockVector<number> > &dst,
+  const InVector                      &src) const
 {
                                   // forward to the correct
                                   // specialization
@@ -546,11 +622,11 @@ copy_to_mg (const MGDoFHandler<dim>             &mg_dof_handler,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferBlock<number>::
-copy_to_mg (const MGDoFHandler<dim>             &,
-            MGLevelObject<BlockVector<number> > &,
-            const InVector                      &,
-           const is_1d<true>                   &) const
+MGTransferBlock<number>::copy_to_mg (
+  const MGDoFHandler<dim>             &,
+  MGLevelObject<BlockVector<number> > &,
+  const InVector                      &,
+  const is_1d<true>                   &) const
 {
   Assert(false, ExcNotImplemented());
 }
@@ -560,11 +636,11 @@ copy_to_mg (const MGDoFHandler<dim>             &,
 template <typename number>
 template <int dim, class InVector>
 void
-MGTransferBlock<number>::
-copy_to_mg (const MGDoFHandler<dim>             &mg_dof_handler,
-            MGLevelObject<BlockVector<number> > &dst,
-            const InVector                      &src,
-           const is_1d<false>                  &) const
+MGTransferBlock<number>::copy_to_mg (
+  const MGDoFHandler<dim>             &mg_dof_handler,
+  MGLevelObject<BlockVector<number> > &dst,
+  const InVector                      &src,
+  const is_1d<false>                  &) const
 {
                                   // Make src a real finite element
                                   // function
@@ -582,7 +658,7 @@ copy_to_mg (const MGDoFHandler<dim>             &mg_dof_handler,
   
   dst.clear();
 
-  MGTools::reinit_vector(mg_dof_handler, dst, selected, target_component);
+  MGTools::reinit_vector(mg_dof_handler, dst, selected, mg_target_component);
   
   std::vector<unsigned int> global_dof_indices (dofs_per_cell);
   std::vector<unsigned int> level_dof_indices  (dofs_per_cell);
@@ -646,10 +722,10 @@ copy_to_mg (const MGDoFHandler<dim>             &mg_dof_handler,
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferBlock<number>::
-copy_from_mg (const MGDoFHandler<dim>                   &mg_dof_handler,
-              OutVector                                 &dst,
-              const MGLevelObject<BlockVector<number> > &src) const
+MGTransferBlock<number>::copy_from_mg (
+  const MGDoFHandler<dim>                   &mg_dof_handler,
+  OutVector                                 &dst,
+  const MGLevelObject<BlockVector<number> > &src) const
 {
 
   const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
@@ -693,10 +769,10 @@ copy_from_mg (const MGDoFHandler<dim>                   &mg_dof_handler,
 template <typename number>
 template <int dim, class OutVector>
 void
-MGTransferBlock<number>::
-copy_from_mg_add (const MGDoFHandler<dim>                   &mg_dof_handler,
-                  OutVector                                 &dst,
-                  const MGLevelObject<BlockVector<number> > &src) const
+MGTransferBlock<number>::copy_from_mg_add (
+  const MGDoFHandler<dim>                   &mg_dof_handler,
+  OutVector                                 &dst,
+  const MGLevelObject<BlockVector<number> > &src) const
 {
   const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
   const unsigned int dofs_per_cell = fe.dofs_per_cell;
index 42bc9452af846a7c894e3e07c35f3a090f509b48..d04abdde0c60d57b79aa399e75ba059dfe041773 100644 (file)
@@ -388,10 +388,6 @@ MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
 }
 
 
-//TODO:[GK] Check if this function can be simplified
-// It took quite a lot of try and error to make it work with
-// target_component, so it is more grown than designed
-
 template<int dim, typename number>
 void
 MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
@@ -418,14 +414,12 @@ MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
 
   Assert (selected.size() == target_component.size(),
          ExcDimensionMismatch(selected.size(), target_component.size()));
-
-  std::vector<bool> target_selected(selected.size(), false);
-  for (unsigned int i=0;i<selected.size();++i)
-    if (selected[i]) target_selected[target_component[i]] = true;
   
-  unsigned int n_selected = std::accumulate(target_selected.begin(),
-                                           target_selected.end(),
-                                           0U);
+                                  // Compute the number of blocks needed
+  const unsigned int n_selected
+    = std::accumulate(selected.begin(),
+                     selected.end(),
+                     0U);
   
   std::vector<std::vector<unsigned int> >
     ndofs(mg_dof.get_tria().n_levels(),
@@ -450,6 +444,46 @@ MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
 }
 
 
+template<int dim, typename number>
+void
+MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                        MGLevelObject<Vector<number> >& v,
+                        const std::vector<bool>& selected,
+                       const std::vector<unsigned int>& target_component)
+{
+  const unsigned int ncomp = mg_dof.get_fe().n_components();
+
+  Assert (selected.size() == target_component.size(),
+         ExcDimensionMismatch(selected.size(), target_component.size()));
+  
+                                  // Compute the number of blocks needed
+  if (DEBUG)
+    {
+      const unsigned int n_selected
+       = std::accumulate(selected.begin(),
+                         selected.end(),
+                         0U);
+      Assert(n_selected == 1, ExcDimensionMismatch(n_selected, 1));
+    }
+  
+  unsigned int selected_block = 0;
+  while (!selected[selected_block])
+    ++selected_block;
+  
+  std::vector<std::vector<unsigned int> >
+    ndofs(mg_dof.get_tria().n_levels(),
+         std::vector<unsigned int>(target_component.size()));
+
+  count_dofs_per_component (mg_dof, ndofs, target_component);
+  
+  for (unsigned int level=v.get_minlevel();
+       level<=v.get_maxlevel();++level)
+    {
+      v[level].reinit(ndofs[level][selected_block]);
+    }
+}
+
+
 
 
 
@@ -562,6 +596,17 @@ template void MGTools::reinit_vector<deal_II_dimension> (
   const std::vector<bool>&,
   const std::vector<unsigned int>&);
 
+template void MGTools::reinit_vector<deal_II_dimension> (
+  const MGDoFHandler<deal_II_dimension>&,
+  MGLevelObject<Vector<double> >&,
+  const std::vector<bool>&,
+  const std::vector<unsigned int>&);
+template void MGTools::reinit_vector<deal_II_dimension> (
+  const MGDoFHandler<deal_II_dimension>&,
+  MGLevelObject<Vector<float> >&,
+  const std::vector<bool>&,
+  const std::vector<unsigned int>&);
+
 
 template void MGTools::count_dofs_per_component<deal_II_dimension> (
   const MGDoFHandler<deal_II_dimension>&,
index b10ff85066cbb714d1a7ce66dc0890baa6dbc962..0562a86dbe102ecaf118c565667623a4c0d54c8a 100644 (file)
@@ -34,12 +34,16 @@ void MGTransferBlockBase::build_matrices (
   const MGDoFHandler<dim>& mg_dof,
   const std::vector<bool>& select)
 {
+                                  // Fill target component with
+                                  // standard values (identity) if it
+                                  // is empty
   if (target_component.size() == 0)
     {
       target_component.resize(mg_dof.get_fe().n_components());
       for (unsigned int i=0;i<target_component.size();++i)
        target_component[i] = i;
     } else {
+                                      // otherwise, check it for consistency
       Assert (target_component.size() == mg_dof.get_fe().n_components(),
              ExcDimensionMismatch(target_component.size(),
                                   mg_dof.get_fe().n_components()));
@@ -49,22 +53,77 @@ void MGTransferBlockBase::build_matrices (
          Assert(i<target_component.size(),
                 ExcIndexRange(i,0,target_component.size()));
        }
+    }
+                                  // Do the same for the multilevel
+                                  // components. These may be
+                                  // different.
+  if (mg_target_component.size() == 0)
+    {
+      mg_target_component.resize(mg_dof.get_fe().n_components());
+      for (unsigned int i=0;i<mg_target_component.size();++i)
+       mg_target_component[i] = i;
+    } else {
+      Assert (mg_target_component.size() == mg_dof.get_fe().n_components(),
+             ExcDimensionMismatch(mg_target_component.size(),
+                                  mg_dof.get_fe().n_components()));
+      
+      for (unsigned int i=0;i<mg_target_component.size();++i)
+       {
+         Assert(i<mg_target_component.size(),
+                ExcIndexRange(i,0,mg_target_component.size()));
+       }
     }
   
   const FiniteElement<dim>& fe = mg_dof.get_fe();
 
+                                  // Effective number of components
+                                  // is the maximum entry in
+                                  // mg_target_component. This
+                                  // assumes that the values in that
+                                  // vector don't have holes.
   const unsigned int n_components  =
-    *std::max_element(target_component.begin(), target_component.end()) + 1;
+    *std::max_element(mg_target_component.begin(), mg_target_component.end()) + 1;
   const unsigned int dofs_per_cell = fe.dofs_per_cell;  
   const unsigned int n_levels      = mg_dof.get_tria().n_levels();
   
+                                  // Selected refers to the component
+                                  // in mg_target_component.
   selected = select;
   
   Assert (selected.size() == fe.n_components(),
-         ExcDimensionMismatch(selected.size(), fe.n_components()))
-
+         ExcDimensionMismatch(selected.size(), fe.n_components()));
+  
+                                  // Compute the lengths of all blocks
   sizes.resize(n_levels);
-  MGTools::count_dofs_per_component(mg_dof, sizes, target_component);
+  MGTools::count_dofs_per_component(mg_dof, sizes, mg_target_component);
+
+                                  // Fill some index vectors
+                                  // for later use.
+  mg_component_start = sizes;
+                                  // Compute start indices from sizes
+  for (unsigned int l=0;l<mg_component_start.size();++l)
+    {
+      unsigned int k=0;
+      for (unsigned int i=0;i<mg_component_start[l].size();++i)
+       {
+         const unsigned int t=mg_component_start[l][i];
+         mg_component_start[l][i] = k;
+         k += t;
+       }
+    }
+
+  component_start.resize(target_component.size());
+  DoFTools::count_dofs_per_component(mg_dof, component_start, target_component);
+
+  unsigned int k=0;
+  for (unsigned int i=0;i<component_start.size();++i)
+    {
+      const unsigned int t=component_start[i];
+      component_start[i] = k;
+      k += t;
+    }
+
+// Building the matrices starts here!
   
                                   // reset the size of the array of
                                   // matrices. call resize(0) first,
@@ -181,7 +240,7 @@ void MGTransferBlockBase::build_matrices (
                          = fe.system_to_component_index(i).first;
                        const unsigned int jcomp
                          = fe.system_to_component_index(j).first;
-                       if ((icomp==jcomp) && selected[target_component[icomp]])
+                       if ((icomp==jcomp) && selected[mg_target_component[icomp]])
                          prolongation_sparsities[level]->add(dof_indices_child[i],
                                                              dof_indices_parent[j]);
                      };
@@ -217,7 +276,7 @@ void MGTransferBlockBase::build_matrices (
                      {
                        const unsigned int icomp = fe.system_to_component_index(i).first;
                        const unsigned int jcomp = fe.system_to_component_index(j).first;
-                       if ((icomp==jcomp) && selected[target_component[icomp]])
+                       if ((icomp==jcomp) && selected[mg_target_component[icomp]])
                          prolongation_matrices[level]->set(dof_indices_child[i],
                                                            dof_indices_parent[j],
                                                            prolongation(i,j));
@@ -225,31 +284,6 @@ void MGTransferBlockBase::build_matrices (
              };
          };
     };
-                                  // Finally, fill some index vectors
-                                  // for later use.
-  mg_component_start = sizes;
-                                  // Compute start indices from sizes
-  for (unsigned int l=0;l<mg_component_start.size();++l)
-    {
-      unsigned int k=0;
-      for (unsigned int i=0;i<mg_component_start[l].size();++i)
-       {
-         const unsigned int t=mg_component_start[l][i];
-         mg_component_start[l][i] = k;
-         k += t;
-       }
-    }
-
-  component_start.resize(n_components);
-  DoFTools::count_dofs_per_component(mg_dof, component_start);
-
-  unsigned int k=0;
-  for (unsigned int i=0;i<component_start.size();++i)
-    {
-      const unsigned int t=component_start[i];
-      component_start[i] = k;
-      k += t;
-    }
 #if defined(DEAL_PREFER_MATRIX_EZ) && 1 == 0
   deallog.push("Transfer");
   for (unsigned int level=0;level<n_levels-1; ++level)
@@ -271,12 +305,14 @@ template <int dim>
 void MGTransferBlock<number>::build_matrices (
   const MGDoFHandler<dim> &mg_dof,
   std::vector<bool> select,
-  const std::vector<unsigned int>& t_component)
+  const std::vector<unsigned int>& t_component,
+  const std::vector<unsigned int>& mg_t_component)
 {
   if (select.size() == 0)
     select = std::vector<bool> (mg_dof.get_fe().n_components(), true);
 
   target_component = t_component;
+  mg_target_component = mg_t_component;
 
   MGTransferBlockBase::build_matrices (mg_dof, select);
 }
@@ -287,15 +323,19 @@ template <int dim>
 void MGTransferSelect<number>::build_matrices (
   const MGDoFHandler<dim> &mg_dof,
   unsigned int select,
-  const std::vector<unsigned int>& t_component)
+  unsigned int mg_select,
+  const std::vector<unsigned int>& t_component,
+  const std::vector<unsigned int>& mg_t_component)
 {
   unsigned int ncomp = mg_dof.get_fe().n_components();
   
-  selected = select;
+  selected_component = select;
   std::vector<bool> s(ncomp, false);
-  s[select] = true;
-
+  s[mg_select] = true;
+//TODO:[GK] Assert that select and mg_select access the same component
+// or figure out how to use different components
   target_component = t_component;
+  mg_target_component = mg_t_component;
   
   MGTransferBlockBase::build_matrices (mg_dof, s);
 }
@@ -312,7 +352,7 @@ void MGTransferBlock<float>::build_matrices<deal_II_dimension>
 template
 void MGTransferSelect<float>::build_matrices<deal_II_dimension>
 (const MGDoFHandler<deal_II_dimension> &mg_dof,
- unsigned int, const std::vector<unsigned int>&);
+ unsigned int, unsigned int, const std::vector<unsigned int>&);
 
 template
 void MGTransferBlock<double>::build_matrices<deal_II_dimension>
@@ -322,7 +362,7 @@ void MGTransferBlock<double>::build_matrices<deal_II_dimension>
 template
 void MGTransferSelect<double>::build_matrices<deal_II_dimension>
 (const MGDoFHandler<deal_II_dimension> &mg_dof,
- unsigned int, const std::vector<unsigned int>&);
+ unsigned int, unsigned int, const std::vector<unsigned int>&);
 
 
 template void
index 093099414af8ac9777d1a69c0c07464780dd5a8d..70ba223c790e9886b1fe20b1adffae598168cc96 100644 (file)
@@ -185,7 +185,8 @@ void MGTransferSelect<number>::prolongate (
   Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
          ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
 
-      prolongation_matrices[to_level-1]->block(selected, selected)
+      prolongation_matrices[to_level-1]->block(selected_component,
+                                              selected_component)
        .vmult (dst, src);
 }
 
@@ -199,7 +200,8 @@ void MGTransferSelect<number>::restrict_and_add (
   Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
          ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
 
-  prolongation_matrices[from_level-1]->block(selected, selected)
+  prolongation_matrices[from_level-1]->block(selected_component,
+                                            selected_component)
     .Tvmult_add (dst, src);
 }
 
index 77353a49e2b5d1194750079fc05c9b54b526a726..6eba0bc4ff1ec4bbcab724bd5fc1cfddb50e9bdb 100644 (file)
@@ -18,6 +18,7 @@
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <dofs/dof_renumbering.h>
+#include <dofs/dof_tools.h>
 #include <fe/fe_dgp.h>
 #include <fe/fe_dgq.h>
 #include <fe/fe_q.h>
 #include <multigrid/mg_dof_handler.h>
 #include <multigrid/mg_transfer.h>
 #include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_level_object.h>
 
 #include <fstream>
 #include <iostream>
+#include <iomanip>
 #include <algorithm>
 
 
@@ -67,6 +70,7 @@ void check_simple(const FiniteElement<dim>& fe)
 
 template <int dim>
 void check_select(const FiniteElement<dim>& fe,
+                 unsigned int selected,
                  std::vector<unsigned int> target_component,
                  std::vector<unsigned int> mg_target_component)
 {
@@ -82,24 +86,24 @@ void check_select(const FiniteElement<dim>& fe,
   for (unsigned int l=0;l<tr.n_levels();++l)
     DoFRenumbering::component_wise(mgdof, l, mg_target_component);
   
-  std::vector<std::vector<unsigned int> > ndofs(mgdof.get_tria().n_levels());
-  MGTools::count_dofs_per_component(mgdof, ndofs, mg_target_component);
-  for (unsigned int l=0;l<ndofs.size();++l)
+  std::vector<std::vector<unsigned int> > mg_ndofs(mgdof.get_tria().n_levels());
+  MGTools::count_dofs_per_component(mgdof, mg_ndofs, mg_target_component);
+  for (unsigned int l=0;l<mg_ndofs.size();++l)
     {
       deallog << "Level " << l << " dofs:";
-      for (unsigned int i=0;i<ndofs[l].size();++i)
-       deallog << ' ' << ndofs[l][i];
+      for (unsigned int i=0;i<mg_ndofs[l].size();++i)
+       deallog << ' ' << mg_ndofs[l][i];
       deallog << std::endl;
     }
   
   
   MGTransferSelect<double> transfer;
-  transfer.build_matrices(mgdof, 0,
-                         mg_target_component);
+  transfer.build_matrices(mgdof, selected, selected,
+                         target_component, mg_target_component);
 
-  Vector<double> u2(ndofs[2][0]);
-  Vector<double> u1(ndofs[1][0]);
-  Vector<double> u0(ndofs[0][0]);
+  Vector<double> u2(mg_ndofs[2][0]);
+  Vector<double> u1(mg_ndofs[1][0]);
+  Vector<double> u0(mg_ndofs[0][0]);
 
   u0 = 1;
   transfer.prolongate(1,u1,u0);
@@ -114,6 +118,23 @@ void check_select(const FiniteElement<dim>& fe,
   deallog << "\tu1 " << u1.l2_norm()
          << "\tu0 " << u0.l2_norm()
          << std::endl;
+
+                                  // Check copy to mg and back
+  std::vector<unsigned int> ndofs (target_component.size());
+  DoFTools::count_dofs_per_component(mgdof, ndofs, target_component);
+  BlockVector<double> u;
+  u.reinit (ndofs);
+  for (unsigned int i=0;i<u.size();++i)
+    u(i) = i;
+  
+  MGLevelObject<Vector<double> > v;
+  v.resize(2,2);
+  v[2].reinit(mg_ndofs[2][selected]);
+
+  transfer.copy_to_mg(mgdof, v, u);
+  for (unsigned int i=0; i<v[2].size();++i)
+    deallog << ' ' << (int) v[2](i);
+  deallog << std::endl;
 }
 
 int main()
@@ -121,7 +142,7 @@ int main()
   std::ofstream logfile("transfer.output");
   logfile.precision(3);
   deallog.attach(logfile);
-  deallog.depth_console(100);
+  deallog.depth_console(0);
   
 //  check_simple (FE_DGP<2>(0));
 //  check_simple (FE_DGP<2>(1));
@@ -140,7 +161,14 @@ int main()
       v3[i] = i/2;
     }
   
-  check_select (FESystem<2>(FE_DGQ<2>(1), 4), v1, v1);
-  check_select (FESystem<2>(FE_DGQ<2>(1), 4), v2, v2);
-  check_select (FESystem<2>(FE_DGQ<2>(1), 4), v3, v3);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 0, v1, v1);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 1, v1, v1);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 0, v2, v2);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 0, v3, v3);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 1, v3, v3);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 0, v1, v3);
+  check_select (FESystem<2>(FE_DGQ<2>(1), 4), 1, v1, v3);
 }
+
+
+

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.