]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
multigrid cleanup branch merged
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Nov 2001 11:55:28 +0000 (11:55 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 29 Nov 2001 11:55:28 +0000 (11:55 +0000)
git-svn-id: https://svn.dealii.org/trunk@5307 0785d39b-7218-0410-832d-ea1e28bc413d

14 files changed:
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_smoother.h
deal.II/deal.II/include/multigrid/mg_smoother.templates.h
deal.II/deal.II/include/multigrid/mg_tools.templates.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/mg_transfer.h [new file with mode: 0644]
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/source/dofs/dof_renumbering.cc
deal.II/deal.II/source/multigrid/mg_base.cc
deal.II/deal.II/source/multigrid/mg_dof_tools.cc
deal.II/deal.II/source/multigrid/mg_smoother.cc
deal.II/deal.II/source/multigrid/multigrid.all_dimensions.cc
deal.II/deal.II/source/multigrid/multigrid.cc
tests/multigrid/mg1.cc

index 3c99b73188e49a9bc9e24326f7f2a2f89c980695..00e12ab33ef7656d0b634be21da7f1185e37bd5a 100644 (file)
 #include <base/subscriptor.h>
 #include <base/smartpointer.h>
 #include <lac/vector.h>
-#include <lac/sparsity_pattern.h>
-
 #include <vector>
 
-class MGSmootherBase;
-
 
 /**
  * An array with an object for each level.  The purpose of this class
@@ -105,39 +101,6 @@ class MGLevelObject : public Subscriptor
 };
 
 
-/**
- * Abstract base class for coarse grid solvers.  The interface of a
- * function call operator is defined to execute coarse grid solution
- * in a derived class.
- *
- * @author Guido Kanschat, 1999
- */
-class MGCoarseGridSolver : public Subscriptor
-{
-  public:
-                                    /**
-                                     * Virtual destructor. Does
-                                     * nothing in particular, but
-                                     * needs to be declared anyhow.
-                                     */
-    virtual ~MGCoarseGridSolver();
-    
-                                    /**
-                                     * Coarse grid solving method.
-                                     * This is only the interface for
-                                     * a function defined in a derived
-                                     * class like
-                                     * @p{MGCoarseGridLACIteration}.
-                                     *
-                                     * Note that the information
-                                     * about the matrix is removed to
-                                     * that class.
-                                     */
-    virtual void operator() (const unsigned int    level,
-                            Vector<double>       &dst,
-                            const Vector<double> &src) const = 0;
-};
-
 
 /**
  * Coarse grid solver using LAC iterative methods.
@@ -150,8 +113,8 @@ class MGCoarseGridSolver : public Subscriptor
  *
  * @author Guido Kanschat, 1999
  */
-template<class SOLVER, class MATRIX, class PRECOND>
-class MGCoarseGridLACIteration :  public MGCoarseGridSolver
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR = Vector<double> >
+class MGCoarseGridLACIteration :  public Subscriptor
 {
   public:
                                     /**
@@ -171,9 +134,9 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
                                      * matrix, vectors and
                                      * preconditioner.
                                      */
-    virtual void operator() (const unsigned int   level,
-                            Vector<double>       &dst,
-                            const Vector<double> &src) const;
+    void operator() (const unsigned int   level,
+                    VECTOR       &dst,
+                    const VECTOR &src) const;
   private:
                                     /**
                                      * Reference to the solver.
@@ -200,303 +163,63 @@ class MGCoarseGridLACIteration :  public MGCoarseGridSolver
  *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGTransferBase : public Subscriptor
-{
-  public:
-                                    /**
-                                     * Destructor. Does nothing here, but
-                                     * needs to be declared virtual anyway.
-                                     */
-    virtual ~MGTransferBase();
-
-                                    /**
-                                     * Prolongate a vector from level
-                                     * @p{to_level-1} to level
-                                     * @p{to_level}. The previous
-                                     * content of @p{dst} is
-                                     * overwritten.
-                                     *
-                                     * @p{src} is assumed to be a vector with
-                                     * as many elements as there are degrees
-                                     * of freedom on the coarser level of
-                                     * the two involved levels, while @p{src}
-                                     * shall have as many elements as there
-                                     * are degrees of freedom on the finer
-                                     * level.
-                                     */
-    virtual void prolongate (const unsigned int    to_level,
-                            Vector<double>       &dst,
-                            const Vector<double> &src) const = 0;
-
-                                    /**
-                                     * Restrict a vector from level
-                                     * @p{from_level} to level
-                                     * @p{from_level-1} and add this
-                                     * restriction to
-                                     * @p{dst}. Obviously, if the
-                                     * refined region on level
-                                     * @p{from_level} is smaller than
-                                     * that on level @p{from_level-1},
-                                     * some degrees of freedom in
-                                     * @p{dst} are not covered and will
-                                     * not be altered. For the other
-                                     * degress of freedom, the result
-                                     * of the restriction is added.
-                                     *
-                                     * @p{src} is assumed to be a vector with
-                                     * as many elements as there are degrees
-                                     * of freedom on the finer level of
-                                     * the two involved levels, while @p{src}
-                                     * shall have as many elements as there
-                                     * are degrees of freedom on the coarser
-                                     * level.
-                                     */
-    virtual void restrict_and_add (const unsigned int    from_level,
-                                  Vector<double>       &dst,
-                                  const Vector<double> &src) const = 0;
-};
-
-
-/**
- * Basic class for preconditioning by multigrid.
- *
- * The functionality of the multigrid method is restricted to defect
- * correction. It is @em{not} iterative and the start solution is
- * always zero. Since by this $u^E_l$ and $u^A_l$ (see report on
- * multigrid) are always zero, restriction is simplified a lot and
- * maybe even the seam condition on grids is oblivious. Still, I am
- * not sure that these restrictions on the class might cause numerical
- * inefficiencies.
- *
- * The function @p{precondition} is the actual multigrid method and
- * makes use of several operations to be implemented in derived
- * classes. It takes a defect in @p{src} and the result is the multigrid
- * preconditioned defect in @p{dst}.
- *
- * @author Guido Kanschat, 1999
- */
-class MGBase : public Subscriptor
-{
-  private:
-                                    /**
-                                     * Copy constructor. Made private
-                                     * to prevent use.
-                                     */
-    MGBase(const MGBase&);
-
-                                    /**
-                                     * Copy operator. Made private
-                                     * to prevent use.
-                                     */
-    const MGBase& operator=(const MGBase&);
-
-  public:
-                                    /**
-                                     * Constructor, subject to change.
-                                     */
-    MGBase (const MGTransferBase &transfer,
-           const unsigned int    minlevel,
-           const unsigned int    maxlevel);
-    
-                                    /**
-                                     * Virtual destructor.
-                                     */
-    virtual ~MGBase();
-    
-                                    /**
-                                     * Execute one step of the
-                                     * v-cycle algorithm.  This
-                                     * function assumes, that the
-                                     * vector @p{d} is properly filled
-                                     * with the residual in the outer
-                                     * defect correction
-                                     * scheme. After execution of
-                                     * @p{vcycle()}, the result is in
-                                     * the vector @p{s}. We propose to
-                                     * write functions @p{copy_from_mg}
-                                     * and @p{copy_to_mg} in derived
-                                     * classes of @p{MGBase} to access
-                                     * @p{d} and @p{s}.
-                                     *
-                                     * The actual work for this
-                                     * function is done in
-                                     * @p{level_mgstep}.
-                                     */
-    void vcycle(const MGSmootherBase     &pre_smooth,
-               const MGSmootherBase     &post_smooth,
-               const MGCoarseGridSolver &cgs);
-
-                                    /**
-                                     * Exception.
-                                     */
-    DeclException2(ExcSwitchedLevels, int, int,
-                  << "minlevel and maxlevel switched, should be: "
-                  << arg1 << "<=" << arg2);
-    
-  protected:
-                                    /**
-                                     * Highest level of cells.
-                                     */
-    unsigned int maxlevel;
-
-                                    /**
-                                     * Level for coarse grid solution.
-                                     */
-    unsigned int minlevel;
-    
-                                    /**
-                                     * Auxiliary vector.
-                                     */
-    MGLevelObject<Vector<double> > defect;
-
-                                    /**
-                                     * Auxiliary vector.
-                                     */
-    MGLevelObject<Vector<double> > solution;
-    
-                                    /**
-                                     * Auxiliary vector.
-                                     */
-    MGLevelObject<Vector<double> > t;    
-
-                                    /**
-                                     * Prolongation and restriction object.
-                                     */
-    SmartPointer<const MGTransferBase> transfer;
-
-                                    /**
-                                     * The actual v-cycle multigrid method.
-                                     * @p{level} is the level the
-                                     * function should work on. It
-                                     * will usually be called for the
-                                     * highest level from outside,
-                                     * but will then call itself
-                                     * recursively for @p{level-1},
-                                     * unless we are on @p{minlevel}
-                                     * where instead of recursing
-                                     * deeper, the coarse grid solver
-                                     * is used to solve the matrix of
-                                     * this level.
-                                     */
-    void level_mgstep (const unsigned int        level,
-                      const MGSmootherBase     &pre_smooth,
-                      const MGSmootherBase     &post_smooth,
-                      const MGCoarseGridSolver &cgs);  
-
-                                    /**
-                                     * Apply negative @p{vmult} on all
-                                     * cells of a level.
-                                     * This is implemented in a
-                                     * derived class.
-                                     */
-    virtual void level_vmult (const unsigned int    level,
-                             Vector<double>       &dst,
-                             const Vector<double> &src,
-                             const Vector<double> &rhs) = 0;
-
-                                    /**
-                                     * Additional multiplication on a
-                                     * refinement edge.  Using DGFEM,
-                                     * the application of the global
-                                     * matrix to a function on the
-                                     * fine level produces results on
-                                     * the surrounding cells of the
-                                     * coarse level. Therefore,
-                                     * additionally to the fine level
-                                     * matrix, we need an operator
-                                     * from the fine level to the
-                                     * coarse level, taking care of
-                                     * this multiplication.
-                                     *
-                                     * Like @ref{level_vmult}, it is
-                                     * expected to add the negative
-                                     * product of the matrix and
-                                     * @p{src} to @p{dst}. Here,
-                                     * @p{src} is a vector on the
-                                     * fine level and @p{dst} is a
-                                     * vector on the coarse level.
-                                     *
-                                     * This function has an empty
-                                     * implementation here and must
-                                     * be overloaded in a multigrid
-                                     * class for discontinuous
-                                     * methods.
-                                     */
-    virtual void edge_vmult (const unsigned int    level,
-                            Vector<double>       &dst,
-                            const Vector<double> &src);
-
-                                    /**
-                                     * Print a level vector using
-                                     * @ref{DoFHandler}.
-                                     */
-    virtual void print_vector (const unsigned int level,
-                              const Vector<double>& v,
-                              const char* name) const = 0;
-    
-};
-
-
-/**
- * Multi-level preconditioner.
- * Here, we collect all information needed for multi-level preconditioning
- * and provide the standard interface for LAC iterative methods.
- *
- * The template parameter class @p{MG} is required to inherit @p{MGBase}.
- * Furthermore, it needs functions @p{void copy_to_mg(const VECTOR&)}
- * to store @p{src} in the right hand side of the multi-level method and
- * @p{void copy_from_mg(VECTOR&)} to store the result of the v-cycle in @p{dst}.
- *
- * @author Guido Kanschat, 1999
- */
-template<class MG, class VECTOR = Vector<double> >
-class PreconditionMG : public Subscriptor
-{
-  public:
-                                    /**
-                                     * Constructor.
-                                     * Arguments are the multigrid object,
-                                     * pre-smoother, post-smoother and
-                                     * coarse grid solver.
-                                     */
-    PreconditionMG(MG                       &mg,
-                  const MGSmootherBase     &pre,
-                  const MGSmootherBase     &post,
-                  const MGCoarseGridSolver &coarse);
-    
-                                    /**
-                                     * Preconditioning operator,
-                                     * calling the @p{vcycle} function
-                                     * of the @p{MG} object passed to
-                                     * the constructor.
-                                     *
-                                     * This is the operator used by
-                                     * LAC iterative solvers.
-                                     */
-    void vmult (VECTOR       &dst,
-               const VECTOR &src) const;
-    
-  private:
-                                    /**
-                                     * The multigrid object.
-                                     */
-    SmartPointer<MG> multigrid;
-    
-                                    /**
-                                     * The pre-smoothing object.
-                                     */
-    SmartPointer<const MGSmootherBase> pre;
+//  template <class VECTOR>
+//  class MGTransferBase : public Subscriptor
+//  {
+//    public:
+//                                  /**
+//                                   * Destructor. Does nothing here, but
+//                                   * needs to be declared virtual anyway.
+//                                   */
+//      virtual ~MGTransferBase();
+
+//                                  /**
+//                                   * Prolongate a vector from level
+//                                   * @p{to_level-1} to level
+//                                   * @p{to_level}. The previous
+//                                   * content of @p{dst} is
+//                                   * overwritten.
+//                                   *
+//                                   * @p{src} is assumed to be a vector with
+//                                   * as many elements as there are degrees
+//                                   * of freedom on the coarser level of
+//                                   * the two involved levels, while @p{src}
+//                                   * shall have as many elements as there
+//                                   * are degrees of freedom on the finer
+//                                   * level.
+//                                   */
+//      virtual void prolongate (const unsigned int to_level,
+//                          VECTOR&            dst,
+//                          const VECTOR&      src) const = 0;
+
+//                                  /**
+//                                   * Restrict a vector from level
+//                                   * @p{from_level} to level
+//                                   * @p{from_level-1} and add this
+//                                   * restriction to
+//                                   * @p{dst}. Obviously, if the
+//                                   * refined region on level
+//                                   * @p{from_level} is smaller than
+//                                   * that on level @p{from_level-1},
+//                                   * some degrees of freedom in
+//                                   * @p{dst} are not covered and will
+//                                   * not be altered. For the other
+//                                   * degress of freedom, the result
+//                                   * of the restriction is added.
+//                                   *
+//                                   * @p{src} is assumed to be a vector with
+//                                   * as many elements as there are degrees
+//                                   * of freedom on the finer level of
+//                                   * the two involved levels, while @p{src}
+//                                   * shall have as many elements as there
+//                                   * are degrees of freedom on the coarser
+//                                   * level.
+//                                   */
+//      virtual void restrict_and_add (const unsigned int from_level,
+//                                VECTOR&            dst,
+//                                const VECTOR&      src) const = 0;
+//  };
 
-                                    /**
-                                     * The post-smoothing object.
-                                     */
-    SmartPointer<const MGSmootherBase> post;
-    
-                                    /**
-                                     * The coarse grid solver.
-                                     */
-    SmartPointer<const MGCoarseGridSolver> coarse;
-};
 
 
 /* ------------------------------------------------------------------- */
@@ -574,10 +297,11 @@ MGLevelObject<Object>::get_maxlevel () const
 /* ------------------ Functions for MGCoarseGridLACIteration ------------ */
 
 
-template<class SOLVER, class MATRIX, class PRECOND>
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::MGCoarseGridLACIteration(SOLVER        &s,
-                                                                         const MATRIX  &m,
-                                                                         const PRECOND &p)
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR>
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND, VECTOR>
+::MGCoarseGridLACIteration(SOLVER& s,
+                          const MATRIX  &m,
+                          const PRECOND &p)
                :
                solver(s),
                matrix(&m),
@@ -585,43 +309,17 @@ MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::MGCoarseGridLACIteration(SOLVER
 {};
 
 
-template<class SOLVER, class MATRIX, class PRECOND>
+template<class SOLVER, class MATRIX, class PRECOND, class VECTOR>
 void
-MGCoarseGridLACIteration<SOLVER,MATRIX,PRECOND>::operator() (const unsigned int    /* level */,
-                                                            Vector<double>       &dst,
-                                                            const Vector<double> &src) const
+MGCoarseGridLACIteration<SOLVER, MATRIX, PRECOND, VECTOR>
+::operator() (const unsigned int    /* level */,
+             VECTOR       &dst,
+             const VECTOR &src) const
 {
   solver.solve(*matrix, dst, src, precondition);
 }
 
 
-/* ------------------------------------------------------------------- */
-
-
-template<class MG, class VECTOR>
-PreconditionMG<MG, VECTOR>::PreconditionMG(MG                       &mg,
-                                          const MGSmootherBase     &pre,
-                                          const MGSmootherBase     &post,
-                                          const MGCoarseGridSolver &coarse)
-               :
-               multigrid(&mg),
-               pre(&pre),
-               post(&post),
-               coarse(&coarse)
-{}
-
-
-template<class MG, class VECTOR>
-void
-PreconditionMG<MG,VECTOR>::vmult (VECTOR       &dst,
-                                 const VECTOR &src) const
-{
-  multigrid->copy_to_mg(src);
-  multigrid->vcycle(*pre, *post, *coarse);
-  multigrid->copy_from_mg(dst);
-}
-
-
 /*----------------------------   mgbase.h     ---------------------------*/
 
 #endif
index 561874100b80fb2571e813fdb746120cd7e06160..bd3802e9dc8c8feca8467bc66acae0457549b05b 100644 (file)
 #ifndef __deal2__mg_dof_tools_h
 #define __deal2__mg_dof_tools_h
 
+template <class Object> class MGLevelObject;
+template <int dim> class MGDoFHandler;
+template <typename number> class Vector;
+template <typename number> class BlockVector;
+template <class number> class FullMatrix;
 
 #include <base/config.h>
 
  * for more information.
  *
  * All member functions are static, so there is no need to create an
- * object of class @p{MGDoFTools}.
+ * object of class @p{MGTools}.
  *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999, 2000
  */
-class MGDoFTools 
+class MGTools
 {
   public:
                                     /**
@@ -47,7 +52,7 @@ class MGDoFTools
                                      * hanging nodes here, since only
                                      * one level is considered.
                                      */
-    template <int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_sparsity_pattern (const MGDoFHandler<dim> &dof_handler,
                           SparsityPattern         &sparsity,
@@ -59,7 +64,7 @@ class MGDoFTools
                                      * @see make_sparsity_pattern
                                      * $see DoFTools
                                      */
-    template <int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof_handler,
                                SparsityPattern         &sparsity,
@@ -72,7 +77,7 @@ class MGDoFTools
                                      * @see{make_flux_sparsity_pattern}
                                      * @see{DoFTools}
                                      */
-    template <int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof_handler,
                                     SparsityPattern         &sparsity,
@@ -95,7 +100,7 @@ class MGDoFTools
                                      * for the couplings occuring in
                                      * fluxes.
                                      */
-    template <int dim>
+    template <int dim, class SparsityPattern>
     static void
     make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
                                SparsityPattern       &sparsity,
@@ -103,6 +108,63 @@ class MGDoFTools
                                const FullMatrix<double>& int_mask,
                                const FullMatrix<double>& flux_mask);
 
+                                    /**
+                                     * Count the dofs component-wise
+                                     * on each level.
+                                     *
+                                     * Result is a vector containing
+                                     * for each level a vector
+                                     * containing the number of dofs
+                                     * for each component (access is
+                                     * @p{result[level][component]}).
+                                     */
+    template <int dim>
+      static void count_dofs_per_component (const MGDoFHandler<dim>& mg_dof,
+                                           std::vector<std::vector<unsigned int> >& result);
+    
+    
+                                    /**
+                                     * Ajust vectors on all levels to
+                                     * correct size.  Here, we just
+                                     * count the numbers of degrees
+                                     * of freedom on each level and
+                                     * @p{reinit} each level vector
+                                     * to this length.
+                                     */
+    template <int dim, typename number>
+      static void
+      reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                    MGLevelObject<Vector<number> >& vector);
+
+                                    /**
+                                     * Adjust block-vectors on all
+                                     * levels to correct size.  Count
+                                     * the numbers of degrees of
+                                     * freedom on each level
+                                     * component-wise. Then, assign
+                                     * each block of @p{vector} the
+                                     * corresponding size.
+                                     *
+                                     * The boolean field @p{selected}
+                                     * allows restricting this
+                                     * operation to certain
+                                     * components. In this case,
+                                     * @p{vector} will only have as
+                                     * many blocks as there are true
+                                     * values in @p{selected} (no
+                                     * blocks of length zero are
+                                     * padded in).
+                                     *
+                                     * Degrees of freedom must be
+                                     * sorted by component in order
+                                     * to obtain reasonable results
+                                     * from this function.
+                                     */
+    template <int dim, typename number>
+      static void
+      reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                    MGLevelObject<BlockVector<number> >& vector,
+                    const vector<bool>& selected);
 };
 
 
index 715e0e35f139d78cf447d6969c9727420e7b672f..0c9fbeae82dd6642923068cdcc4fd1ef03f76e81 100644 (file)
 
 template <int dim> class MGDoFHandler;
 
-/**
- * Abstract base class for multigrid smoothers.
- * In fact, this class only provides the interface of the smoothing function.
- * Using @p{deal.II} grid handling, @p{MGSmoother} is a good point to start.
- *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
- */
-class MGSmootherBase :  public Subscriptor 
-{  
-  public:
-                                    /**
-                                     * Virtual destructor. Does
-                                     * nothing in particular, but
-                                     * needs to be declared anyhow.
-                                     */
-    virtual ~MGSmootherBase();
-    
-                                    /**
-                                     * Smooth the residual of @p{u} on
-                                     * the given level. If $S$ is the
-                                     * smoothing operator, then this
-                                     * function should do the
-                                     * following operation:
-                                     * $u += S (rhs - Au)$, where @p{u} and
-                                     * @p{rhs} are the input parameters.
-                                     *
-                                     * This function should keep the
-                                     * interior level boundary
-                                     * values, so you may want to
-                                     * call @p{set_zero_interior_boundary}
-                                     * in the derived class
-                                     * @p{MGSmoother} somewhere in your
-                                     * derived function, or another
-                                     * function doing similar things.
-                                     */
-    virtual void smooth (const unsigned int    level,
-                        Vector<double>       &u,
-                        const Vector<double> &rhs) const = 0;
-};
-
-
 /**
  * Smoother doing nothing. This class is not useful for many applications other
  * than for testing some multigrid procedures. Also some applications might
@@ -71,12 +30,12 @@ class MGSmootherBase :  public Subscriptor
  *
  * @author Guido Kanschat, 1999
  */
-class MGSmootherIdentity : public MGSmootherBase
+class MGSmootherIdentity : public Subscriptor
 {
   public:
                                     /**
                                      * Implementation of the
-                                     * interface in @p{MGSmootherBase}.
+                                     * interface for @p{Multigrid}.
                                      * This function does nothing,
                                      * which by comparison with the
                                      * definition of this function
@@ -84,33 +43,33 @@ class MGSmootherIdentity : public MGSmootherBase
                                      * operator equals the null
                                      * operator.
                                      */
-    virtual void smooth (const unsigned int   level,
-                        Vector<double>       &u,
-                        const Vector<double> &rhs) const;
+  template <class VECTOR>
+  void smooth (const unsigned int level,
+              VECTOR&            u,
+              const VECTOR&      rhs) const;
 };
 
 
 /**
- * Base class for multigrid smoothers. It declares the interface
- * to smoothers by inheriting @p{MGSmootherBase} and implements some functionality for setting
- * the values
- * of vectors at interior boundaries (i.e. boundaries between differing
- * levels of the triangulation) to zero, by building a list of these degrees
- * of freedom's indices at construction time.
+ * Base class for multigrid smoothers. It implements some
+ * functionality for setting the values of vectors at interior
+ * boundaries (i.e. boundaries between differing levels of the
+ * triangulation) to zero, by building a list of these degrees of
+ * freedom's indices at construction time.
  *
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-class MGSmoother : public MGSmootherBase
+class MGSmootherContinuous : public Subscriptor
 {
   private:
                                     /**
                                      * Default constructor. Made private
                                      * to prevent it being called, which
                                      * is necessary since this could
-                                     * circumpass the set-up of the list
+                                     * circumvent the set-up of the list
                                      * if interior boundary dofs.
                                      */
-    MGSmoother ();
+    MGSmootherContinuous ();
     
   public:
 
@@ -122,15 +81,19 @@ class MGSmoother : public MGSmootherBase
                                      * needed by the function
                                      * @p{set_zero_interior_boundaries}.
                                      *
-                                     * Since this function is implemented
-                                     * a bit different in 1d (there are no
-                                     * faces of cells, just vertices),
-                                     * there are actually two sets of
-                                     * constructors, namely this one for 1d
-                                     * and the following one for all other
-                                     * dimensions.
+                                     * Since this function is
+                                     * implemented a bit different in
+                                     * 1d (there are no faces of
+                                     * cells, just vertices), there
+                                     * are actually two sets of
+                                     * constructors, namely this one
+                                     * for 1d and the following one
+                                     * for all other
+                                     * dimensions. Really amusing
+                                     * about this text is, that there
+                                     * is no 1d implementation.
                                      */
-    MGSmoother (const MGDoFHandler<1> &mg_dof,
+    MGSmootherContinuous (const MGDoFHandler<1> &mg_dof,
                unsigned int           steps);
 
                                     /**
@@ -145,7 +108,7 @@ class MGSmoother : public MGSmootherBase
                                      * steps to be executed by @p{smooth}.
                                      */
     template <int dim>
-    MGSmoother (const MGDoFHandler<dim> &mg_dof,
+    MGSmootherContinuous (const MGDoFHandler<dim> &mg_dof,
                unsigned int             steps);    
 
                                     /**
@@ -158,12 +121,15 @@ class MGSmoother : public MGSmootherBase
                                      * has no interior boundaries, this
                                      * function does nothing in this case.
                                      */
-    void set_zero_interior_boundary (const unsigned int  level,
-                                    Vector<double>      &u) const;
+    template <class VECTOR>
+    void set_zero_interior_boundary (const unsigned int level,
+                                    VECTOR&            u) const;
+    
                                     /**
                                      * Modify the number of smoothing steps.
                                      */
     void set_steps(unsigned int steps);
+    
                                     /**
                                      * How many steps should be used?
                                      */
@@ -174,6 +140,7 @@ class MGSmoother : public MGSmootherBase
                                      * Number of smoothing steps.
                                      */
     unsigned int steps;
+    
                                     /**
                                      * For each level, we store a list of
                                      * degree of freedom indices which are
@@ -191,27 +158,36 @@ class MGSmoother : public MGSmootherBase
 
 
 /**
- * Implementation of a smoother using matrix builtin relaxation methods.
+ * A smoother using matrix builtin relaxation methods.
+ *
+ * Additionally to smoothing with the matrix built-in relaxation
+ * scheme, hanging nodes of continuous finite elements are handled
+ * properly (at least in a future version).
+ *
+ * The library contains instantiation for @p{SparseMatrix<.>} and
+ * @p{Vector<.>}, where the template arguments are all combinations of
+ * @p{float} and @p{double}. Additional instantiations may be created
+ * by including the file mg_smoother.templates.h.
  * 
  * @author Wolfgang Bangerth, Guido Kanschat, 1999
  */
-template<typename number>
-class MGSmootherRelaxation : public MGSmoother
+template<class MATRIX, class VECTOR>
+class MGSmootherRelaxation : public MGSmootherContinuous
 {
   public:
                                     /**
                                      * Type of the smoothing
                                      * function of the matrix.
                                      */
-    typedef void
-    (SparseMatrix<number>::* function_ptr)(Vector<double>&, const Vector<double>&,
-                                          typename SparseMatrix<number>::value_type) const;
+    typedef void (MATRIX::* function_ptr)(VECTOR&,
+                                         const VECTOR&,
+                                         typename MATRIX::value_type) const;
     
                                     /**
                                      * Constructor.
                                      * The constructor uses an @p{MGDoFHandler} to initialize
                                      * data structures for interior level boundary handling
-                                     * in @p{MGSmoother}.
+                                     * in @p{MGSmootherContinuous}.
                                      *
                                      * Furthermore, it takes a pointer to the
                                      * level matrices and their smoothing function.
@@ -225,25 +201,26 @@ class MGSmootherRelaxation : public MGSmoother
                                      */
 
     template<int dim>
-    MGSmootherRelaxation(const MGDoFHandler<dim>                       &mg_dof,
-                        const MGLevelObject<SparseMatrix<number> > &matrix,
-                        const function_ptr                             function,
-                        const unsigned int                             steps,
-                        const double                                   omega = 1.);
+    MGSmootherRelaxation(const MGDoFHandler<dim>&     mg_dof,
+                        const MGLevelObject<MATRIX>& matrix,
+                        const function_ptr           function,
+                        const unsigned int           steps,
+                        const double                 omega = 1.);
     
                                     /**
-                                     * Implementation of the interface in @p{MGSmootherBase}.
-                                     * We use the SOR method in @p{SparseMatrix} for the real work
-                                     * and find a way to keep the boundary values.
+                                     * We use the relaxation method in
+                                     * @p{MATRIX} for the real
+                                     * work and find a way to keep
+                                     * the boundary values.
                                      */
-    virtual void smooth (const unsigned int   level,
-                        Vector<double>       &u,
-                        const Vector<double> &rhs) const;
+    void smooth (const unsigned int level,
+                VECTOR&            u,
+                const VECTOR&      rhs) const;
   private:
                                     /**
                                      * Pointer to the matrices.
                                      */
-    SmartPointer<const MGLevelObject<SparseMatrix<number> > > matrix;
+    SmartPointer<const MGLevelObject<MATRIX> > matrix;
     
                                     /**
                                      * Pointer to the relaxation function.
@@ -269,16 +246,23 @@ class MGSmootherRelaxation : public MGSmoother
 
 /* ------------------------------- Inline functions -------------------------- */
 
+template <class VECTOR>
+inline void
+MGSmootherIdentity::smooth (const unsigned int, VECTOR&, const VECTOR&) const
+{};
+
+/*----------------------------------------------------------------------*/
+
 inline
 void
-MGSmoother::set_steps(unsigned int i)
+MGSmootherContinuous::set_steps(unsigned int i)
 {
   steps = i;
 }
 
 inline
 unsigned int
-MGSmoother::get_steps() const
+MGSmootherContinuous::get_steps() const
 {
   return steps;
 }
index 6d3393085433a5384e297c20d8f60ca77a0adba7..846cd969ed37e9f56cc03c878aa9f2a0b6804dc5 100644 (file)
@@ -1,26 +1,28 @@
 // $Id$
 
 
-template<typename number>
+template<class MATRIX, class VECTOR>
 template<int dim>
-MGSmootherRelaxation<number>::MGSmootherRelaxation (const MGDoFHandler<dim>                       &mg_dof,
-                                                   const MGLevelObject<SparseMatrix<number> > &matrix,
-                                                   const function_ptr                             relaxation,
-                                                   const unsigned int                             steps,
-                                                   const double                                   omega)
+MGSmootherRelaxation<MATRIX, VECTOR>
+::MGSmootherRelaxation (const MGDoFHandler<dim>&     mg_dof,
+                       const MGLevelObject<MATRIX>& matrix,
+                       const function_ptr           relaxation,
+                       const unsigned int           steps,
+                       const double                 omega)
                :
-               MGSmoother(mg_dof, steps),
+               MGSmootherContinuous(mg_dof, steps),
                matrix(&matrix),
                relaxation(relaxation),
                omega(omega)
 {};
 
 
-template<typename number>
+template<class MATRIX, class VECTOR>
 void
-MGSmootherRelaxation<number>::smooth (const unsigned int   level,
-                                     Vector<double>       &u,
-                                     const Vector<double> &rhs) const
+MGSmootherRelaxation<MATRIX, VECTOR>
+::smooth (const unsigned int level,
+         VECTOR&            u,
+         const VECTOR&      rhs) const
 {
   h1.reinit(u);
   h2.reinit(u);
diff --git a/deal.II/deal.II/include/multigrid/mg_tools.templates.h b/deal.II/deal.II/include/multigrid/mg_tools.templates.h
new file mode 100644 (file)
index 0000000..ca1aac5
--- /dev/null
@@ -0,0 +1,471 @@
+//----------------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------------
+
+#ifndef __deal2__mg_tools_templates_h
+#define __deal2__mg_tools_templates_h
+
+#include <dofs/dof_constraints.h>
+#include <numerics/data_out.h>
+#include <multigrid/multigrid.h>
+#include <algorithm>
+#include <fstream>
+
+#include <lac/sparse_matrix.h>
+
+
+//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces
+#if deal_II_dimension == 1
+
+template <int dim, class InVector>
+void
+MGTransferPrebuilt::copy_to_mg (const MGDoFHandler<dim>&,
+                               MGLevelObject<Vector<double> >&,
+                               const InVector&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+#else
+
+
+template <int dim, class InVector>
+void
+MGTransferPrebuilt::copy_to_mg (const MGDoFHandler<dim>& mg_dof_handler,
+                               MGLevelObject<Vector<double> >& dst,
+                               const InVector& src) const
+{
+
+                                  // Make src a real finite element function
+//  InVector src = osrc;
+//  constraints->distribute(src);
+  const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
+  const unsigned int dofs_per_face = mg_dof_handler.get_fe().dofs_per_face;
+
+                                  // set the elements of the vectors
+                                  // on all levels to zero
+  unsigned int minlevel = dst.get_minlevel();
+  unsigned int maxlevel = dst.get_maxlevel();
+  
+  dst.clear();
+
+  for (unsigned int l=minlevel;l<=maxlevel;++l)
+    dst[l].reinit(sizes[l]);
+  
+  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);
+
+                                  // traverse the grid top-down
+                                  // (i.e. starting with the most
+                                  // refined grid). this way, we can
+                                  // always get that part of one
+                                  // level of the output vector which
+                                  // corresponds to a region which is
+                                  // more refined, by restriction of
+                                  // the respective vector on the
+                                  // next finer level, which we then
+                                  // already have built.
+  for (int level=maxlevel; level>=static_cast<int>(minlevel); --level)
+    {
+      typename MGDoFHandler<dim>::active_cell_iterator
+       level_cell = mg_dof_handler.begin_active(level);
+      const typename MGDoFHandler<dim>::active_cell_iterator
+       level_end  = mg_dof_handler.end_active(level);
+
+//TODO:[?] Treat hanging nodes properly
+// The single-level vector is not an FE-function, because the values at
+// hanging nodes are set to zero. This should be treated before the restriction.
+
+                                      // Compute coarse level right hand side
+                                      // by restricting from fine level.
+      for (; level_cell!=level_end; ++level_cell)
+       {
+         DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+                                          // get the dof numbers of
+                                          // this cell for the global
+                                          // and the level-wise
+                                          // numbering
+         global_cell.get_dof_indices(global_dof_indices);
+         level_cell->get_mg_dof_indices (level_dof_indices);
+
+                                          // transfer the global
+                                          // defect in the vector
+                                          // into the level-wise one
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           dst[level](level_dof_indices[i]) = src(global_dof_indices[i]);
+
+         for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
+           {
+             const typename MGDoFHandler<dim>::face_iterator
+               face = level_cell->face(face_n);
+             if (face->has_children())
+               {
+                 face->get_mg_dof_indices(level_face_indices);
+
+
+                                                  // Delete values on refinement edge,
+                                                  // since restriction will add them again.
+                 for (unsigned int i=0; i<dofs_per_face; ++i)
+                   dst[level](level_face_indices[i])
+                     = 0.;
+               };
+           };
+       };
+                                      // for that part of the level
+                                      // which is further refined:
+                                      // get the defect by
+                                      // restriction of the defect on
+                                      // one level higher
+      if (static_cast<unsigned int>(level) < maxlevel)
+       {
+         restrict_and_add (level+1, dst[level], dst[level+1]);
+       }
+    };
+};
+
+#endif
+
+
+template <int dim, class OutVector>
+void
+MGTransferPrebuilt::copy_from_mg(const MGDoFHandler<dim>& mg_dof_handler,
+                                OutVector &dst,
+                                const MGLevelObject<Vector<double> > &src) const
+{
+  const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
+
+  std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+  std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+  typename MGDoFHandler<dim>::active_cell_iterator
+    level_cell = mg_dof_handler.begin_active();
+  const typename MGDoFHandler<dim>::active_cell_iterator
+    endc = mg_dof_handler.end();
+
+                                  // traverse all cells and copy the
+                                  // data appropriately to the output
+                                  // vector
+
+                                  // Is the level monotonuosly increasing?
+
+  for (; level_cell != endc; ++level_cell)
+    {
+      DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+      const unsigned int level = level_cell->level();
+      
+                                      // get the dof numbers of
+                                      // this cell for the global
+                                      // and the level-wise
+                                      // numbering
+      global_cell.get_dof_indices (global_dof_indices);
+      level_cell->get_mg_dof_indices(level_dof_indices);
+
+                                      // copy level-wise data to
+                                      // global vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       dst(global_dof_indices[i]) = src[level](level_dof_indices[i]);
+    };
+
+                                  // clear constrained nodes
+//TODO:  constraints->set_zero(dst);
+}
+
+template <int dim, class OutVector>
+void
+MGTransferPrebuilt::copy_from_mg_add(const MGDoFHandler<dim>& mg_dof_handler,
+                                    OutVector &dst,
+                                    const MGLevelObject<Vector<double> > &src) const
+{
+  const unsigned int dofs_per_cell = mg_dof_handler.get_fe().dofs_per_cell;
+
+  std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+  std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+  typename MGDoFHandler<dim>::active_cell_iterator
+    level_cell = mg_dof_handler.begin_active();
+  const typename MGDoFHandler<dim>::active_cell_iterator
+    endc = mg_dof_handler.end();
+
+                                  // traverse all cells and copy the
+                                  // data appropriately to the output
+                                  // vector
+
+                                  // Is the level monotonuosly increasing?
+
+  for (; level_cell != endc; ++level_cell)
+    {
+      DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+      const unsigned int level = level_cell->level();
+      
+                                      // get the dof numbers of
+                                      // this cell for the global
+                                      // and the level-wise
+                                      // numbering
+      global_cell.get_dof_indices (global_dof_indices);
+      level_cell->get_mg_dof_indices(level_dof_indices);
+
+//TODO: Probably wrong for continuous elements
+
+                                      // copy level-wise data to
+                                      // global vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       dst(global_dof_indices[i]) += src[level](level_dof_indices[i]);
+    };
+
+                                  // clear constrained nodes
+//TODO:  constraints->set_zero(dst);
+}
+
+
+//TODO:[?] This function needs to be specially implemented, since in 2d mode we use faces
+#if deal_II_dimension == 1
+
+template <int dim, class InVector>
+void
+MGTransferSelect::copy_to_mg (const MGDoFHandler<dim>&,
+                               MGLevelObject<Vector<double> >&,
+                               const InVector&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+#else
+
+
+template <int dim, class InVector>
+void
+MGTransferSelect::copy_to_mg (const MGDoFHandler<dim>& mg_dof_handler,
+                               MGLevelObject<Vector<double> >& dst,
+                               const InVector& osrc) const
+{
+                                  // Make src a real finite element function
+  InVector src = osrc;
+//  constraints->distribute(src);
+
+  const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int dofs_per_face = fe.dofs_per_face;
+
+                                  // set the elements of the vectors
+                                  // on all levels to zero
+  unsigned int minlevel = dst.get_minlevel();
+  unsigned int maxlevel = dst.get_maxlevel();
+  
+  dst.clear();
+  
+  for (unsigned int l=minlevel;l<=maxlevel;++l)
+    dst[l].reinit(sizes[l][selected]);
+  
+  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
+                                  // always get that part of one
+                                  // level of the output vector which
+                                  // corresponds to a region which is
+                                  // more refined, by restriction of
+                                  // the respective vector on the
+                                  // next finer level, which we then
+                                  // already have built.
+  for (int level=maxlevel; level>=static_cast<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
+       level_end  = mg_dof_handler.end_active(level);
+
+//TODO:[?] Treat hanging nodes properly
+// The single-level vector is not an FE-function, because the values at
+// hanging nodes are set to zero. This should be treated before the restriction.
+
+                                      // Compute coarse level right hand side
+                                      // by restricting from fine level.
+      for (; level_cell!=level_end; ++level_cell)
+       {
+         DoFObjectAccessor<dim, dim>& global_cell = *level_cell;
+                                          // get the dof numbers of
+                                          // this cell for the global
+                                          // and the level-wise
+                                          // numbering
+         global_cell.get_dof_indices(global_dof_indices);
+         level_cell->get_mg_dof_indices (level_dof_indices);
+
+                                          // transfer the global
+                                          // defect in the vector
+                                          // into the level-wise one
+         for (unsigned int i=0; i<dofs_per_cell; ++i)
+           {
+             if (fe.system_to_component_index(i).first == selected)
+               dst[level](level_dof_indices[i] - level_start)
+                 = src(global_dof_indices[i] - start);
+           }
+         
+         for (unsigned int face_n=0; face_n<GeometryInfo<dim>::faces_per_cell; ++face_n)
+           {
+             const typename MGDoFHandler<dim>::face_iterator
+               face = level_cell->face(face_n);
+             if (face->has_children())
+               {
+                 face->get_mg_dof_indices(level_face_indices);
+
+
+                                                  // Delete values on refinement edge,
+                                                  // since restriction will add them again.
+//               for (unsigned int i=0; i<dofs_per_face; ++i)
+//                 dst[level](level_face_indices[i])
+//                   = 0.;
+               };
+           };
+       };
+                                      // for that part of the level
+                                      // which is further refined:
+                                      // get the defect by
+                                      // restriction of the defect on
+                                      // one level higher
+      if (static_cast<unsigned int>(level) < maxlevel)
+       {
+         restrict_and_add (level+1, dst[level], dst[level+1]);
+       }
+    };
+};
+
+#endif
+
+
+template <int dim, class OutVector>
+void
+MGTransferSelect::copy_from_mg(const MGDoFHandler<dim>& mg_dof_handler,
+                                OutVector &dst,
+                                const MGLevelObject<Vector<double> > &src) const
+{
+
+  const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+  std::vector<unsigned int> global_dof_indices (dofs_per_cell);
+  std::vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+  typename MGDoFHandler<dim>::active_cell_iterator
+    level_cell = mg_dof_handler.begin_active();
+  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
+
+                                  // Is the level monotonuosly increasing?
+
+  for (; level_cell != endc; ++level_cell)
+    {
+       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
+                                      // numbering
+      global_cell.get_dof_indices (global_dof_indices);
+      level_cell->get_mg_dof_indices(level_dof_indices);
+
+                                      // copy level-wise data to
+                                      // global vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       if (fe.system_to_component_index(i).first == selected)
+         dst(global_dof_indices[i]-start)
+           = src[level](level_dof_indices[i]-level_start);
+    };
+
+                                  // clear constrained nodes
+//TODO:  constraints->set_zero(dst);
+}
+
+template <int dim, class OutVector>
+void
+MGTransferSelect::copy_from_mg_add(const MGDoFHandler<dim>& mg_dof_handler,
+                                    OutVector &dst,
+                                    const MGLevelObject<Vector<double> > &src) const
+{
+
+  const FiniteElement<dim>& fe = mg_dof_handler.get_fe();
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;
+
+  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
+    endc = mg_dof_handler.end();
+
+                                  // traverse all cells and copy the
+                                  // data appropriately to the output
+                                  // vector
+
+                                  // Is the level monotonuosly increasing?
+
+  for (; level_cell != endc; ++level_cell)
+    {
+       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
+                                      // numbering
+      global_cell.get_dof_indices (global_dof_indices);
+      level_cell->get_mg_dof_indices(level_dof_indices);
+
+//TODO: Probably wrong for continuous elements
+
+                                      // copy level-wise data to
+                                      // global vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       if (fe.system_to_component_index(i).first == selected)
+         dst(global_dof_indices[i]-start)
+           += src[level](level_dof_indices[i]-level_start);
+    };
+
+                                  // clear constrained nodes
+//TODO:  constraints->set_zero(dst);
+}
+
+
+#endif
diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h
new file mode 100644 (file)
index 0000000..1aeda70
--- /dev/null
@@ -0,0 +1,450 @@
+//----------------------------------------------------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2001 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------------------------------------------------
+#ifndef __deal2__mg_transfer_h
+#define __deal2__mg_transfer_h
+
+
+#include <lac/sparsity_pattern.h>
+#include <lac/block_vector.h>
+#include <lac/block_sparsity_pattern.h>
+#include <multigrid/mg_base.h>
+
+template <int dim> class MGDoFHandler;
+
+
+/**
+ * Implementation of the @p{MGTransferBase} interface for which the transfer
+ * operations are prebuilt upon construction of the object of this class as
+ * matrices. This is the fast way, since it only needs to build the operation
+ * once by looping over all cells and storing the result in a matrix for
+ * each level, but requires additional memory.
+ *
+ * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ */
+class MGTransferPrebuilt : public Subscriptor // MGTransferBase<Vector<double> > 
+{
+  public:
+                                    /**
+                                     * Destructor.
+                                     */
+    virtual ~MGTransferPrebuilt ();
+    
+                                    /**
+                                     * Actually build the prolongation
+                                     * matrices for each level.
+                                     */
+    template <int dim>
+    void build_matrices (const MGDoFHandler<dim> &mg_dof);
+
+                                    /**
+                                     * Prolongate a vector from level
+                                     * @p{to_level-1} to level
+                                     * @p{to_level}. The previous
+                                     * content of @p{dst} is
+                                     * overwritten.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the coarser level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the finer
+                                     * level.
+                                     */
+    virtual void prolongate (const unsigned int    to_level,
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const;
+
+                                    /**
+                                     * Restrict a vector from level
+                                     * @p{from_level} to level
+                                     * @p{from_level-1} and add this
+                                     * restriction to
+                                     * @p{dst}. Obviously, if the
+                                     * refined region on level
+                                     * @p{from_level} is smaller than
+                                     * that on level @p{from_level-1},
+                                     * some degrees of freedom in
+                                     * @p{dst} are not covered and will
+                                     * not be altered. For the other
+                                     * degress of freedom, the result
+                                     * of the restriction is added.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the finer level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the coarser
+                                     * level.
+                                     */
+    virtual void restrict_and_add (const unsigned int    from_level,
+                                  Vector<double>       &dst,
+                                  const Vector<double> &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, class InVector>
+    void
+    copy_to_mg (const MGDoFHandler<dim>& mg_dof,
+               MGLevelObject<Vector<double> > &dst,
+               const InVector &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<double>}. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom
+                                     * are set to zero.
+                                     */
+    template<int dim, class OutVector>
+    void
+    copy_from_mg (const MGDoFHandler<dim>& mg_dof,
+                 OutVector &dst,
+                 const MGLevelObject<Vector<double> > &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, class OutVector>
+    void
+    copy_from_mg_add (const MGDoFHandler<dim>& mg_dof,
+                     OutVector &dst,
+                     const MGLevelObject<Vector<double> > &src) const;
+
+  private:
+
+                                  /**
+                                   * Selected component.
+                                   */
+    unsigned int selected;
+
+                                  /**
+                                   * Sizes of the multi-level vectors.
+                                   */
+    std::vector<unsigned int> sizes;
+
+
+    std::vector<SparsityPattern>   prolongation_sparsities;
+
+                                    /**
+                                     * The actual prolongation matrix.
+                                     * column indices belong to the
+                                     * dof indices of the mother cell,
+                                     * i.e. the coarse level.
+                                     * while row indices belong to the
+                                     * child cell, i.e. the fine level.
+                                     */
+    std::vector<SparseMatrix<double> > prolongation_matrices;
+};
+
+
+/**
+ * Implementation of matrix generation for @ref{MGTransferBlock} and @p{MGTransferSelected}.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferBlockBase
+{
+  protected:  
+                                    /**
+                                     * Actually build the prolongation
+                                     * matrices for each level.
+                                     *
+                                     * If a field selected is given,
+                                     * only matrices for these
+                                     * components are to be built. By
+                                     * default, all matrices are
+                                     * built.
+                                     */
+    template <int dim>
+    void build_matrices (const MGDoFHandler<dim> &mg_dof,
+                        std::vector<bool> selected);
+
+                                  /**
+                                   * Flag of selected components.
+                                   */
+    std::vector<bool> selected;
+
+                                  /**
+                                   * Sizes of the multi-level vectors.
+                                   */
+    std::vector<std::vector<unsigned int> > sizes;
+  
+                                  /**
+                                   * Start index of each component.
+                                   */
+    std::vector<unsigned int> component_start;
+  
+                                  /**
+                                   * Start index of eah component on
+                                   * all levels.
+                                   */
+    std::vector<std::vector<unsigned int> > mg_component_start;
+
+  private:
+
+    std::vector<BlockSparsityPattern>   prolongation_sparsities;
+
+  protected:
+  
+                                    /**
+                                     * The actual prolongation matrix.
+                                     * column indices belong to the
+                                     * dof indices of the mother cell,
+                                     * i.e. the coarse level.
+                                     * while row indices belong to the
+                                     * child cell, i.e. the fine level.
+                                     */
+    std::vector<BlockSparseMatrix<double> > prolongation_matrices;
+};
+
+
+
+/**
+ * Implementation of the @p{MGTransferBase} interface for block
+ * matrices and block vectors.  In addition to the functionality of
+ * @ref{MGTransferPrebuilt}, the operation may be restricted to
+ * certain blocks of the vector.
+ *
+ * If the restricted mode is chosen, block vectors used in the
+ * transfer routines may only have as many components as there are
+ * @p{true}s in the selected-field.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferBlock : public Subscriptor, // MGTransferBase<BlockVector<double> >,
+                       private MGTransferBlockBase
+{
+  public:
+                                    /**
+                                     * Destructor.
+                                     */
+    virtual ~MGTransferBlock ();
+    
+                                    /**
+                                     * Build the prolongation
+                                     * matrices for each level.
+                                     *
+                                     * If a field selected is given,
+                                     * only matrices for these
+                                     * components are to be built. By
+                                     * default, all matrices are
+                                     * built.
+                                     *
+                                     * This function is a front-end
+                                     * for the same function in
+                                     * @ref{MGTransferBlockBase}.
+                                     */
+    template <int dim>
+    void build_matrices (const MGDoFHandler<dim> &mg_dof,
+                        std::vector<bool> selected = std::vector<bool>());
+
+                                    /**
+                                     * Prolongate a vector from level
+                                     * @p{to_level-1} to level
+                                     * @p{to_level}. The previous
+                                     * content of @p{dst} is
+                                     * overwritten.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the coarser level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the finer
+                                     * level.
+                                     */
+    virtual void prolongate (const unsigned int    to_level,
+                            BlockVector<double>       &dst,
+                            const BlockVector<double> &src) const;
+
+                                    /**
+                                     * Restrict a vector from level
+                                     * @p{from_level} to level
+                                     * @p{from_level-1} and add this
+                                     * restriction to
+                                     * @p{dst}. Obviously, if the
+                                     * refined region on level
+                                     * @p{from_level} is smaller than
+                                     * that on level @p{from_level-1},
+                                     * some degrees of freedom in
+                                     * @p{dst} are not covered and will
+                                     * not be altered. For the other
+                                     * degress of freedom, the result
+                                     * of the restriction is added.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the finer level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the coarser
+                                     * level.
+                                     */
+    virtual void restrict_and_add (const unsigned int    from_level,
+                                  BlockVector<double>       &dst,
+                                  const BlockVector<double> &src) const;
+
+  protected: 
+};
+
+
+
+/**
+ * Implementation of the @p{MGTransferBase} interface for block
+ * matrices and simple vectors. This class uses @ref{MGTransferBlock}
+ * selecting a single component. The transfer operators themselves are
+ * implemented for simple vectors again.
+ *
+ * @author Guido Kanschat, 2001
+ */
+class MGTransferSelect : public Subscriptor, // MGTransferBase<Vector<double> >,
+                        private MGTransferBlockBase
+{
+  public:
+                                    /**
+                                     * Destructor.
+                                     */
+    virtual ~MGTransferSelect ();
+    
+                                    /**
+                                     * Actually build the prolongation
+                                     * matrices for each level.
+                                     *
+                                     * Select the component you want
+                                     * to apply multigrid to.
+                                     *
+                                     * This function is a front-end
+                                     * for the same function in
+                                     * @ref{MGTransferBlockBase}.
+                                     */
+    template <int dim>
+    void build_matrices (const MGDoFHandler<dim> &mg_dof,
+                        unsigned int selected);
+
+                                    /**
+                                     * Prolongate a vector from level
+                                     * @p{to_level-1} to level
+                                     * @p{to_level}. The previous
+                                     * content of @p{dst} is
+                                     * overwritten.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the coarser level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the finer
+                                     * level.
+                                     */
+    virtual void prolongate (const unsigned int    to_level,
+                            Vector<double>       &dst,
+                            const Vector<double> &src) const;
+
+                                    /**
+                                     * Restrict a vector from level
+                                     * @p{from_level} to level
+                                     * @p{from_level-1} and add this
+                                     * restriction to
+                                     * @p{dst}. Obviously, if the
+                                     * refined region on level
+                                     * @p{from_level} is smaller than
+                                     * that on level @p{from_level-1},
+                                     * some degrees of freedom in
+                                     * @p{dst} are not covered and will
+                                     * not be altered. For the other
+                                     * degress of freedom, the result
+                                     * of the restriction is added.
+                                     *
+                                     * @p{src} is assumed to be a vector with
+                                     * as many elements as there are degrees
+                                     * of freedom on the finer level of
+                                     * the two involved levels, while @p{src}
+                                     * shall have as many elements as there
+                                     * are degrees of freedom on the coarser
+                                     * level.
+                                     */
+    virtual void restrict_and_add (const unsigned int    from_level,
+                                  Vector<double>       &dst,
+                                  const Vector<double> &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, class InVector>
+    void
+    copy_to_mg (const MGDoFHandler<dim>& mg_dof,
+               MGLevelObject<Vector<double> > &dst,
+               const InVector &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<double>}. In order to
+                                     * keep the result consistent,
+                                     * constrained degrees of freedom
+                                     * are set to zero.
+                                     */
+    template<int dim, class OutVector>
+    void
+    copy_from_mg (const MGDoFHandler<dim>& mg_dof,
+                 OutVector &dst,
+                 const MGLevelObject<Vector<double> > &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, class OutVector>
+    void
+    copy_from_mg_add (const MGDoFHandler<dim>& mg_dof,
+                     OutVector &dst,
+                     const MGLevelObject<Vector<double> > &src) const;
+
+  private:
+
+                                  /**
+                                   * Selected component.
+                                   */
+  unsigned int selected;
+};
+
+
+
+#endif
index 02acbd1cf254edfe62579055e61dbfb08bd42c7e..babcf2cc3cba4ae8b2b9880860c6b21bfab7f292 100644 (file)
 
 
 /**
- * Implementation of multigrid with matrices.
- * While @p{MGBase} was only an abstract framework for the V-cycle,
- * here we have the implementation of the pure virtual functions defined there.
- * Furthermore, level information is obtained from a triangulation object.
+ * Implementation of the multigrid method.
+ *
+ * The function actually performing a multi-level cycle,
+ * @p{level_mgstep}, as well as the function @p{vcycle}, calling it,
+ * require several helper classes handed over as template parameters.
+ * These classes have to meet the following requirements:
+ *
+ * @p{MATRIX} is a matrix as specified for LAC-solvers, that is, it
+ * has a function
+ * \begin{verbatim}
+ *   void MATRIX::vmult(VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * performing the matrix vector product @p{x=Ab}.
+ *
+ * @p{SMOOTH} is a class with a function
+ * \begin{verbatim}
+ *   void SMOOTH::smooth (unsigned int level, VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * modifying the vector @p{x} such that the residual @{b-Ax}
+ * on level @p{l} is smoothened. Refer to @ref{MGSmootherRelaxation}
+ * for an example.
+ *
+ * @p{COARSE} has an
+ * \begin{verbatim}
+ *   void COARSE::operator() (VECTOR& x, const VECTOR& b)
+ * \end{verbatim}
+ * returning the solution to the system @p{Ax=b} on the coarsest level
+ * in @p{x}.
  * 
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ * @author Guido Kanschat, 1999, 2001
  */
-template <int dim>
-class Multigrid : public MGBase
+template <class VECTOR>
+class Multigrid : public Subscriptor
 {
   public:
+  typedef VECTOR vector_type;
+  typedef const VECTOR const_vector_type;
+  
                                     /**
-                                     * Constructor. Take the object
-                                     * describing the multilevel
-                                     * hierarchy of degrees of
-                                     * freedom, an object describing
-                                     * constraints to degrees of
-                                     * freedom on the global grid
-                                     * (for example for hanging
-                                     * nodes), a set of matrices
-                                     * which describe the equation
-                                     * discretized on each level,
-                                     * their respective sparsity
-                                     * patterns, and an object
-                                     * mediating the transfer of
-                                     * solutions between different
-                                     * levels.
+                                     * Constructor. The
+                                     * @p{MGDoFHandler} is used to
+                                     * determine the highest possible
+                                     * level. @p{transfer} is an
+                                     * object performing prolongation
+                                     * and restriction.
+                                     *
+                                     * The V-cycle will start on
+                                     * level @p{maxlevel} and goes
+                                     * down to level @p{minlevel},
+                                     * where the coarse grid solver
+                                     * will be used.
                                      *
                                      * This function already
                                      * initializes the vectors which
@@ -63,189 +87,439 @@ class Multigrid : public MGBase
                                      * therefore create objects of
                                      * this type as late as possible.
                                      */
-//TODO:[GK] meaning of minlevel, maxlevel?
-    Multigrid(const MGDoFHandler<dim>                       &mg_dof_handler,
-             const ConstraintMatrix                        &constraints,
-             const MGLevelObject<SparsityPattern>       &level_sparsities,
-             const MGLevelObject<SparseMatrix<double> > &level_matrices,
-             const MGTransferBase                          &transfer,
-             const unsigned int                             minlevel = 0,
-             const unsigned int                             maxlevel = 1000000);
+  template <int dim>
+    Multigrid(const MGDoFHandler<dim>& mg_dof_handler,
+             const unsigned int            minlevel = 0,
+             const unsigned int            maxlevel = 1000000);
     
+
                                     /**
-                                     * Transfer from a vector on the
-                                     * global grid to vectors defined
-                                     * on each of the levels
-                                     * separately, i.a. an @p{MGVector}.
+                                     * Execute one step of the
+                                     * v-cycle algorithm.  This
+                                     * function assumes, that the
+                                     * vector @p{defect} is properly
+                                     * filled with the residual in
+                                     * the outer defect correction
+                                     * scheme (usually performed by
+                                     * @ref{PreconditionMG}). After
+                                     * execution of @p{vcycle()}, the
+                                     * result is in the vector
+                                     * @p{solution}. See
+                                     * @p{copy_*_mg} in class
+                                     * @p{MGTools} if you want to use
+                                     * these vectors yourself.
+                                     *
+                                     * The actual work for this
+                                     * function is done in
+                                     * @p{level_mgstep}.
                                      */
-    template<typename number>
-    void copy_to_mg (const Vector<number> &src);
+    template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+    void vcycle(const MGLevelObject<MATRIX>& matrix,
+               const TRANSFER&              transfer,
+               const SMOOTHER&              pre_smooth,
+               const SMOOTHER&              post_smooth,
+               const COARSE&                cgs);
 
                                     /**
-                                     * Transfer from multi-level vector to
-                                     * normal vector.
-                                     *
-                                     * Copies data from active
-                                     * portions of an MGVector into
-                                     * the respective positions of a
-                                     * @p{Vector<double>}. In order to
-                                     * keep the result consistent,
-                                     * constrained degrees of freedom
-                                     * are set to zero.
+                                     * Print a level vector using
+                                     * @ref{DoFHandler}.
                                      */
-    template<typename number>
-    void copy_from_mg (Vector<number> &dst) const;
+/*      virtual void print_vector (const unsigned int level, */
+/*                            const VECTOR& v, */
+/*                            const char* name) const; */
+    
 
+  private:
+    
+                                    /**
+                                     * The actual v-cycle multigrid method.
+                                     * @p{level} is the level the
+                                     * function should work on. It
+                                     * will usually be called for the
+                                     * highest level from outside,
+                                     * but will then call itself
+                                     * recursively for @p{level-1},
+                                     * unless we are on @p{minlevel}
+                                     * where instead of recursing
+                                     * deeper, the coarse grid solver
+                                     * is used to solve the matrix of
+                                     * this level.
+                                     */
+    template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+    void level_mgstep (const unsigned int           level,
+                      const MGLevelObject<MATRIX>& matrix,
+                      const TRANSFER&              transfer,
+                      const SMOOTHER&              pre_smooth,
+                      const SMOOTHER&              post_smooth,
+                      const COARSE&                cgs);
                                     /**
-                                     * Negative @p{vmult} on a level.
-                                     * @see MGBase.
+                                     * Level for coarse grid solution.
                                      */
-//TODO:[GK] what does this function do?    
-    virtual void level_vmult (const unsigned int    level,
-                             Vector<double>       &dest,
-                             const Vector<double> &src,
-                             const Vector<double> &rhs);
+    unsigned int minlevel;
 
                                     /**
-                                     * Print a level vector using
-                                     * @ref{DoFHandler}.
+                                     * Highest level of cells.
                                      */
-    virtual void print_vector (const unsigned int level,
-                              const Vector<double>& v,
-                              const char* name) const;
+    unsigned int maxlevel;
     
                                     /**
-                                     * Read-only access to level matrices.
+                                     * Auxiliary vector. Contains the
+                                     * defect to be corrected before
+                                     * the multigrid step.
                                      */
-    const SparseMatrix<double>& get_matrix (const unsigned int level) const;
+    MGLevelObject<VECTOR> defect;
 
-  private:
                                     /**
-                                     * Associated @p{MGDoFHandler}.
+                                     * Auxiliary vector. Contains the
+                                     * updated solution after the
+                                     * multigrid step.
                                      */
-    SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
-
-
-/**
- * Sparsity patterns for each level.
- */
-    SmartPointer<const MGLevelObject<SparsityPattern> > level_sparsities;
+    MGLevelObject<VECTOR> solution;
     
                                     /**
-                                     * Matrices for each level. The
-                                     * matrices are prepared by the
-                                     * constructor of @p{Multigrid} and
-                                     * can be accessed for
-                                     * assembling.
+                                     * Auxiliary vector.
                                      */
-    SmartPointer<const MGLevelObject<SparseMatrix<double> > > level_matrices;
+    MGLevelObject<VECTOR> t;    
+
+
+  private:
+    
 
                                     /**
-                                     * Pointer to the object holding
-                                     * constraints.
+                                     * Exception.
                                      */
-    SmartPointer<const ConstraintMatrix>                 constraints;
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
+    template<int dim, class MATRIX, class VECTOR2, class TRANSFER, class SMOOTHER, class COARSE> friend class PreconditionMG;
 };
 
 
 /**
- * Implementation of the @p{MGTransferBase} interface for which the transfer
- * operations are prebuilt upon construction of the object of this class as
- * matrices. This is the fast way, since it only needs to build the operation
- * once by looping over all cells and storing the result in a matrix for
- * each level, but requires additional memory.
+ * Multi-level preconditioner.
+ * Here, we collect all information needed for multi-level preconditioning
+ * and provide the standard interface for LAC iterative methods.
+ *
+ * The template parameter class @p{MG} is required to inherit @p{MGBase}.
+ * Furthermore, it needs functions @p{void copy_to_mg(const VECTOR&)}
+ * to store @p{src} in the right hand side of the multi-level method and
+ * @p{void copy_from_mg(VECTOR&)} to store the result of the v-cycle in @p{dst}.
  *
- * @author Wolfgang Bangerth, Guido Kanschat, 1999
+ * @author Guido Kanschat, 1999
  */
-class MGTransferPrebuilt : public MGTransferBase 
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+class PreconditionMG : public Subscriptor
 {
   public:
                                     /**
-                                     * Destructor.
+                                     * Constructor.
+                                     * Arguments are the multigrid object,
+                                     * pre-smoother, post-smoother and
+                                     * coarse grid solver.
                                      */
-    virtual ~MGTransferPrebuilt ();
+    PreconditionMG(const MGDoFHandler<dim>&     mg_dof,
+                  Multigrid<VECTOR>&           mg,
+                  const MGLevelObject<MATRIX>& matrix,
+                  const TRANSFER&              transfer,
+                  const SMOOTHER&              pre,
+                  const SMOOTHER&              post,
+                  const COARSE&                coarse);
     
                                     /**
-                                     * Actually build the prolongation
-                                     * matrices for each level.
+                                     * Preconditioning operator.
+                                     * Calls the @p{vcycle} function
+                                     * of the @p{MG} object passed to
+                                     * the constructor.
+                                     *
+                                     * This is the operator used by
+                                     * LAC iterative solvers.
                                      */
-    template <int dim>
-    void build_matrices (const MGDoFHandler<dim> &mg_dof);
-
+    void vmult (VECTOR       &dst,
+               const VECTOR &src) const;
+    
+                                    /**
+                                     * Preconditioning operator.
+                                     * Calls the @p{vcycle} function
+                                     * of the @p{MG} object passed to
+                                     * the constructor.
+                                     */
+  void vmult_add (VECTOR       &dst,
+                 const VECTOR &src) const;
+    
                                     /**
-                                     * Prolongate a vector from level
-                                     * @p{to_level-1} to level
-                                     * @p{to_level}. The previous
-                                     * content of @p{dst} is
-                                     * overwritten.
+                                     * Tranposed preconditioning operator.
                                      *
-                                     * @p{src} is assumed to be a vector with
-                                     * as many elements as there are degrees
-                                     * of freedom on the coarser level of
-                                     * the two involved levels, while @p{src}
-                                     * shall have as many elements as there
-                                     * are degrees of freedom on the finer
-                                     * level.
-                                     */
-    virtual void prolongate (const unsigned int    to_level,
-                            Vector<double>       &dst,
-                            const Vector<double> &src) const;
-
-                                    /**
-                                     * Restrict a vector from level
-                                     * @p{from_level} to level
-                                     * @p{from_level-1} and add this
-                                     * restriction to
-                                     * @p{dst}. Obviously, if the
-                                     * refined region on level
-                                     * @p{from_level} is smaller than
-                                     * that on level @p{from_level-1},
-                                     * some degrees of freedom in
-                                     * @p{dst} are not covered and will
-                                     * not be altered. For the other
-                                     * degress of freedom, the result
-                                     * of the restriction is added.
+                                     * Not implemented, but the
+                                     * definition may be needed.
+                                     */
+    void Tvmult (VECTOR       &dst,
+               const VECTOR &src) const;
+    
+                                    /**
+                                     * Tranposed preconditioning operator.
                                      *
-                                     * @p{src} is assumed to be a vector with
-                                     * as many elements as there are degrees
-                                     * of freedom on the finer level of
-                                     * the two involved levels, while @p{src}
-                                     * shall have as many elements as there
-                                     * are degrees of freedom on the coarser
-                                     * level.
-                                     */
-    virtual void restrict_and_add (const unsigned int    from_level,
-                                  Vector<double>       &dst,
-                                  const Vector<double> &src) const;
-
+                                     * Not implemented, but the
+                                     * definition may be needed.
+                                     */
+    void Tvmult_add (VECTOR       &dst,
+                    const VECTOR &src) const;
+    
   private:
+                                    /**
+                                     * Associated @p{MGDoFHandler}.
+                                     */
+    SmartPointer<const MGDoFHandler<dim> > mg_dof_handler;
+
+                                    /**
+                                     * The multigrid object.
+                                     */
+    SmartPointer<Multigrid<VECTOR> > multigrid;
+    
+                                    /**
+                                     * Matrices for each level. The
+                                     * matrices are prepared by the
+                                     * constructor of @p{Multigrid} and
+                                     * can be accessed for
+                                     * assembling.
+                                     */
+    SmartPointer<const MGLevelObject<MATRIX> > level_matrices;
 
-    std::vector<SparsityPattern>   prolongation_sparsities;
+                                  /**
+                                   * Object for grid tranfer.
+                                   */
+    SmartPointer<const TRANSFER> transfer;
+  
+                                    /**
+                                     * The pre-smoothing object.
+                                     */
+    SmartPointer<const SMOOTHER> pre;
 
                                     /**
-                                     * The actual prolongation matrix.
-                                     * column indices belong to the
-                                     * dof indices of the mother cell,
-                                     * i.e. the coarse level.
-                                     * while row indices belong to the
-                                     * child cell, i.e. the fine level.
+                                     * The post-smoothing object.
+                                     */
+    SmartPointer<const SMOOTHER> post;
+    
+                                    /**
+                                     * The coarse grid solver.
                                      */
-    std::vector<SparseMatrix<double> > prolongation_matrices;
+    SmartPointer<const COARSE> coarse;
 };
 
 
+
+
 /* --------------------------- inline functions --------------------- */
 
 
-template<int dim>
-inline
-const SparseMatrix<double>&
-Multigrid<dim>::get_matrix (const unsigned int level) const
+template <class VECTOR>
+template <int dim>
+Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
+                             const unsigned int min_level,
+                             const unsigned int max_level)
+               :
+               minlevel(min_level),
+                                    maxlevel(std::min(mg_dof_handler.get_tria().n_levels()-1,
+                                                      max_level)),
+                                    defect(minlevel,maxlevel),
+                                    solution(minlevel,maxlevel),
+                                    t(minlevel,maxlevel)
+{
+};
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+Multigrid<VECTOR>::level_mgstep(const unsigned int        level,
+                               const MGLevelObject<MATRIX>& matrix,
+                               const TRANSFER&              transfer,
+                               const SMOOTHER     &pre_smooth,
+                               const SMOOTHER     &post_smooth,
+                               const COARSE &coarse_grid_solver)
 {
-  Assert((level>=minlevel) && (level<maxlevel),
-        ExcIndexRange(level, minlevel, maxlevel));
+#ifdef MG_DEBUG
+  char *name = new char[100];
+  sprintf(name, "MG%d-defect",level);
+  print_vector(level, defect[level], name);
+#endif
+
+  solution[level] = 0.;
+  
+  if (level == minlevel)
+    {
+      coarse_grid_solver(level, solution[level], defect[level]);
+#ifdef MG_DEBUG
+      sprintf(name, "MG%d-solution",level);
+      print_vector(level, solution[level], name);
+#endif
+      return;
+    }
+
+                          // smoothing of the residual by modifying s
+  pre_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-pre",level);
+  print_vector(level, solution[level], name);
+#endif
+  
+                                  // t = -A*solution[level]
+  matrix[level].vmult(t[level], solution[level]);
+  t[level].scale(-1);
+  
+                                  // make t rhs of lower level
+                                  // The non-refined parts of the
+                                  // coarse-level defect already contain
+                                  // the global defect, the refined parts
+                                  // its restriction.
+  for (unsigned int l = level;l>minlevel;--l)
+    {
+      t[l-1] = 0.;
+      transfer.restrict_and_add (l, t[l-1], t[l]);
+      defect[l-1] += t[l-1];
+    }
+
+                                  // add additional DG contribution
+//  edge_vmult(level, defect[level-1], defect[level]);
+  
+                                  // do recursion
+  level_mgstep(level-1, matrix, transfer, pre_smooth, post_smooth, coarse_grid_solver);
+
+                                  // reset size of the auxiliary
+                                  // vector, since it has been
+                                  // resized in the recursive call to
+                                  // level_mgstep directly above
+  t[level] = 0.;
+
+                                  // do coarse grid correction
+
+  transfer.prolongate(level, t[level], solution[level-1]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-cgc",level);
+  print_vector(level, t, name);
+#endif
+
+  solution[level] += t[level];
   
-  return (*level_matrices)[level];
+                                  // post-smoothing
+  post_smooth.smooth(level, solution[level], defect[level]);
+
+#ifdef MG_DEBUG
+  sprintf(name, "MG%d-post",level);
+  print_vector(level, solution[level], name);
+
+  delete[] name;
+#endif
+}
+
+
+template <class VECTOR>
+template <class MATRIX, class TRANSFER, class SMOOTHER, class COARSE>
+void
+Multigrid<VECTOR>::vcycle(const MGLevelObject<MATRIX>& matrix,
+                         const TRANSFER&              transfer,
+                         const SMOOTHER     &pre_smooth,
+                         const SMOOTHER     &post_smooth,
+                         const COARSE &coarse_grid_solver)
+{
+                                  // The defect vector has been
+                                  // initialized by copy_to_mg. Now
+                                  // adjust the other vectors.
+  solution.resize(minlevel, maxlevel);
+  t.resize(minlevel, maxlevel);
+  
+  for (unsigned int level=minlevel; level<=maxlevel;++level)
+    {
+      solution[level].reinit(defect[level]);
+      t[level].reinit(defect[level]);
+    }
+
+  level_mgstep (maxlevel, matrix, transfer, pre_smooth, post_smooth, coarse_grid_solver);
+//  abort ();
+}
+
+
+/* ------------------------------------------------------------------- */
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>
+::PreconditionMG(const MGDoFHandler<dim>&      mg_dof_handler,
+                Multigrid<VECTOR>& mg,
+                const MGLevelObject<MATRIX>& matrix,
+                const TRANSFER&              transfer,
+                const SMOOTHER&              pre,
+                const SMOOTHER&              post,
+                const COARSE&                coarse)
+               :
+  mg_dof_handler(&mg_dof_handler),
+  multigrid(&mg),
+  level_matrices(&matrix),
+  transfer(&transfer),
+  pre(&pre),
+  post(&post),
+  coarse(&coarse)
+{}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult (
+  VECTOR& dst,
+  const VECTOR& src) const
+{
+  transfer->copy_to_mg(*mg_dof_handler,
+                      multigrid->defect,
+                      src);
+  
+  multigrid->vcycle(*level_matrices,
+                   *transfer,
+                   *pre, *post, *coarse);
+
+  transfer->copy_from_mg(*mg_dof_handler,
+                        dst,
+                        multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::vmult_add (
+  VECTOR& dst,
+  const VECTOR& src) const
+{
+  transfer->copy_to_mg(*mg_dof_handler,
+                      multigrid->defect,
+                      src);
+  
+  multigrid->vcycle(*level_matrices,
+                   *transfer,
+                   *pre, *post, *coarse);
+
+  transfer->copy_from_mg_add(*mg_dof_handler,
+                            dst,
+                            multigrid->solution);
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult (
+  VECTOR&,
+  const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
+}
+
+
+template<int dim, class MATRIX, class VECTOR, class TRANSFER, class SMOOTHER, class COARSE>
+void
+PreconditionMG<dim, MATRIX, VECTOR, TRANSFER, SMOOTHER, COARSE>::Tvmult_add (
+  VECTOR&,
+  const VECTOR&) const
+{
+  Assert(false, ExcNotImplemented());
 }
 
 
index 33484dea773dab7a4654a5687890e17dc4c6c4fa..ccf50ff2bbcdde8e1b0ba2d9fa4336e91272f80b 100644 (file)
@@ -276,7 +276,7 @@ void DoFRenumbering::Cuthill_McKee (MGDoFHandler<dim>               &dof_handler
                                   // make the connection graph
   SparsityPattern sparsity (dof_handler.n_dofs(level),
                            dof_handler.max_couplings_between_dofs());
-  MGDoFTools::make_sparsity_pattern (dof_handler, sparsity, level);
+  MGTools::make_sparsity_pattern (dof_handler, sparsity, level);
     
   const unsigned int n_dofs = sparsity.n_rows();
                                   // store the new dof numbers; invalid_dof_index means
@@ -531,6 +531,103 @@ void DoFRenumbering::component_wise (DoFHandler<dim>                 &dof_handle
 
 
 
+template <int dim>
+void DoFRenumbering::component_wise (MGDoFHandler<dim>& dof_handler,
+                                    unsigned int level,
+                                    const std::vector<unsigned int> &component_order_arg)
+{
+  const unsigned int dofs_per_cell = dof_handler.get_fe().dofs_per_cell;
+
+                                  // do nothing if the FE has only
+                                  // one component
+  if (dof_handler.get_fe().n_components() == 1)
+    return;
+  
+  std::vector<unsigned int> component_order (component_order_arg);
+  if (component_order.size() == 0)
+    for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+      component_order.push_back (i);
+
+                                  // check whether the component list has
+                                  // the right length and contains all
+                                  // component numbers
+  Assert (component_order.size() == dof_handler.get_fe().n_components(),
+         ExcInvalidComponentOrder());
+  for (unsigned int i=0; i<dof_handler.get_fe().n_components(); ++i)
+    Assert (std::find (component_order.begin(), component_order.end(), i)
+           != component_order.end(),
+           ExcInvalidComponentOrder ());
+
+                                  // vector to hold the dof indices on
+                                  // the cell we visit at a time
+  std::vector<unsigned int> local_dof_indices(dofs_per_cell);
+                                  // prebuilt list to which component
+                                  // a given dof on a cell belongs
+  std::vector<unsigned int> component_list (dofs_per_cell);
+  for (unsigned int i=0; i<component_list.size(); ++i)
+    component_list[i] = dof_handler.get_fe().system_to_component_index(i).first;
+                                  // set up a map where for each component
+                                  // the respective degrees of freedom are
+                                  // collected.
+                                  //
+                                  // note that this map is sorted by component
+                                  // but that within each component it is NOT
+                                  // sorted by dof index. note also that some
+                                  // dof indices are entered multiply, so we
+                                  // will have to take care of that
+  std::vector<std::vector<unsigned int> > component_to_dof_map (dof_handler.get_fe().n_components());
+  for (typename MGDoFHandler<dim>::cell_iterator cell=dof_handler.begin(level);
+       cell!=dof_handler.end(level); ++cell)
+    {
+                                      // on each cell: get dof indices
+                                      // and insert them into the global
+                                      // list using their component
+      cell->get_mg_dof_indices (local_dof_indices);
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       component_to_dof_map[component_list[i]].push_back (local_dof_indices[i]);
+    };
+  
+                                  // now we've got all indices sorted into
+                                  // buckets labelled with their component
+                                  // number. we've only got to traverse this
+                                  // list and assign the new indices
+                                  //
+                                  // however, we first want to sort the
+                                  // indices entered into the buckets to
+                                  // preserve the order within each component
+                                  // and during this also remove duplicate
+                                  // entries
+  for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+    {
+      std::sort (component_to_dof_map[component].begin(),
+                component_to_dof_map[component].end());
+      component_to_dof_map[component].erase (std::unique (component_to_dof_map[component].begin(),
+                                                         component_to_dof_map[component].end()),
+                                            component_to_dof_map[component].end());
+    };
+  
+  unsigned int next_free_index = 0;
+  std::vector<unsigned int> new_indices (dof_handler.n_dofs(level),
+                                        DoFHandler<dim>::invalid_dof_index);
+  for (unsigned int component=0; component<dof_handler.get_fe().n_components(); ++component)
+    {
+      const typename std::vector<unsigned int>::const_iterator
+       begin_of_component = component_to_dof_map[component].begin(),
+       end_of_component   = component_to_dof_map[component].end();
+      
+      for (typename std::vector<unsigned int>::const_iterator dof_index = begin_of_component;
+          dof_index != end_of_component; ++dof_index)
+         new_indices[*dof_index] = next_free_index++;
+    };
+
+//  Assert (next_free_index == dof_handler.n_dofs(level),
+//       ExcInternalError());
+
+  dof_handler.renumber_dofs (level, new_indices);
+};
+
+
+
 template <int dim>
 void
 DoFRenumbering::sort_selected_dofs_back (DoFHandler<dim>         &dof_handler,
@@ -775,6 +872,11 @@ template
 void DoFRenumbering::component_wise (DoFHandler<deal_II_dimension>&,
                                     const std::vector<unsigned int>&);
 
+template
+void DoFRenumbering::component_wise (MGDoFHandler<deal_II_dimension>&,
+                                    unsigned int,
+                                    const std::vector<unsigned int>&);
+
 
 template
 void DoFRenumbering::cell_wise_dg (DoFHandler<deal_II_dimension>&,
index 8b50079c181b7ff9acca517c896a054a4d172454..aa8d42fc470b441ef5b4bc0371cbdb8b47c95a2c 100644 (file)
 
 #include <multigrid/mg_base.h>
 #include <multigrid/mg_smoother.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 #include <iostream>
 #include <cmath>
 
 
-MGBase::~MGBase ()
-{}
-
-
-MGBase::MGBase(const MGTransferBase &transfer,
-              const unsigned        minlevel,
-              const unsigned        maxlevel)
-               :
-               maxlevel(maxlevel),
-               minlevel(minlevel),
-               defect(minlevel,maxlevel),
-               solution(minlevel,maxlevel),
-               t(minlevel,maxlevel),
-               transfer(&transfer)
-{
-  Assert(minlevel <= maxlevel,
-        ExcSwitchedLevels(minlevel, maxlevel));
-}
-
-
-void
-MGBase::vcycle(const MGSmootherBase     &pre_smooth,
-              const MGSmootherBase     &post_smooth,
-              const MGCoarseGridSolver &coarse_grid_solver)
-{
-  level_mgstep (maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-//  abort ();
-}
-
-
-void
-MGBase::level_mgstep(const unsigned int        level,
-                    const MGSmootherBase     &pre_smooth,
-                    const MGSmootherBase     &post_smooth,
-                    const MGCoarseGridSolver &coarse_grid_solver)
-{
-#ifdef MG_DEBUG
-  char *name = new char[100];
-  sprintf(name, "MG%d-defect",level);
-  print_vector(level, defect[level], name);
-#endif
-
-  solution[level] = 0.;
-  
-  if (level == minlevel)
-    {
-      coarse_grid_solver(level, solution[level], defect[level]);
-#ifdef MG_DEBUG
-      sprintf(name, "MG%d-solution",level);
-      print_vector(level, solution[level], name);
-#endif
-      return;
-    }
-
-                          // smoothing of the residual by modifying s
-  pre_smooth.smooth(level, solution[level], defect[level]);
-
-#ifdef MG_DEBUG
-  sprintf(name, "MG%d-pre",level);
-  print_vector(level, solution[level], name);
-#endif
-  
-  t[level] = 0.;
-
-                                  // t = -A*solution[level]
-  level_vmult(level, t[level], solution[level], defect[level]);
-  
-                                  // make t rhs of lower level
-                                  // The non-refined parts of the
-                                  // coarse-level defect already contain
-                                  // the global defect, the refined parts
-                                  // its restriction.
-  for (unsigned int l = level;l>minlevel;--l)
-    {
-      t[l-1] = 0.;
-      transfer->restrict_and_add (l, t[l-1], t[l]);
-      defect[l-1] += t[l-1];
-    }
-
-                                  // add additional DG contribution
-//  edge_vmult(level, defect[level-1], defect[level]);
-  
-                                  // do recursion
-  level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
-
-                                  // reset size of the auxiliary
-                                  // vector, since it has been
-                                  // resized in the recursive call to
-                                  // level_mgstep directly above
-  t[level] = 0.;
-
-                                  // do coarse grid correction
-
-  transfer->prolongate(level, t[level], solution[level-1]);
-
-#ifdef MG_DEBUG
-  sprintf(name, "MG%d-cgc",level);
-  print_vector(level, t, name);
-#endif
-
-  solution[level] += t[level];
-  
-                                  // post-smoothing
-  post_smooth.smooth(level, solution[level], defect[level]);
-
-#ifdef MG_DEBUG
-  sprintf(name, "MG%d-post",level);
-  print_vector(level, solution[level], name);
-
-  delete[] name;
-#endif
-}
-
-
-void
-MGBase::edge_vmult (const unsigned int,
-                   Vector<double>&,
-                   const Vector<double>&)
-{}
-
 //////////////////////////////////////////////////////////////////////
 
-MGCoarseGridSolver::~MGCoarseGridSolver()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
+//  template <class VECTOR>
+//  MGTransferBase<VECTOR>::~MGTransferBase()
+//  {};
 
-MGTransferBase::~MGTransferBase()
-{};
 
+//  // Explicit instantiations
 
+//  template class MGTransferBase<Vector<double> >;
+//  template class MGTransferBase<Vector<float> >;
+//  template class MGTransferBase<BlockVector<double> >;
+//  template class MGTransferBase<BlockVector<float> >;
index 671bb3fb308d52809221c122cfa4be13e656918f..e84883e623c27ffba4b7a935f98db5200989b3b4 100644 (file)
 //----------------------------  mg_dof_tools.cc  ---------------------------
 
 
+#include <base/multithread_info.h>
+#include <base/thread_management.h>
 #include <lac/sparsity_pattern.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
+#include <lac/block_sparsity_pattern.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/sparsity_pattern.h>
+#include <lac/block_vector.h>
 #include <grid/tria.h>
 #include <grid/tria_iterator.h>
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
 #include <multigrid/mg_dof_tools.h>
+#include <multigrid/mg_base.h>
+#include <dofs/dof_tools.h>
 #include <fe/fe.h>
 
+//#include <multigrid/mg_tools.templates.h>
 
-template <int dim>
-void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<dim> &dof,
+#include <algorithm>
+#include <numeric>
+
+template <int dim, class SparsityPattern>
+void MGTools::make_sparsity_pattern (const MGDoFHandler<dim> &dof,
                                        SparsityPattern         &sparsity,
                                        const unsigned int       level)
 {
@@ -50,9 +62,9 @@ void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<dim> &dof,
 
 
 
-template<int dim>
+template<int dim, class SparsityPattern>
 void
-MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
                                        SparsityPattern       &sparsity,
                                        const unsigned int level)
 {
@@ -106,9 +118,9 @@ MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
 
 
 
-template<int dim>
+template<int dim, class SparsityPattern>
 void
-MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
                                             SparsityPattern       &sparsity,
                                             const unsigned int level)
 {
@@ -164,9 +176,10 @@ MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<dim> &dof,
 
 
 
-template<int dim>
+//TODO: Replace FullMatrix by vector2d<bool>
+template<int dim, class SparsityPattern>
 void
-MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
                                        SparsityPattern       &sparsity,
                                        const unsigned int level,
                                        const FullMatrix<double>& int_mask,
@@ -262,23 +275,281 @@ MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<dim> &dof,
 
 
 
+template <int dim>
+void
+MGTools::count_dofs_per_component (const MGDoFHandler<dim>& dof_handler,
+                                     std::vector<std::vector<unsigned int> >& result)
+{
+  const unsigned int nlevels = dof_handler.get_tria().n_levels();
+  
+  Assert (result.size() == nlevels,
+         ExcDimensionMismatch(result.size(), nlevels));
 
-// explicit instantiations
-template void MGDoFTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
-                                                SparsityPattern &,
-                                                const unsigned int);
-template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
-                                                     SparsityPattern &,
-                                                     const unsigned int);
+  const unsigned int n_components = dof_handler.get_fe().n_components();
+  for (unsigned int l=0;l<nlevels;++l)
+    {
+      result[l].resize (n_components);
+      
+                                      // special case for only one
+                                      // component. treat this first
+                                      // since it does not require any
+                                      // computations
+      if (n_components == 1)
+       {
+         result[l][0] = dof_handler.n_dofs(l);
+       } else {
+                                          // 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::ThreadManager thread_manager;
+         for (unsigned int i=0; i<n_components; ++i)
+           {
+             void (*fun_ptr) (const unsigned int       level,
+                              const MGDoFHandler<dim>    &,
+                              const std::vector<bool>    &,
+                              std::vector<bool>          &)
+               = &DoFTools::template extract_level_dofs<dim>;
+             component_select[i][i] = true;
+             Threads::spawn (thread_manager,
+                             Threads::encapsulate (fun_ptr)
+                             .collect_args (l, dof_handler, component_select[i],
+                                            dofs_in_component[i]));
+           };
+         thread_manager.wait();
+         
+                                          // next count what we got
+         for (unsigned int i=0; i<n_components; ++i)
+             result[l][i] = std::count(dofs_in_component[i].begin(),
+                                       dofs_in_component[i].end(),
+                                       true);
+         
+                                          // finally sanity check
+         Assert (std::accumulate (result[l].begin(),
+                                  result[l].end(), 0U)
+                 ==
+                 dof_handler.n_dofs(l),
+                 ExcInternalError());
+       }
+    }
+}
+
+
+
+template<int dim, typename number>
+void
+MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                          MGLevelObject<Vector<number> >& v)
+{
+  for (unsigned int level=v.get_minlevel();
+       level<=v.get_maxlevel();++level)
+    {
+      unsigned int n = mg_dof.n_dofs (level);
+      v[level].reinit(n);
+    }
+
+}
+
+
+
+
+
+template<int dim, typename number>
+void
+MGTools::reinit_vector (const MGDoFHandler<dim>& mg_dof,
+                          MGLevelObject<BlockVector<number> >& v,
+                          const std::vector<bool>& selected)
+{
+  const unsigned int ncomp = mg_dof.get_fe().n_components();
+  
+  if (selected.size() == 0)
+    {
+      selected.resize(ncomp);
+      fill_n (selected.begin(), ncomp, true);
+    }
+
+  Assert (selected.size() == ncomp,
+         ExcDimensionMismatch(selected.size(), ncomp));
 
-template void MGDoFTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
-                                                          SparsityPattern &,
-                                                          const unsigned int);
+  unsigned int n_selected = std::accumulate (selected.begin(),
+                                            selected.end(),
+                                            0U);
+                            
+  std::vector<std::vector<unsigned int> >
+    ndofs(mg_dof.n_levels(), ncomp);
+
+  count_dofs_per_component (mg_dof, ndofs);
+  
+  for (unsigned int level=v.get_minlevel();
+       level<=v.get_maxlevel();++level)
+    {
+      v[level].reinit(n_selected);
+      unsigned int k=0;
+      for (unsigned int i=0;i<ncomp;++i)
+       if (selected[i])
+         v[level].block(k++).reinit(ndofs[level][i]);
+    }
+}
+
+
+
+
+
+// explicit instantiations
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                               SparsityPattern &,
+                               const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    SparsityPattern &,
+                                    const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+                                         SparsityPattern &,
+                                         const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                               CompressedSparsityPattern &,
+                               const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    CompressedSparsityPattern &,
+                                    const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+                                         CompressedSparsityPattern &,
+                                         const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                               BlockSparsityPattern &,
+                               const unsigned int);
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    BlockSparsityPattern &,
+                                    const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+                                         BlockSparsityPattern &,
+                                         const unsigned int);
+
+template void
+MGTools::make_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                               CompressedBlockSparsityPattern &,
+                               const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    CompressedBlockSparsityPattern &,
+                                    const unsigned int);
+
+template void
+MGTools::make_flux_sparsity_pattern_edge (const MGDoFHandler<deal_II_dimension> &,
+                                         CompressedBlockSparsityPattern &,
+                                         const unsigned int);
 
 #if deal_II_dimension > 1
-template void MGDoFTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
-                                                     SparsityPattern &,
-                                                     const unsigned int,
-                                                     const FullMatrix<double>&,
-                                                     const FullMatrix<double>&);
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    SparsityPattern &,
+                                    const unsigned int,
+                                    const FullMatrix<double>&,
+                                    const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    CompressedSparsityPattern &,
+                                    const unsigned int,
+                                    const FullMatrix<double>&,
+                                    const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    BlockSparsityPattern &,
+                                    const unsigned int,
+                                    const FullMatrix<double>&,
+                                    const FullMatrix<double>&);
+
+template void
+MGTools::make_flux_sparsity_pattern (const MGDoFHandler<deal_II_dimension> &,
+                                    CompressedBlockSparsityPattern &,
+                                    const unsigned int,
+                                    const FullMatrix<double>&,
+                                    const FullMatrix<double>&);
+
 #endif
+
+template void MGTools::reinit_vector (const MGDoFHandler<deal_II_dimension>&,
+                                     MGLevelObject<Vector<double> >&);
+template void MGTools::reinit_vector (const MGDoFHandler<deal_II_dimension>&,
+                                     MGLevelObject<Vector<float> >&);
+
+
+//  template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+//                                const MGTransferBase<Vector<double> >&,
+//                                MGLevelObject<Vector<double> >&,
+//                                const Vector<double>&);
+
+//  template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+//                                 Vector<double>&,
+//                                 const MGLevelObject<Vector<double> >&);
+
+//  template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+//                                     Vector<double>&,
+//                                     const MGLevelObject<Vector<double> >&);
+
+//  template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+//                                const MGTransferBase<Vector<double> >&,
+//                                MGLevelObject<Vector<double> >&,
+//                                const Vector<float>&);
+
+//  template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+//                                 Vector<double>&,
+//                                 const MGLevelObject<Vector<float> >&);
+
+//  template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+//                                     Vector<double>&,
+//                                     const MGLevelObject<Vector<float> >&);
+
+//  template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+//                                const MGTransferBase<Vector<double> >&,
+//                                MGLevelObject<Vector<double> >&,
+//                                const BlockVector<double>&);
+
+//  template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+//                                 BlockVector<double>&,
+//                                 const MGLevelObject<Vector<double> >&);
+
+//  template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+//                                     BlockVector<double>&,
+//                                     const MGLevelObject<Vector<double> >&);
+
+//  template void MGTools::copy_to_mg (const MGDoFHandler<deal_II_dimension>&,
+//                                const MGTransferBase<Vector<double> >&,
+//                                MGLevelObject<Vector<double> >&,
+//                                const BlockVector<float>&);
+
+//  template void MGTools::copy_from_mg(const MGDoFHandler<deal_II_dimension>&,
+//                                 BlockVector<double>&,
+//                                 const MGLevelObject<Vector<float> >&);
+
+//  template void MGTools::copy_from_mg_add(const MGDoFHandler<deal_II_dimension>&,
+//                                     BlockVector<double>&,
+//                                     const MGLevelObject<Vector<float> >&);
+
+template void MGTools::count_dofs_per_component (const MGDoFHandler<deal_II_dimension>&,
+                                                std::vector<std::vector<unsigned int> >&);
+
index 9fd6e980537fe1c0df7969f96189362db13d519b..fb505e35e19abcac427230d8f2cadbe9f99fd9f4 100644 (file)
 #include <multigrid/mg_dof_accessor.h>
 #include <grid/tria_iterator.h>
 #include <fe/fe.h>
-#include <multigrid/multigrid.h>
+//#include <multigrid/multigrid.h>
 #include <multigrid/mg_smoother.h>
 #include <multigrid/mg_smoother.templates.h>
 #include <lac/vector.h>
+#include <lac/block_vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 
 #include <algorithm>
 
 //////////////////////////////////////////////////////////////////////
 
 
-MGSmootherBase::~MGSmootherBase()
-{};
-
-
-//////////////////////////////////////////////////////////////////////
-
-
-void
-MGSmootherIdentity::smooth (const unsigned int,
-                           Vector<double>       &,
-                           const Vector<double> &) const
-{};
-
 
 #if deal_II_dimension == 1
 
-MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps)
+MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<1> &/*mg_dof*/,
+                                           unsigned int steps)
                :
                steps(steps)
 {
@@ -58,7 +49,8 @@ MGSmoother::MGSmoother (const MGDoFHandler<1> &/*mg_dof*/, unsigned int steps)
 #if deal_II_dimension > 1
 
 template <int dim>
-MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
+MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<dim> &mg_dof,
+                                           unsigned int steps)
                :
                steps(steps)
 {
@@ -121,9 +113,10 @@ MGSmoother::MGSmoother (const MGDoFHandler<dim> &mg_dof, unsigned int steps)
 
 //////////////////////////////////////////////////////////////////////
 
+template <class VECTOR>
 void
-MGSmoother::set_zero_interior_boundary (const unsigned int  level,
-                                       Vector<double>      &u) const
+MGSmootherContinuous::set_zero_interior_boundary (const unsigned int level,
+                                                 VECTOR&            u) const
 {
   if (level==0)
     return;
@@ -140,22 +133,89 @@ MGSmoother::set_zero_interior_boundary (const unsigned int  level,
 // don't do the following instantiation in 1d, since there is a specialized
 // function there
 #if deal_II_dimension > 1
-template MGSmoother::MGSmoother (const MGDoFHandler<deal_II_dimension>&, unsigned int);
+template MGSmootherContinuous::MGSmootherContinuous (const MGDoFHandler<deal_II_dimension>&, unsigned int);
 #endif
 
 template
-MGSmootherRelaxation<float>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>        &mg_dof,
-                      const MGLevelObject<SparseMatrix<float> > &matrix,
-                      const function_ptr                            relaxation,
-                      const unsigned int                            steps,
-                      const double                                  omega);
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+                                            Vector<double>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+                                            Vector<float>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+                                            BlockVector<double>&) const;
+
+template
+void MGSmootherContinuous::set_zero_interior_boundary (const unsigned int,
+                                            BlockVector<float>&) const;
+
+
+  
+template
+MGSmootherRelaxation<SparseMatrix<float>, Vector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<SparseMatrix<float> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<float>, Vector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<SparseMatrix<float> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<double>, Vector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<SparseMatrix<double> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+template
+MGSmootherRelaxation<SparseMatrix<double>, Vector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<SparseMatrix<double> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<float>, BlockVector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<BlockSparseMatrix<float> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<float>, BlockVector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<BlockSparseMatrix<float> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
+
+template
+MGSmootherRelaxation<BlockSparseMatrix<double>, BlockVector<float> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<BlockSparseMatrix<double> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
 
 template
-MGSmootherRelaxation<double>
-::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>         &mg_dof,
-                      const MGLevelObject<SparseMatrix<double> > &matrix,
-                      const function_ptr                             relaxation,
-                      const unsigned int                             steps,
-                      const double                                   omega);
+MGSmootherRelaxation<BlockSparseMatrix<double>, BlockVector<double> >
+::MGSmootherRelaxation(const MGDoFHandler<deal_II_dimension>&,
+                      const MGLevelObject<BlockSparseMatrix<double> >&,
+                      const function_ptr,
+                      const unsigned int,
+                      const double);
 
index be32a19fa508b975963303c482ffc322a1442991..879b8e7c32a4827adcfc71e7d4f2a3ae9fc23255 100644 (file)
 
 
 
-#include <grid/tria.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
-#include <grid/tria_iterator.h>
-#include <fe/fe.h>
-#include <multigrid/multigrid.h>
-#include <multigrid/multigrid.templates.h>
-#include <multigrid/mg_smoother.h>
 #include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
+#include <multigrid/mg_transfer.h>
 
 
 MGTransferPrebuilt::~MGTransferPrebuilt () 
@@ -50,3 +45,74 @@ void MGTransferPrebuilt::restrict_and_add (const unsigned int   from_level,
 };
 
 
+
+
+
+MGTransferBlock::~MGTransferBlock () 
+{};
+
+
+void MGTransferBlock::prolongate (const unsigned int   to_level,
+                                    BlockVector<double>       &dst,
+                                    const BlockVector<double> &src) const 
+{
+  Assert ((to_level >= 1) && (to_level<=prolongation_matrices.size()),
+         ExcIndexRange (to_level, 1, prolongation_matrices.size()+1));
+
+  unsigned int k=0;
+  for (unsigned int b=0; b<src.n_blocks();++b)
+    {
+      if (!selected[k])
+       ++k;
+      prolongation_matrices[to_level-1].block(k,k).vmult (dst.block(b), src.block(b));
+      ++k;
+    }
+};
+
+
+void MGTransferBlock::restrict_and_add (const unsigned int   from_level,
+                                          BlockVector<double>       &dst,
+                                          const BlockVector<double> &src) const 
+{
+  Assert ((from_level >= 1) && (from_level<=prolongation_matrices.size()),
+         ExcIndexRange (from_level, 1, prolongation_matrices.size()+1));
+
+  unsigned int k=0;
+  for (unsigned int b=0; b<src.n_blocks();++b)
+    {
+      if (!selected[k])
+       ++k;
+      prolongation_matrices[from_level-1].block(k,k).Tvmult_add (dst.block(b), src.block(b));
+      ++k;
+    }
+};
+
+
+
+
+MGTransferSelect::~MGTransferSelect () 
+{};
+
+
+void MGTransferSelect::prolongate (const unsigned int   to_level,
+                                    Vector<double>       &dst,
+                                    const Vector<double> &src) const 
+{
+  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)
+       .vmult (dst, src);
+};
+
+
+void MGTransferSelect::restrict_and_add (const unsigned int   from_level,
+                                          Vector<double>       &dst,
+                                          const Vector<double> &src) const 
+{
+  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)
+    .Tvmult_add (dst, src);
+};
index 6861fcd10012edfce0096f75406c4139060ca0c1..8eeb927ee41867995b6163f4ff74250b71e93027 100644 (file)
@@ -1,4 +1,4 @@
-//----------------------------  multigrid.cc  ---------------------------
+//-----------------------------------------------------------------------
 //    $Id$
 //    Version: $Name$
 //
@@ -9,52 +9,20 @@
 //    to the file deal.II/doc/license.html for the  text  and
 //    further information on this license.
 //
-//----------------------------  multigrid.cc  ---------------------------
+//-----------------------------------------------------------------------
 
 
+#include <lac/vector.h>
+#include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 #include <grid/tria.h>
-#include <dofs/dof_constraints.h>
-#include <multigrid/mg_dof_handler.h>
-#include <multigrid/mg_dof_accessor.h>
 #include <grid/tria_iterator.h>
+#include <dofs/dof_tools.h>
 #include <fe/fe.h>
-#include <multigrid/multigrid.h>
-#include <multigrid/multigrid.templates.h>
-#include <multigrid/mg_smoother.h>
-#include <lac/vector.h>
-
-
-/* --------------------------------- MG ------------------------------------ */
-
-template <int dim>
-Multigrid<dim>::Multigrid (const MGDoFHandler<dim>                       &mg_dof_handler,
-                          const ConstraintMatrix                        &constraints,
-                          const MGLevelObject<SparsityPattern>       &level_sparsities,
-                          const MGLevelObject<SparseMatrix<double> > &level_matrices,
-                          const MGTransferBase                          &transfer,
-                          const unsigned int                             minlevel,
-                          const unsigned int                             maxlevel)
-               :
-               MGBase(transfer, minlevel,
-                      std::min(mg_dof_handler.get_tria().n_levels()-1,
-                          maxlevel)),
-               mg_dof_handler(&mg_dof_handler),
-               level_sparsities(&level_sparsities),
-               level_matrices(&level_matrices),
-               constraints(&constraints)
-{};
-
-
-template <int dim>
-void
-Multigrid<dim>::level_vmult (const unsigned int    level,
-                            Vector<double>       &result,
-                            const Vector<double> &u,
-                            const Vector<double> &/* rhs */)
-{
-  (*level_matrices)[level].vmult(result,u);
-  result.scale(-1.);
-}
+#include <multigrid/mg_dof_handler.h>
+#include <multigrid/mg_dof_accessor.h>
+#include <multigrid/mg_transfer.h>
+#include <multigrid/mg_dof_tools.h>
 
 
 /* ----------------------------- MGTransferPrebuilt ------------------------ */
@@ -66,18 +34,20 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
   const unsigned int n_levels      = mg_dof.get_tria().n_levels();
   const unsigned int dofs_per_cell = mg_dof.get_fe().dofs_per_cell;
   
+  sizes.resize(n_levels);
+  for (unsigned int l=0;l<n_levels;++l)
+    sizes[l] = mg_dof.n_dofs(l);
+  
                                   // reset the size of the array of
                                   // matrices
-  prolongation_sparsities.clear ();
-  prolongation_matrices.clear ();
-  prolongation_sparsities.reserve (n_levels-1);
-  prolongation_matrices.reserve (n_levels-1);
+  prolongation_sparsities.resize (n_levels-1);
+  prolongation_matrices.resize (n_levels-1);
 
 
-// two fields which will store the
+                                  // two fields which will store the
                                   // indices of the multigrid dofs
                                   // for a cell and one of its children
-  std::vector<unsigned int> dof_indices_mother (dofs_per_cell);
+  std::vector<unsigned int> dof_indices_parent (dofs_per_cell);
   std::vector<unsigned int> dof_indices_child (dofs_per_cell);
   
                                   // for each level: first build the sparsity
@@ -87,24 +57,23 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
                                   // level which have children
   for (unsigned int level=0; level<n_levels-1; ++level)
     {
-      prolongation_sparsities.push_back (SparsityPattern());
                                       // reset the dimension of the structure.
                                       // note that for the number of entries
-                                      // per row, the number of mother dofs
+                                      // per row, the number of parent dofs
                                       // coupling to a child dof is
                                       // necessary. this, of course, is the
                                       // number of degrees of freedom per
                                       // cell
-      prolongation_sparsities.back().reinit (mg_dof.n_dofs(level+1),
-                                            mg_dof.n_dofs(level),
+      prolongation_sparsities[level].reinit (sizes[level+1],
+                                             sizes[level],
 //TODO:[WB,GK] evil hack, must be corrected!
-                                            dofs_per_cell+1);
+                                             dofs_per_cell+1);
       
       for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
           cell != mg_dof.end(level); ++cell)
        if (cell->has_children())
          {
-           cell->get_mg_dof_indices (dof_indices_mother);
+           cell->get_mg_dof_indices (dof_indices_parent);
 
            for (unsigned int child=0;
                 child<GeometryInfo<dim>::children_per_cell; ++child)
@@ -119,20 +88,19 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
 
                                                 // now tag the entries in the
                                                 // matrix which will be used
-                                                // for this pair of mother/child
+                                                // for this pair of parent/child
                for (unsigned int i=0; i<dofs_per_cell; ++i)
                  for (unsigned int j=0; j<dofs_per_cell; ++j)
                    if (prolongation(i,j) != 0)
                      {
                        prolongation_sparsities[level].add
                          (dof_indices_child[i],
-                          dof_indices_mother[j]);
+                          dof_indices_parent[j]);
                      };
              };
          };
       prolongation_sparsities[level].compress ();
 
-      prolongation_matrices.push_back (SparseMatrix<double>());
       prolongation_matrices[level].reinit (prolongation_sparsities[level]);
 
                                       // now actually build the matrices
@@ -140,7 +108,7 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
           cell != mg_dof.end(level); ++cell)
        if (cell->has_children())
          {
-           cell->get_mg_dof_indices (dof_indices_mother);
+           cell->get_mg_dof_indices (dof_indices_parent);
 
            for (unsigned int child=0;
                 child<GeometryInfo<dim>::children_per_cell; ++child)
@@ -161,7 +129,7 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
                      {
                        prolongation_matrices[level].set
                          (dof_indices_child[i],
-                          dof_indices_mother[j],
+                          dof_indices_parent[j],
                           prolongation(i,j));
                      };
              };
@@ -170,21 +138,205 @@ void MGTransferPrebuilt::build_matrices (const MGDoFHandler<dim> &mg_dof)
 };
 
 
-// explicit instatations
-template class Multigrid<deal_II_dimension>;
-template
-void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
+/* ----------------------------- MGTransferBlock ------------------------ */
 
-template
-void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<float>& src);
 
-template
-void Multigrid<deal_II_dimension>::copy_to_mg (const Vector<double>& src);
+template <int dim>
+void MGTransferBlockBase::build_matrices (const MGDoFHandler<dim> &mg_dof,
+                                     std::vector<bool> select) 
+{
+  const FiniteElement<dim>& fe = mg_dof.get_fe();
+  const unsigned int n_components  = fe.n_components();
+  const unsigned int dofs_per_cell = fe.dofs_per_cell;  
+  const unsigned int n_levels      = mg_dof.get_tria().n_levels();
+  
+  selected = select;
+  
+  Assert (selected.size() == n_components,
+         ExcDimensionMismatch(selected.size(), n_components))
 
-template
-void Multigrid<deal_II_dimension>::copy_from_mg (Vector<float>& src) const;
+  sizes.resize(n_levels);
+  MGTools::count_dofs_per_component(mg_dof, sizes);
+  
+                                  // reset the size of the array of
+                                  // matrices
+  prolongation_sparsities.resize (n_levels-1);
+  prolongation_matrices.resize (n_levels-1);
+
+                                  // two fields which will store the
+                                  // indices of the multigrid dofs
+                                  // for a cell and one of its children
+  std::vector<unsigned int> dof_indices_parent (dofs_per_cell);
+  std::vector<unsigned int> dof_indices_child (dofs_per_cell);
+  
+                                  // for each level: first build the
+                                  // sparsity pattern of the matrices
+                                  // and then build the matrices
+                                  // themselves. note that we only
+                                  // need to take care of cells on
+                                  // the coarser level which have
+                                  // children
+  for (unsigned int level=0; level<n_levels-1; ++level)
+    {
+                                      // reset the dimension of the
+                                      // structure.  note that for
+                                      // the number of entries per
+                                      // row, the number of parent
+                                      // dofs coupling to a child dof
+                                      // is necessary. this, is the
+                                      // number of degrees of freedom
+                                      // per cell
+      prolongation_sparsities[level].reinit (n_components, n_components);
+      for (unsigned int i=0;i<n_components;++i)
+       for (unsigned int j=0;j<n_components;++j)
+         if (i==j)
+           prolongation_sparsities[level].block(i,j).reinit(
+             sizes[level+1][i],
+             sizes[level][j],
+//TODO:[GK] Split this by component to save memory
+             dofs_per_cell);
+         else
+           prolongation_sparsities[level].block(i,j).reinit(
+             sizes[level+1][i],
+             sizes[level][j],
+             0);
+
+      prolongation_sparsities[level].collect_sizes();
+      
+      for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
+          cell != mg_dof.end(level); ++cell)
+       if (cell->has_children())
+         {
+           cell->get_mg_dof_indices (dof_indices_parent);
+
+           for (unsigned int child=0;
+                child<GeometryInfo<dim>::children_per_cell; ++child)
+             {
+                                                // set an alias to the
+                                                // prolongation matrix for
+                                                // this child
+               const FullMatrix<double> &prolongation
+                 = mg_dof.get_fe().prolongate(child);
+           
+               cell->child(child)->get_mg_dof_indices (dof_indices_child);
+
+                                                // now tag the entries in the
+                                                // matrix which will be used
+                                                // for this pair of parent/child
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   if (prolongation(i,j) != 0)
+                     {
+                       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[icomp])
+                         prolongation_sparsities[level].add
+                           (dof_indices_child[i],
+                            dof_indices_parent[j]);
+                     };
+             };
+         };
+      prolongation_sparsities[level].compress ();
+
+      prolongation_matrices[level].reinit (prolongation_sparsities[level]);
+
+                                      // now actually build the matrices
+      for (typename MGDoFHandler<dim>::cell_iterator cell=mg_dof.begin(level);
+          cell != mg_dof.end(level); ++cell)
+       if (cell->has_children())
+         {
+           cell->get_mg_dof_indices (dof_indices_parent);
+
+           for (unsigned int child=0;
+                child<GeometryInfo<dim>::children_per_cell; ++child)
+             {
+                                                // set an alias to the
+                                                // prolongation matrix for
+                                                // this child
+               const FullMatrix<double> &prolongation
+                 = mg_dof.get_fe().prolongate(child);
+           
+               cell->child(child)->get_mg_dof_indices (dof_indices_child);
+
+                                                // now set the entries in the
+                                                // matrix
+               for (unsigned int i=0; i<dofs_per_cell; ++i)
+                 for (unsigned int j=0; j<dofs_per_cell; ++j)
+                   if (prolongation(i,j) != 0)
+                     {
+                       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[icomp])
+                         prolongation_matrices[level].set
+                           (dof_indices_child[i],
+                            dof_indices_parent[j],
+                            prolongation(i,j));
+                     };
+             };
+         };
+    };
+                                  // 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;
+    }
+  
+};
+
+
+template <int dim>
+void MGTransferBlock::build_matrices (const MGDoFHandler<dim> &mg_dof,
+                                         std::vector<bool> select) 
+{
+  if (select.size() == 0)
+    select = std::vector<bool> (mg_dof.get_fe().n_components(), true);
+
+  MGTransferBlockBase::build_matrices (mg_dof, select);
+}
+
+
+template <int dim>
+void MGTransferSelect::build_matrices (const MGDoFHandler<dim> &mg_dof,
+                                        unsigned int select)
+{
+  selected = select;
+  std::vector<bool> s(mg_dof.get_fe().n_components(), false);
+  s[select] = true;
+
+  MGTransferBlockBase::build_matrices (mg_dof, s);
+}
+
+
+
+// explicit instatations
 
 template
-void Multigrid<deal_II_dimension>::copy_from_mg (Vector<double>& src) const;
+void MGTransferPrebuilt::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof);
 
+template
+void MGTransferBlock::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof,
+                                     std::vector<bool>);
 
+template
+void MGTransferSelect::build_matrices (const MGDoFHandler<deal_II_dimension> &mg_dof,
+                                      unsigned int);
index 453a58ca1ed8673f28be77bc172d92238dfc2ce7..fa4c7cf818d0d7b24db31377d666ca785bc0fcba 100644 (file)
@@ -135,7 +135,7 @@ int main()
            {
              mgstruct[i].reinit(mgdof.n_dofs(i), mgdof.n_dofs(i),
                                 mgdof.max_couplings_between_dofs());
-             MGDoFTools::make_sparsity_pattern(mgdof, mgstruct[i], i);
+             MGTools::make_sparsity_pattern(mgdof, mgstruct[i], i);
              mgstruct[i].compress();
              
              mgA[i].reinit(mgstruct[i]);

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.