From 62ade8b06fb6263eab3343ad3441da4eefe9dba3 Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Tue, 8 Mar 2005 15:42:11 +0000 Subject: [PATCH] Multigrid prepared for W- and F-cycles 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 | 9 +- deal.II/deal.II/include/multigrid/mg_coarse.h | 8 +- .../include/multigrid/mg_dof_accessor.h | 6 +- .../include/multigrid/mg_dof_handler.h | 11 +- .../deal.II/include/multigrid/mg_dof_tools.h | 5 +- .../include/multigrid/mg_level_object.h | 8 +- deal.II/deal.II/include/multigrid/mg_matrix.h | 6 +- .../deal.II/include/multigrid/mg_smoother.h | 4 + .../deal.II/include/multigrid/mg_transfer.h | 7 +- deal.II/deal.II/include/multigrid/multigrid.h | 265 +++++++++++------- .../include/multigrid/multigrid.templates.h | 199 ++++++++++++- 11 files changed, 392 insertions(+), 136 deletions(-) diff --git a/deal.II/deal.II/include/multigrid/mg_base.h b/deal.II/deal.II/include/multigrid/mg_base.h index 80bcedc2bd..a771e59089 100644 --- a/deal.II/deal.II/include/multigrid/mg_base.h +++ b/deal.II/deal.II/include/multigrid/mg_base.h @@ -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 - +/*!@addtogroup mg */ +/*@{*/ /** @@ -226,9 +227,7 @@ class MGSmootherBase : public Subscriptor const VECTOR& rhs) const = 0; }; +/*@}*/ -/*---------------------------- mgbase.h ---------------------------*/ - #endif -/*---------------------------- mgbase.h ---------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_coarse.h b/deal.II/deal.II/include/multigrid/mg_coarse.h index 10a8e4d4d4..8804b74b3a 100644 --- a/deal.II/deal.II/include/multigrid/mg_coarse.h +++ b/deal.II/deal.II/include/multigrid/mg_coarse.h @@ -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 #include +/*!@addtogroup mg */ +/*@{*/ /** * Coarse grid solver using LAC iterative methods. @@ -158,8 +160,9 @@ class MGCoarseGridHouseholder : public MGCoarseGridBase mutable VECTOR aux; }; +/*@}*/ - +/// @if NoDoc /* ------------------ Functions for MGCoarseGridLACIteration ------------ */ @@ -273,5 +276,6 @@ MGCoarseGridHouseholder::operator() ( work.least_squares(dst, aux); } +/// @endif #endif diff --git a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h index 912df74557..8ce17013e7 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_accessor.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_accessor.h @@ -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 class MGDoFObjectAccessor<1, dim>; template class MGDoFObjectAccessor<2, dim>; template 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 { DeclException0 (ExcNotUsefulForThisDimension); }; +/*@}*/ + template<> MGDoFCellAccessor<1>::face_iterator MGDoFCellAccessor<1>::face (const unsigned int i) const; diff --git a/deal.II/deal.II/include/multigrid/mg_dof_handler.h b/deal.II/deal.II/include/multigrid/mg_dof_handler.h index 0904913095..2ff32ac04f 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_handler.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_handler.h @@ -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 #include @@ -27,9 +24,10 @@ template class MGDoFObjectAccessor<1, dim>; template class MGDoFObjectAccessor<2, dim>; template 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 template friend class MGDoFObjectAccessor; }; +/*@}*/ /* ----------------------- Inline functions of MGDoFHandler -------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_dof_tools.h b/deal.II/deal.II/include/multigrid/mg_dof_tools.h index 7d0a3f894e..746b3ac80d 100644 --- a/deal.II/deal.II/include/multigrid/mg_dof_tools.h +++ b/deal.II/deal.II/include/multigrid/mg_dof_tools.h @@ -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 class Vector; template class BlockVector; template class FullMatrix; +/*!@addtogroup mg */ +/*@{*/ /** * This is a collection of functions operating on, and manipulating @@ -220,5 +222,6 @@ class MGTools std::vector >& cached_sizes); }; +/*@}*/ #endif diff --git a/deal.II/deal.II/include/multigrid/mg_level_object.h b/deal.II/deal.II/include/multigrid/mg_level_object.h index b721513f71..95d65659b1 100644 --- a/deal.II/deal.II/include/multigrid/mg_level_object.h +++ b/deal.II/deal.II/include/multigrid/mg_level_object.h @@ -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 +/*!@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 > objects; }; +/*@}*/ /* ------------------------------------------------------------------- */ @@ -215,7 +218,4 @@ MGLevelObject::memory_consumption () const } -/*----------------------- mg_level_object.h -----------------------*/ - #endif -/*----------------------- mg_level_object.h -----------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_matrix.h b/deal.II/deal.II/include/multigrid/mg_matrix.h index ca97f82f38..60f77b4405 100644 --- a/deal.II/deal.II/include/multigrid/mg_matrix.h +++ b/deal.II/deal.II/include/multigrid/mg_matrix.h @@ -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 #include +/*!@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 > }; +/*@}*/ /*----------------------------------------------------------------------*/ diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index fb65a10a1a..763f54be27 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -29,6 +29,9 @@ template 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 VectorMemory& mem; }; +/*@}*/ /* ------------------------------- Inline functions -------------------------- */ diff --git a/deal.II/deal.II/include/multigrid/mg_transfer.h b/deal.II/deal.II/include/multigrid/mg_transfer.h index 1e34181d94..652b4f6eaa 100644 --- a/deal.II/deal.II/include/multigrid/mg_transfer.h +++ b/deal.II/deal.II/include/multigrid/mg_transfer.h @@ -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 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 >, unsigned int mg_selected_component; }; +/*@}*/ + //----------------------------------------------------------------------// template inline void diff --git a/deal.II/deal.II/include/multigrid/multigrid.h b/deal.II/deal.II/include/multigrid/multigrid.h index ad5746611c..1e92fbe06f 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.h +++ b/deal.II/deal.II/include/multigrid/multigrid.h @@ -25,144 +25,210 @@ #include - -//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 * discontinuous finite elements 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: - * - * MATRIX 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 x=Ab. + * 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(). * - * SMOOTH is a class with a function - * @code - * void SMOOTH::smooth (unsigned int level, VECTOR& x, const VECTOR& b) - * @endcode - * modifying the vector x such that the residual b-Ax - * on level l 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. * - * COARSE has an - * @code - * void COARSE::operator() (VECTOR& x, const VECTOR& b) - * @endcode - * returning the solution to the system Ax=b on the coarsest level - * in x. - * - * @author Guido Kanschat, 1999 - 2004 + * @author Guido Kanschat, 1999 - 2005 */ template 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. 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 - * 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 + template Multigrid(const MGDoFHandler& mg_dof_handler, const MGMatrixBase& matrix, const MGCoarseGridBase& coarse, const MGTransferBase& transfer, const MGSmootherBase& pre_smooth, const MGSmootherBase& 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& matrix, + const MGCoarseGridBase& coarse, + const MGTransferBase& transfer, + const MGSmootherBase& pre_smooth, + const MGSmootherBase& 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 vcycle(), the - * result is in the vector - * @p solution. See + * vcycle(), the result is in the + * multilevel vector + * #solution. See * copy_*_mg 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& edge_down, const MGMatrixBase& 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. + * level 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 level-1, + * 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. + * level 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 level-1, + * 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 defect; /** - * Auxiliary vector. Contains the - * updated solution after the + * The solution update after the * multigrid step. */ MGLevelObject solution; + private: /** * Auxiliary vector. */ MGLevelObject t; - private: /** * The matrix for each level. */ @@ -230,11 +297,12 @@ class Multigrid : public Subscriptor SmartPointer > 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 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 void copy_to_mg(const VECTOR&) * to store @p src in the right hand side of the multi-level method and * void copy_from_mg(VECTOR&) to store the result of the v-cycle in @p dst. @@ -324,13 +391,15 @@ class PreconditionMG : public Subscriptor */ SmartPointer > multigrid; - /** - * Object for grid tranfer. - */ + /** + * Object for grid tranfer. + */ SmartPointer transfer; }; +/*@}*/ +/// @if NoDoc /* --------------------------- inline functions --------------------- */ @@ -342,12 +411,11 @@ Multigrid::Multigrid (const MGDoFHandler& mg_dof_handler, const MGTransferBase& transfer, const MGSmootherBase& pre_smooth, const MGSmootherBase& post_smooth, - const unsigned int min_level, - const unsigned int max_level) + Multigrid::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::Multigrid (const MGDoFHandler& mg_dof_handler, pre_smooth(&pre_smooth), post_smooth(&post_smooth), edge_down(0), - edge_up(0) + edge_up(0), + debug(0) {} -template -void -Multigrid::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::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::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::Tvmult_add ( Assert(false, ExcNotImplemented()); } +/// @endif #endif diff --git a/deal.II/deal.II/include/multigrid/multigrid.templates.h b/deal.II/deal.II/include/multigrid/multigrid.templates.h index 810421da40..c6cc454391 100644 --- a/deal.II/deal.II/include/multigrid/multigrid.templates.h +++ b/deal.II/deal.II/include/multigrid/multigrid.templates.h @@ -20,6 +20,77 @@ using namespace std; + +template +Multigrid::Multigrid (const unsigned int minlevel, + const unsigned int maxlevel, + const MGMatrixBase& matrix, + const MGCoarseGridBase& coarse, + const MGTransferBase& transfer, + const MGSmootherBase& pre_smooth, + const MGSmootherBase& post_smooth, + Multigrid::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 +void +Multigrid::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 +void +Multigrid::set_maxlevel(unsigned int l) +{ + Assert (l <= maxlevel, ExcIndexRange(l,minlevel,maxlevel+1)); + Assert (l >= minlevel, ExcIndexRange(l,minlevel,maxlevel+1)); + maxlevel = l; +} + + +template +void +Multigrid::set_cycle(Multigrid::Cycle c) +{ + cycle_type = c; +} + + +template +void +Multigrid::set_debug(unsigned int d) +{ + debug = d; +} + + template void Multigrid::set_edge_matrices (const MGMatrixBase& down, @@ -36,25 +107,99 @@ template void Multigrid::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 +void +Multigrid::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::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::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::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 void -Multigrid::vcycle() +Multigrid::cycle() { // The defect vector has been // initialized by copy_to_mg. Now @@ -114,9 +263,37 @@ Multigrid::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 +void +Multigrid::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 (); } -- 2.39.5