]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Multigrid prepared for W- and F-cycles
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 8 Mar 2005 15:42:11 +0000 (15:42 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Tue, 8 Mar 2005 15:42:11 +0000 (15:42 +0000)
Put multigrid classes in multigrid module

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

deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/mg_coarse.h
deal.II/deal.II/include/multigrid/mg_dof_accessor.h
deal.II/deal.II/include/multigrid/mg_dof_handler.h
deal.II/deal.II/include/multigrid/mg_dof_tools.h
deal.II/deal.II/include/multigrid/mg_level_object.h
deal.II/deal.II/include/multigrid/mg_matrix.h
deal.II/deal.II/include/multigrid/mg_smoother.h
deal.II/deal.II/include/multigrid/mg_transfer.h
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h

index 80bcedc2bd740359a323d8e0dee2850dc5fdf8a9..a771e5908945e842c671ff1fc868209541b8f5c3 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -24,7 +24,8 @@
 #include <lac/vector.h>
 
 
-
+/*!@addtogroup mg */
+/*@{*/
 
 
 /**
@@ -226,9 +227,7 @@ class MGSmootherBase : public Subscriptor
                       const VECTOR&      rhs) const = 0;  
 };
 
+/*@}*/
 
 
-/*----------------------------   mgbase.h     ---------------------------*/
-
 #endif
-/*----------------------------   mgbase.h     ---------------------------*/
index 10a8e4d4d475ecf658c65eff78afcc0f6fd2c01a..8804b74b3ad3eafe57dd6c0e1b18de26aa4297aa 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -17,6 +17,8 @@
 #include <lac/full_matrix.h>
 #include <multigrid/mg_base.h>
 
+/*!@addtogroup mg */
+/*@{*/
 
 /**
  * Coarse grid solver using LAC iterative methods.
@@ -158,8 +160,9 @@ class MGCoarseGridHouseholder : public MGCoarseGridBase<VECTOR>
     mutable VECTOR aux;
 };
 
+/*@}*/
   
-
+/// @if NoDoc
 /* ------------------ Functions for MGCoarseGridLACIteration ------------ */
 
 
@@ -273,5 +276,6 @@ MGCoarseGridHouseholder<number, VECTOR>::operator() (
   work.least_squares(dst, aux);
 }
 
+/// @endif
 
 #endif
index 912df745573728faf7d4149fe37f8b1dcad4fb81..8ce17013e7af4d790ed94ccc54a816c497904dad 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -25,6 +25,8 @@ template <int dim>              class MGDoFObjectAccessor<1, dim>;
 template <int dim>              class MGDoFObjectAccessor<2, dim>;
 template <int dim>              class MGDoFObjectAccessor<3, dim>;
 
+/*!@addtogroup mg */
+/*@{*/
 
 /**
  * Define the basis for accessors to the degrees of freedom for
@@ -681,6 +683,8 @@ class MGDoFCellAccessor :  public MGDoFObjectAccessor<dim, dim> {
     DeclException0 (ExcNotUsefulForThisDimension);
 };
 
+/*@}*/
+
 
 template<> MGDoFCellAccessor<1>::face_iterator
 MGDoFCellAccessor<1>::face (const unsigned int i) const;
index 0904913095d1f496006e8265f7f7eecf24e4ad8e..2ff32ac04f339dc5f8be67ebff5885666e997ec8 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -14,9 +14,6 @@
 #define __deal2__mg_dof_handler_h
 
 
-/*----------------------------   mg_dof.h     ---------------------------*/
-
-
 #include <base/config.h>
 #include <dofs/dof_handler.h>
 
@@ -27,9 +24,10 @@ template <int dim>              class MGDoFObjectAccessor<1, dim>;
 template <int dim>              class MGDoFObjectAccessor<2, dim>;
 template <int dim>              class MGDoFObjectAccessor<3, dim>;
 
+/*!@addtogroup mg */
+/*@{*/
+
 /**
- * @brief Iterators for MDDoFHandler
- *
  * Depending on the dimension, this class defines the iterator types for
  * MGDoFHandler objects.
  *
@@ -1165,6 +1163,7 @@ class MGDoFHandler : public DoFHandler<dim>
     template <int dim1, int dim2> friend class MGDoFObjectAccessor;
 };
 
+/*@}*/
 
 /* ----------------------- Inline functions of MGDoFHandler -------------------*/
 
index 7d0a3f894eb9a39186e8197b1fd0e7244c9d0af1..746b3ac80d0530e985a0a2b71e7eb6da26075b05 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -24,6 +24,8 @@ template <typename number> class Vector;
 template <typename number> class BlockVector;
 template <class number> class FullMatrix;
 
+/*!@addtogroup mg */
+/*@{*/
 
 /**
  * This is a collection of functions operating on, and manipulating
@@ -220,5 +222,6 @@ class MGTools
                     std::vector<std::vector<unsigned int> >& cached_sizes);
 };
 
+/*@}*/
 
 #endif
index b721513f71a2afa63db60630605bef0c728564bf..95d65659b170cb3d9383540c01b73c8725a85951 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -19,6 +19,8 @@
 
 #include <boost/shared_ptr.hpp>
 
+/*!@addtogroup mg */
+/*@{*/
 
 /**
  * An array with an object for each level.  The purpose of this class
@@ -112,6 +114,7 @@ class MGLevelObject : public Subscriptor
     std::vector<boost::shared_ptr<Object> > objects;
 };
 
+/*@}*/
 
 /* ------------------------------------------------------------------- */
 
@@ -215,7 +218,4 @@ MGLevelObject<Object>::memory_consumption () const
 }
 
 
-/*-----------------------   mg_level_object.h     -----------------------*/
-
 #endif
-/*-----------------------   mg_level_object.h     -----------------------*/
index ca97f82f388bd98274926db2dfa6fc105b1875d4..60f77b44057c272b9d8e1c4fd24b0b4340301f0a 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 2003 by the deal.II authors
+//    Copyright (C) 2002, 2003, 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -17,6 +17,9 @@
 #include <multigrid/mg_base.h>
 #include <multigrid/mg_level_object.h>
 
+/*!@addtogroup mg */
+/*@{*/
+
 /**
  * Multilevel matrix. This class implements the interface defined by
  * MGMatrixBase, using MGLevelObject of an arbitrary
@@ -180,6 +183,7 @@ class MGMatrixSelect : public MGMatrixBase<Vector<number> >
     
 };
 
+/*@}*/
 
 /*----------------------------------------------------------------------*/
 
index fb65a10a1ae51409f90fd75b2144a545e73abd54..763f54be27222eae5d2e116d8e2334c229abb2fa 100644 (file)
@@ -29,6 +29,9 @@ template <int dim> class MGDoFHandler;
  * MGSmootherBase is defined in mg_base.h
  */
 
+/*!@addtogroup mg */
+/*@{*/
+
 /**
  * Smoother doing nothing. This class is not useful for many applications other
  * than for testing some multigrid procedures. Also some applications might
@@ -366,6 +369,7 @@ class MGSmootherRelaxation : public MGSmootherBase<VECTOR>
     VectorMemory<VECTOR>& mem;
 };
 
+/*@}*/
 
 /* ------------------------------- Inline functions -------------------------- */
 
index 1e34181d9422288cf0c413895cbf0ee8c8a0439b..652b4f6eaae8b64e851e092ba96cf97787a74e06 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors
+//    Copyright (C) 1998 - 2005 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
@@ -37,6 +37,9 @@ template <int dim> class MGDoFHandler;
  * MGTransferBase is defined in mg_base.h
  */
 
+/*!@addtogroup mg */
+/*@{*/
+
 /**
  * Implementation of the MGTransferBase interface for which the transfer
  * operations are prebuilt upon construction of the object of this class as
@@ -811,6 +814,8 @@ class MGTransferSelect : public MGTransferBase<Vector<number> >,
     unsigned int mg_selected_component;
 };
 
+/*@}*/
+
 //----------------------------------------------------------------------//
 template <typename number>
 inline void
index ad5746611c4c2aa792106390e3350df2e4c8eae7..1e92fbe06f46a832ea7807601ef7805371ac70b4 100644 (file)
 
 #include <vector>
 
-
-//TODO:[GK] Cleanup
-//  move definitions of virtual destructors to mg_base.templates.h
-//
+/*!@addtogroup mg */
+/*@{*/
 
 /**
  * Implementation of the multigrid method.
  *
- * Warning: multigrid on locally refined meshes only works with
+ * @warning multigrid on locally refined meshes only works with
  * <b>discontinuous finite elements</b> right now. It is not clear,
  * whether the paradigm of local smoothing we use is applicable to
  * continuous elements with hanging nodes; in fact, most people you
  * meet on conferences seem to deny this.
  *
- * The function actually performing a multi-level V-cycle,
- * @p level_v_step, 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:
- *
- * <tt>MATRIX</tt> is a matrix as specified for LAC-solvers, that is, it
- * has a function
- * @code
- *   void MATRIX::vmult(VECTOR& x, const VECTOR& b)
- * @endcode
- * performing the matrix vector product <i>x=Ab</i>.
+ * The function which starts a multigrid cycle on the finest level is
+ * cycle(). Depending on the cycle type chosen with the constructor
+ * (see enum Cycle), this function triggers one of the cycles
+ * level_v_step(), level_w_step() or level_f_step().
  *
- * <tt>SMOOTH</tt> is a class with a function
- * @code
- *   void SMOOTH::smooth (unsigned int level, VECTOR& x, const VECTOR& b)
- * @endcode
- * modifying the vector <i>x</i> such that the residual <i>b-Ax</i>
- * on level <i>l</i> is smoothened. Refer to MGSmootherRelaxation
- * for an example.
+ * @note The interface of this class is still very clumsy. In
+ * particular, you will have to set up quite a few auxiliary objects
+ * before you can use it. Unfortunately, it seems that this can be
+ * avoided only be restricting the flexibility of this class in an
+ * unacceptable way.
  *
- * <tt>COARSE</tt> has an
- * @code
- *   void COARSE::operator() (VECTOR& x, const VECTOR& b)
- * @endcode
- * returning the solution to the system <tt>Ax=b</tt> on the coarsest level
- * in <i>x</i>.
- * 
- * @author Guido Kanschat, 1999 - 2004
+ * @author Guido Kanschat, 1999 - 2005
  */
 template <class VECTOR>
 class Multigrid : public Subscriptor
 {
   public:
-  typedef VECTOR vector_type;
-  typedef const VECTOR const_vector_type;
+                                    /**
+                                     * List of implemented cycle types.
+                                     */
+    enum Cycle
+    {
+                                          /// The V-cycle
+         v_cycle,
+                                          /// The W-cycle (not working yet)
+         w_cycle,
+                                          /// The F-cycle (not working yet)
+         f_cycle
+    };
+    
+    typedef VECTOR vector_type;
+    typedef const VECTOR const_vector_type;
   
                                     /**
                                      * Constructor. The
-                                     * @p MGDoFHandler is used to
+                                     * MGDoFHandler is used to
                                      * determine the highest possible
-                                     * level. @p transfer is an
+                                     * level. <tt>transfer</tt> 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
-                                     * will be used later on in the
+                                     * will be used later in the
                                      * course of the
                                      * computations. You should
                                      * therefore create objects of
                                      * this type as late as possible.
                                      */
-  template <int dim>
+    template <int dim>
     Multigrid(const MGDoFHandler<dim>& mg_dof_handler,
              const MGMatrixBase<VECTOR>& matrix,
              const MGCoarseGridBase<VECTOR>& coarse,
              const MGTransferBase<VECTOR>& transfer,
              const MGSmootherBase<VECTOR>& pre_smooth,
              const MGSmootherBase<VECTOR>& post_smooth,
-             const unsigned int            minlevel = 0,
-             const unsigned int            maxlevel = 1000000);
-
+             Cycle cycle = v_cycle);
+                                    /**
+                                     * Experimental constructor for
+                                     * cases in which no MGDoFHandler
+                                     * is available.
+                                     *
+                                     * @warning Not intended for general use.
+                                     */
+    Multigrid(const unsigned int minlevel,
+             const unsigned int maxlevel,
+             const MGMatrixBase<VECTOR>& matrix,
+             const MGCoarseGridBase<VECTOR>& coarse,
+             const MGTransferBase<VECTOR>& transfer,
+             const MGSmootherBase<VECTOR>& pre_smooth,
+             const MGSmootherBase<VECTOR>& post_smooth,
+             Cycle cycle = v_cycle);
 
                                     /**
                                      * Reinit this class according to
-                                     * @p minlevel and @p maxlevel.
+                                     * #minlevel and #maxlevel.
                                      */
     void reinit (const unsigned int minlevel,
                 const unsigned int maxlevel);
 
+                                    /**
+                                     * Execute one multigrid
+                                     * cycle. The type of cycle is
+                                     * selected by the constructor
+                                     * argument cycle. See the enum
+                                     * Cycle for available types.
+                                     */
+    void cycle ();
+    
                                     /**
                                      * Execute one step of the
-                                     * v-cycle algorithm.  This
+                                     * 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
+                                     * multilevel vector #defect is
+                                     * filled with the residual of an
+                                     * outer defect correction
+                                     * scheme. This is usually taken
+                                     * care of by
                                      * PreconditionMG). After
-                                     * execution of <tt>vcycle()</tt>, the
-                                     * result is in the vector
-                                     * @p solution. See
+                                     * vcycle(), the result is in the
+                                     * multilevel vector
+                                     * #solution. See
                                      * <tt>copy_*_mg</tt> in class
-                                     * @p MGTools if you want to use
+                                     * MGTools if you want to use
                                      * these vectors yourself.
                                      *
                                      * The actual work for this
                                      * function is done in
-                                     * @p level_v_step.
+                                     * level_v_step().
                                      */
-    void vcycle();
+    void vcycle ();
 
                                     /**
                                      * Set additional matrices to
                                      * correct residual computation
-                                     * at refinement edges.
+                                     * at refinement edges. These
+                                     * matrices originate from
+                                     * discontinuous Galerkin methods
+                                     * (see FE_DGQ etc.), where they
+                                     * correspond tu the edge fluxes
+                                     * at the refinement edge between
+                                     * two levels.
                                      */
     void set_edge_matrices (const MGMatrixBase<VECTOR>& edge_down,
                            const MGMatrixBase<VECTOR>& edge_up);
 
+                                    /**
+                                     * Set the highest level for
+                                     * which the multilevel method is
+                                     * performed. By default, this is
+                                     * the finest level of the
+                                     * Triangulation; therefore, this
+                                     * function will only accept
+                                     * arguments smaller than the
+                                     * current #maxlevel.
+                                     */
+    void set_maxlevel(unsigned int);
+
+                                    /**
+                                     * Chance #cycle_type used in cycle().
+                                     */
+    void set_cycle(Cycle);
+    
+                                    /**
+                                     * Set the debug level. Higher
+                                     * values will create more
+                                     * debugging output during the
+                                     * multigrid cycles.
+                                     */
+    void set_debug(unsigned int);
+    
+                                    /**
+                                     * Exception.
+                                     */
+    DeclException2(ExcSwitchedLevels, int, int,
+                  << "minlevel and maxlevel switched, should be: "
+                  << arg1 << "<=" << arg2);
   private:
     
                                     /**
-                                     * The actual v-cycle multigrid method.
-                                     * @p level is the level the
-                                     * function should work on. It
+                                     * The V-cycle multigrid method.
+                                     * <tt>level</tt> is the level the
+                                     * function starts 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.
+                                     * recursively for <tt>level-1</tt>,
+                                     * unless we are on #minlevel
+                                     * where the coarse grid solver
+                                     * solves the problem exactly.
                                      */
     void level_v_step (const unsigned int level);
 
+                                    /**
+                                     * The actual W-cycle multigrid method.
+                                     * <tt>level</tt> is the level the
+                                     * function starts on. It
+                                     * will usually be called for the
+                                     * highest level from outside,
+                                     * but will then call itself
+                                     * recursively for <tt>level-1</tt>,
+                                     * unless we are on #minlevel
+                                     * where the coarse grid solver
+                                     * solves the problem exactly.
+                                     */
+    void level_w_step (const unsigned int level);
+
+                                    /**
+                                     * Cycle type performed by the method cycle().
+                                     */
+    Cycle cycle_type;
+    
                                     /**
                                      * Level for coarse grid solution.
                                      */
@@ -172,28 +238,29 @@ class Multigrid : public Subscriptor
                                      * Highest level of cells.
                                      */
     unsigned int maxlevel;
-    
+
+  public:
                                     /**
-                                     * Auxiliary vector. Contains the
-                                     * defect to be corrected before
-                                     * the multigrid step.
+                                     * Input vector for the
+                                     * cycle. Contains the defect of
+                                     * the outer method projected to
+                                     * the multilevel vectors.
                                      */
     MGLevelObject<VECTOR> defect;
 
                                     /**
-                                     * Auxiliary vector. Contains the
-                                     * updated solution after the
+                                     * The solution update after the
                                      * multigrid step.
                                      */
     MGLevelObject<VECTOR> solution;
     
+  private:
                                     /**
                                      * Auxiliary vector.
                                      */
     MGLevelObject<VECTOR> t;    
 
 
-  private:
                                     /**
                                      * The matrix for each level.
                                      */
@@ -230,11 +297,12 @@ class Multigrid : public Subscriptor
     SmartPointer<const MGMatrixBase<VECTOR> > edge_up;
 
                                     /**
-                                     * Exception.
+                                     * Level for debug
+                                     * output. Defaults to zero and
+                                     * can be set by set_debug().
                                      */
-    DeclException2(ExcSwitchedLevels, int, int,
-                  << "minlevel and maxlevel switched, should be: "
-                  << arg1 << "<=" << arg2);
+    unsigned int debug;
+    
     template<int dim, class VECTOR2, class TRANSFER> friend class PreconditionMG;
 };
 
@@ -244,7 +312,6 @@ class Multigrid : public Subscriptor
  * 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 <tt>void copy_to_mg(const VECTOR&)</tt>
  * to store @p src in the right hand side of the multi-level method and
  * <tt>void copy_from_mg(VECTOR&)</tt> to store the result of the v-cycle in @p dst.
@@ -324,13 +391,15 @@ class PreconditionMG : public Subscriptor
                                      */
     SmartPointer<Multigrid<VECTOR> > multigrid;
     
-                                  /**
-                                   * Object for grid tranfer.
-                                   */
+                                    /**
+                                     * Object for grid tranfer.
+                                     */
     SmartPointer<const TRANSFER> transfer;  
 };
 
+/*@}*/
 
+/// @if NoDoc
 /* --------------------------- inline functions --------------------- */
 
 
@@ -342,12 +411,11 @@ Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
                              const MGTransferBase<VECTOR>& transfer,
                              const MGSmootherBase<VECTOR>& pre_smooth,
                              const MGSmootherBase<VECTOR>& post_smooth,
-                             const unsigned int min_level,
-                             const unsigned int max_level)
+                             Multigrid<VECTOR>::Cycle cycle)
                :
-               minlevel(min_level),
-               maxlevel(std::min(mg_dof_handler.get_tria().n_levels()-1,
-                                 max_level)),
+               cycle_type(cycle),
+               minlevel(0),
+               maxlevel(mg_dof_handler.get_tria().n_levels()-1),
                defect(minlevel,maxlevel),
                solution(minlevel,maxlevel),
                t(minlevel,maxlevel),
@@ -357,23 +425,11 @@ Multigrid<VECTOR>::Multigrid (const MGDoFHandler<dim>& mg_dof_handler,
                pre_smooth(&pre_smooth),
                post_smooth(&post_smooth),
                edge_down(0),
-               edge_up(0)
+               edge_up(0),
+               debug(0)
 {}
 
 
-template <class VECTOR>
-void
-Multigrid<VECTOR>::reinit(const unsigned int min_level,
-                         const unsigned int max_level)
-{
-  minlevel=min_level;
-  maxlevel=max_level;
-  defect.resize(minlevel, maxlevel);
-  solution.resize(minlevel, maxlevel);
-  t.resize(minlevel, maxlevel);
-}
-
-  
 /* --------------------------- inline functions --------------------- */
 
 
@@ -405,7 +461,7 @@ PreconditionMG<dim, VECTOR, TRANSFER>::vmult (
   transfer->copy_to_mg(*mg_dof_handler,
                       multigrid->defect,
                       src);
-  multigrid->vcycle();
+  multigrid->cycle();
   transfer->copy_from_mg(*mg_dof_handler,
                         dst,
                         multigrid->solution);
@@ -422,7 +478,7 @@ PreconditionMG<dim, VECTOR, TRANSFER>::vmult_add (
   transfer->copy_to_mg(*mg_dof_handler,
                       multigrid->defect,
                       src);
-  multigrid->vcycle();
+  multigrid->cycle();
   transfer->copy_from_mg_add(*mg_dof_handler,
                             dst,
                             multigrid->solution);
@@ -450,5 +506,6 @@ PreconditionMG<dim, VECTOR, TRANSFER>::Tvmult_add (
   Assert(false, ExcNotImplemented());
 }
 
+/// @endif
 
 #endif
index 810421da40e7602ff23688df2d493c38a1660685..c6cc4543916ef8a22341d017aa5104b1759bd182 100644 (file)
 
 using namespace std;
 
+
+template <class VECTOR>
+Multigrid<VECTOR>::Multigrid (const unsigned int minlevel,
+                             const unsigned int maxlevel,
+                             const MGMatrixBase<VECTOR>& matrix,
+                             const MGCoarseGridBase<VECTOR>& coarse,
+                             const MGTransferBase<VECTOR>& transfer,
+                             const MGSmootherBase<VECTOR>& pre_smooth,
+                             const MGSmootherBase<VECTOR>& post_smooth,
+                             Multigrid<VECTOR>::Cycle cycle)
+               :
+               cycle_type(cycle),
+               minlevel(minlevel),
+               maxlevel(maxlevel),
+               defect(minlevel,maxlevel),
+               solution(minlevel,maxlevel),
+               t(minlevel,maxlevel),
+               matrix(&matrix),
+               coarse(&coarse),
+               transfer(&transfer),
+               pre_smooth(&pre_smooth),
+               post_smooth(&post_smooth),
+               edge_down(0),
+               edge_up(0),
+               debug(0)
+{}
+
+
+
+//TODO: This function cannot work. At least maxlevel should be tested!
+template <class VECTOR>
+void
+Multigrid<VECTOR>::reinit(const unsigned int min_level,
+                         const unsigned int max_level)
+{
+  Assert(false, ExcNotImplemented());
+  
+  minlevel=min_level;
+  maxlevel=max_level;
+  defect.resize(minlevel, maxlevel);
+  solution.resize(minlevel, maxlevel);
+  t.resize(minlevel, maxlevel);
+}
+
+
+template <class VECTOR>
+void
+Multigrid<VECTOR>::set_maxlevel(unsigned int l)
+{
+  Assert (l <= maxlevel, ExcIndexRange(l,minlevel,maxlevel+1));
+  Assert (l >= minlevel, ExcIndexRange(l,minlevel,maxlevel+1));
+  maxlevel = l;
+}
+
+
+template <class VECTOR>
+void
+Multigrid<VECTOR>::set_cycle(Multigrid<VECTOR>::Cycle c)
+{
+  cycle_type = c;
+}
+
+
+template <class VECTOR>
+void
+Multigrid<VECTOR>::set_debug(unsigned int d)
+{
+  debug = d;
+}
+
+
 template <class VECTOR>
 void
 Multigrid<VECTOR>::set_edge_matrices (const MGMatrixBase<VECTOR>& down,
@@ -36,25 +107,99 @@ template <class VECTOR>
 void
 Multigrid<VECTOR>::level_v_step(const unsigned int level)
 {
-  solution[level] = 0.;
+  if (debug>0)
+    deallog << "V-cycle entering level " << level << std::endl;
   
   if (level == minlevel)
     {
       (*coarse)(level, solution[level], defect[level]);
       return;
     }
+  if (debug>1)
+    deallog << "Smoothing on     level " << level << std::endl;
+                                  // smoothing of the residual by
+                                  // modifying s
+  pre_smooth->smooth(level, solution[level], defect[level]);
+
+  if (debug>1)
+    deallog << "Residual on      level " << level << std::endl;
+                                  // t = A*solution[level]
+  matrix->vmult(level, t[level], solution[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.;
+      if (l==level && edge_down != 0)
+       {
+         edge_down->vmult(level, t[level-1], solution[level]);
+       }
+      transfer->restrict_and_add (l, t[l-1], t[l]);
+      defect[l-1] -= t[l-1];
+    }
+
+                                  // do recursion
+  solution[level-1] = 0.;
+  level_v_step(level-1);
+
+                                  // reset size of the auxiliary
+                                  // vector, since it has been
+                                  // resized in the recursive call to
+                                  // level_v_step directly above
+  t[level] = 0.;
+
+                                  // do coarse grid correction
+  transfer->prolongate(level, t[level], solution[level-1]);
+  
+  solution[level] += t[level];
+
+  if (edge_up != 0)
+    {
+      edge_up->Tvmult(level, t[level], solution[level-1]);
+      defect[level] -= t[level];
+    }
+  
+  if (debug>1)
+    deallog << "Smoothing on     level " << level << std::endl;
+                                  // post-smoothing
+  post_smooth->smooth(level, solution[level], defect[level]);
+
+  if (debug>1)
+    deallog << "V-cycle leaving  level " << level << std::endl;
+}
 
-//  deallog << "Pre-smooth " << level << endl;
+
+
+template <class VECTOR>
+void
+Multigrid<VECTOR>::level_w_step(const unsigned int level)
+{
+  if (level == minlevel)
+    {
+      if (debug>0)
+       deallog << "W-cycle solving  level " << level << std::endl;
   
+      (*coarse)(level, solution[level], defect[level]);
+      return;
+    }
+  if (debug>0)
+    deallog << "W-cycle entering level " << level << std::endl;
+  
+  if (debug>1)
+    deallog << "Smoothing on     level " << level << std::endl;
                                   // smoothing of the residual by
                                   // modifying s
   pre_smooth->smooth(level, solution[level], defect[level]);
 
-//  deallog << "vmult " << level << endl;
+  if (debug>1)
+    deallog << "Residual on      level " << level << std::endl;
                                   // t = A*solution[level]
   matrix->vmult(level, t[level], solution[level]);
 
-//  deallog << "restrict " << level << endl;
                                   // make t rhs of lower level The
                                   // non-refined parts of the
                                   // coarse-level defect already
@@ -71,9 +216,10 @@ Multigrid<VECTOR>::level_v_step(const unsigned int level)
       defect[l-1] -= t[l-1];
     }
 
-//    deallog << "recursion " << level << endl;
                                   // do recursion
-  level_v_step(level-1);
+  solution[level-1] = 0.;
+  level_w_step(level-1);
+  level_w_step(level-1);
 
                                   // reset size of the auxiliary
                                   // vector, since it has been
@@ -82,7 +228,6 @@ Multigrid<VECTOR>::level_v_step(const unsigned int level)
   t[level] = 0.;
 
                                   // do coarse grid correction
-//  deallog << "prolongate " << level << endl;
   transfer->prolongate(level, t[level], solution[level-1]);
   
   solution[level] += t[level];
@@ -92,16 +237,20 @@ Multigrid<VECTOR>::level_v_step(const unsigned int level)
       edge_up->Tvmult(level, t[level], solution[level-1]);
       defect[level] -= t[level];
     }
-//  deallog << "Post-smooth " << level << endl;
+  
+  if (debug>1)
+    deallog << "Smoothing on     level " << level << std::endl;
                                   // post-smoothing
   post_smooth->smooth(level, solution[level], defect[level]);
-//  deallog << "ready " << level << endl;
+
+  if (debug>1)
+    deallog << "W-cycle leaving  level " << level << std::endl;
 }
 
 
 template <class VECTOR>
 void
-Multigrid<VECTOR>::vcycle()
+Multigrid<VECTOR>::cycle()
 {
                                   // The defect vector has been
                                   // initialized by copy_to_mg. Now
@@ -114,9 +263,37 @@ Multigrid<VECTOR>::vcycle()
       solution[level].reinit(defect[level]);
       t[level].reinit(defect[level]);
     }
+  
+  switch(cycle_type)
+    {
+      case v_cycle:
+       level_v_step (maxlevel);
+       break;
+      case w_cycle:
+       level_w_step (maxlevel);
+       break;
+      case f_cycle:
+       break;
+    }
+}
+
 
+template <class VECTOR>
+void
+Multigrid<VECTOR>::vcycle()
+{
+                                  // 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_v_step (maxlevel);
-//  abort ();
 }
 
 

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.