From: guido Date: Mon, 6 Jan 2003 19:19:07 +0000 (+0000) Subject: new MGSmootherRelaxation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fe2b942b3659b9c1387f8e16ae795014255ab1d3;p=dealii-svn.git new MGSmootherRelaxation git-svn-id: https://svn.dealii.org/trunk@6874 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index cb1b87311f..131b782eaa 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -15,8 +15,9 @@ #include -#include #include +#include +#include #include template class MGDoFHandler; @@ -44,9 +45,9 @@ class MGSmootherIdentity : public MGSmoother * operator equals the null * operator. */ - virtual void smooth (const unsigned int level, - VECTOR& u, - const VECTOR& rhs) const; + virtual void smooth (const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const; }; @@ -160,89 +161,193 @@ class MGSmootherContinuous : public MGSmoother /** - * A smoother using matrix builtin relaxation methods. + * Smoother using relaxation classes. + * + * This class performs smoothing on each level. The opreation can be + * controlled by several parameters. First, the relaxation parameter + * @p{omega} is used in the underlying relaxation method. @p{steps} is + * the number of relaxation steps on the finest level (on all levels + * if @p{variable} is off). If @p{variable} is @p{true}, the number of + * smoothing steps is doubled on each corser level. This results in a + * method having the complexity of the W-cycle, but saving grid + * transfers. This is the method proposed by Bramble at al. + * + * The option @p{symmetric} switches on alternating between the + * smoother and its transpose in each step as proposed by Bramble. * - * 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). + * @p{transpose} and @p{reverse} are similar in the effect that + * instead of the smoother the transposed is used. Typically, this is + * off for pre-smoothing and on for post-smoothing. While + * @p{transpose} is the true transposed smoothing operation, + * @p{reverse} just uses the transposed of the smoother, but the + * non-transposed matrix-vector multiplication; this can be used to + * invert the direction of the Gauss-Seidel method. + * + * If you are using block matrices, the second @p{initialize} function + * offers the possibility to extract a single block for smoothing. In + * this case, the multigrid method must be used only with the vector + * associated to that sinlge block. * * 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 + * @author Guido Kanschat, 2003 */ -template -class MGSmootherRelaxation : public MGSmootherContinuous +template +class MGSmootherRelaxation : public MGSmoother { public: /** - * Type of the smoothing - * function of the matrix. + * Constructor. Sets memory and smoothing parameters. */ - typedef void (MATRIX::* function_ptr)(VECTOR&, - const VECTOR&, - typename MATRIX::value_type) const; - + MGSmootherRelaxation(VectorMemory& mem, + const unsigned int steps = 1, + const bool variable = false, + const bool symmetric = false, + const bool transpose = false, + const bool reverse = false); + /** - * Constructor. - * The constructor uses an @p{MGDoFHandler} to initialize - * data structures for interior level boundary handling - * in @p{MGSmootherContinuous}. + * Initialize for matrices. The + * parameter @p{matrices} can be + * any object having functions + * @p{get_minlevel()} and + * @p{get_maxlevel()} as well as + * an @p{operator[]} returning a + * reference to @p{MATRIX}. This + * function stores pointers to + * the level matrices and + * initializes the smoothing + * operator for each level. * - * Furthermore, it takes a pointer to the - * level matrices and their smoothing function. - * This function must perform one relaxation step - * like @p{SparseMatrix::SOR_step} does. Do not - * use the preconditioning methods because they apply - * a preconditioning matrix to the residual. + * @p{additional_data} is an + * object of type + * @p{RELAX::AdditionalData} and + * is handed to the + * initialization function of the + * relaxation method. + */ + template + void initialize (const MGMATRIXOBJECT& matrices, + const DATA& additional_data); + + /** + * Initialize for single blocks + * of matrices. The parameter + * @p{matrices} can be any object + * having functions + * @p{get_minlevel()} and + * @p{get_maxlevel()} as well as + * an @p{operator[]} returning a + * reference to a block matrix + * where each block is of type + * @p{MATRIX}. Of this block + * matrix, the block indicated by + * @p{block_row} and + * @p{block_col} is selected on + * each level. This function + * stores pointers to the level + * matrices and initializes the + * smoothing operator for each + * level. * - * The final two parameters are the number of relaxation - * steps and the relaxation parameter. + * @p{additional_data} is an + * object of type + * @p{RELAX::AdditionalData} and + * is handed to the + * initialization function of the + * relaxation method. */ + template + void initialize (const MGMATRIXOBJECT& matrices, + const DATA& additional_data, + const unsigned int block_row, + const unsigned int block_col); - template - MGSmootherRelaxation(const MGDoFHandler& mg_dof, - const MGLevelObject& matrix, - const function_ptr function, - const unsigned int steps, - const double omega = 1.); + /** + * Empty all vectors. + */ + void clear (); /** - * We use the relaxation method in - * @p{MATRIX} for the real - * work and find a way to keep - * the boundary values. + * Modify the number of smoothing + * steps on finest level. + */ + void set_steps (const unsigned int); + + /** + * Switch on/off variable smoothing. + */ + void set_variable (const bool); + + /** + * Switch on/off symmetric smoothing. + */ + void set_symmetric (const bool); + + /** + * Switch on/off transposed. This + * is mutually exclusive with + * @ref{reverse}. + */ + void set_transpose (const bool); + + /** + * Switch on/off reversed. This + * is mutually exclusive with + * @ref{transpose}. + */ + void set_reverse (const bool); + + /** + * The actual smoothing method. */ virtual void smooth (const unsigned int level, VECTOR& u, const VECTOR& rhs) const; - private: + + /** + * Object containing relaxation methods. + */ + MGLevelObject smoothers; + + private: /** * Pointer to the matrices. */ - SmartPointer > matrix; - + MGLevelObject > matrices; + /** - * Pointer to the relaxation function. + * Number of smoothing steps. */ - function_ptr relaxation; + unsigned int steps; /** - * Relaxation parameter. + * Variable smoothing? */ - double omega; - + bool variable; + /** - * Auxiliary vector. + * Symmetric smoothing? */ - mutable VECTOR h1; - + bool symmetric; + + /* + * Transposed? + */ + bool transpose; + + /** + * Reverse? + */ + bool reverse; + /** - * Auxiliary vector. + * Memory for auxiliary vectors. */ - mutable VECTOR h2; + VectorMemory& mem; }; @@ -273,5 +378,173 @@ MGSmootherContinuous::get_steps() const return steps; } +//----------------------------------------------------------------------// + +template +inline +MGSmootherRelaxation::MGSmootherRelaxation( + VectorMemory& mem, + const unsigned int steps, + const bool variable, + const bool symmetric, + const bool transpose, + const bool reverse) + : + steps(steps), + variable(variable), + symmetric(symmetric), + transpose(transpose), + reverse(reverse), + mem(mem) +{} + + +template +inline void +MGSmootherRelaxation::clear () +{ + smoothers.clear(); + matrices.clear(); +} + + +template +template +inline void +MGSmootherRelaxation::initialize ( + const MGMATRIXOBJECT& m, + const DATA& data) +{ + const unsigned int min = m.get_minlevel(); + const unsigned int max = m.get_maxlevel(); + + matrices.resize(min, max); + smoothers.resize(min, max); + + for (unsigned int i=min;i<=max;++i) + { + matrices[i] = &m[i]; + smoothers[i].initialize(m[i], data); + } +} + + +template +template +inline void +MGSmootherRelaxation::initialize ( + const MGMATRIXOBJECT& m, + const DATA& data, + const unsigned int row, + const unsigned int col) +{ + const unsigned int min = m.get_minlevel(); + const unsigned int max = m.get_maxlevel(); + + matrices.resize(min, max); + smoothers.resize(min, max); + + for (unsigned int i=min;i<=max;++i) + { + matrices[i] = &(m[i].block(row, col)); + smoothers[i].initialize(m[i].block(row, col), data); + } +} + + +template +inline void +MGSmootherRelaxation:: +set_steps (const unsigned int s) +{ + steps = s; +} + + +template +inline void +MGSmootherRelaxation:: +set_variable (const bool flag) +{ + variable = flag; +} + + +template +inline void +MGSmootherRelaxation:: +set_symmetric (const bool flag) +{ + symmetric = flag; +} + + +template +inline void +MGSmootherRelaxation:: +set_transpose (const bool flag) +{ + transpose = flag; +} + + +template +inline void +MGSmootherRelaxation:: +set_reverse (const bool flag) +{ + reverse = flag; +} + + +template +inline void +MGSmootherRelaxation::smooth( + const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const +{ + unsigned int maxlevel = matrices.get_maxlevel(); + unsigned int steps2 = steps; + + if (variable) + steps2 *= (1<<(maxlevel-level)); + + Vector* r = mem.alloc(); + Vector* d = mem.alloc(); + r->reinit(u); + d->reinit(u); + + bool T = transpose; + if (symmetric && (steps2 % 2 == 0)) + T = false; +// cerr << 'S' << level; + + for (unsigned int i=0; isadd(-1.,1.,rhs); + smoothers[level].Tvmult(*d, *r); + } else { +// cerr << 'N'; + matrices[level].vmult(*r,u); + r->sadd(-1.,1.,rhs); + smoothers[level].vmult(*d, *r); + } +// cerr << '{' << r->l2_norm() << '}'; + u += *d; + if (symmetric) + T = !T; + } + + mem.free(r); + mem.free(d); +} #endif diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h index fb56334d44..ec7a3cada5 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.templates.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.templates.h @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -17,40 +17,4 @@ #include -template -template -MGSmootherRelaxation -::MGSmootherRelaxation (const MGDoFHandler& mg_dof, - const MGLevelObject& matrix, - const function_ptr relaxation, - const unsigned int steps, - const double omega) - : - MGSmootherContinuous(mg_dof, steps), - matrix(&matrix), - relaxation(relaxation), - omega(omega) -{} - - -template -void -MGSmootherRelaxation -::smooth (const unsigned int level, - VECTOR& u, - const VECTOR& rhs) const -{ - h1.reinit(u); - h2.reinit(u); - for(unsigned i=0; iget_steps(); ++i) - { - (*matrix)[level].residual(h1, u, rhs); - this->set_zero_interior_boundary(level,h1); - ((*matrix)[level].*relaxation)(h2,h1,omega); - this->set_zero_interior_boundary(level,h2); - u.add(h2); - } -} - - #endif diff --git a/deal.II/deal.II/source/multigrid/mg_smoother.cc b/deal.II/deal.II/source/multigrid/mg_smoother.cc index 9ef239b4d1..45bd4a7d9a 100644 --- a/deal.II/deal.II/source/multigrid/mg_smoother.cc +++ b/deal.II/deal.II/source/multigrid/mg_smoother.cc @@ -2,7 +2,7 @@ // $Id$ // Version: $Name$ // -// Copyright (C) 1998, 1999, 2000, 2001, 2002 by the deal.II authors +// Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003 by the deal.II authors // // This file is subject to QPL and may not be distributed // without copyright and license information. Please refer @@ -177,70 +177,3 @@ template void MGSmootherContinuous >::set_zero_interior_boundary ( const unsigned int, BlockVector&) const; - - - -template -MGSmootherRelaxation, Vector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, Vector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, Vector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, Vector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - - -template -MGSmootherRelaxation, BlockVector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, BlockVector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, BlockVector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double); - -template -MGSmootherRelaxation, BlockVector > -::MGSmootherRelaxation(const MGDoFHandler&, - const MGLevelObject >&, - const function_ptr, - const unsigned int, - const double);