From: kanschat Date: Wed, 29 Sep 2010 00:41:59 +0000 (+0000) Subject: add MG namespace for reimplementation X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6094eb64ba7666ee8ddb2676705eadf309d165d7;p=dealii-svn.git add MG namespace for reimplementation git-svn-id: https://svn.dealii.org/trunk@22186 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/multigrid/mg_matrix.h b/deal.II/deal.II/include/multigrid/mg_matrix.h index 274b1eb3ea..4dc6584170 100644 --- a/deal.II/deal.II/include/multigrid/mg_matrix.h +++ b/deal.II/deal.II/include/multigrid/mg_matrix.h @@ -14,16 +14,75 @@ #define __deal2__mg_matrix_h #include +#include #include #include #include +#include DEAL_II_NAMESPACE_OPEN /*!@addtogroup mg */ /*@{*/ +namespace MG +{ + +/** + * Multilevel matrix. This matrix stores an MGLevelObject of + * PointerMatrixBase objects. It implementes the interface defined in + * MGMatrixBase, so that it can be used as a matrix in Multigrid. + * + * @author Guido Kanschat + * @date 2002, 2010 + */ + template > + class Matrix + : public MGMatrixBase + { + public: + /** + * Default constructor for an + * empty object. + */ + Matrix(); + + /** + * Constructor setting up + * pointers to the matrices in + * M by calling initialize(). + */ + template + Matrix(const MGLevelObject& M); + + /** + * Initialize the object such + * that the level + * multiplication uses the + * matrices in M + */ + template + void + initialize(const MGLevelObject& M); + + virtual void vmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const; + virtual void vmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const; + virtual void Tvmult (const unsigned int level, VECTOR& dst, const VECTOR& src) const; + virtual void Tvmult_add (const unsigned int level, VECTOR& dst, const VECTOR& src) const; + + /** + * Memory used by this object. + */ + unsigned int memory_consumption () const; + private: + MGLevelObject > > matrices; + }; + +} + /** + * @deprecated Use the much simpler class MG::Matrix instead. + * * Multilevel matrix. This class implements the interface defined by * MGMatrixBase, using MGLevelObject of an arbitrary * matrix class. @@ -198,6 +257,84 @@ class MGMatrixSelect : public MGMatrixBase > /*----------------------------------------------------------------------*/ +namespace MG +{ + template + template + inline void + Matrix::initialize (const MGLevelObject& p) + { + matrices.resize(p.min_level(), p.max_level()); + for (unsigned int level=p.min_level();level <= p.max_level();++level) + matrices[level] = std_cxx1x::shared_ptr > (new_pointer_matrix_base(p[level], VECTOR())); + } + + + template + template + Matrix::Matrix (const MGLevelObject& p) + { + initialize(p); + } + + + template + Matrix::Matrix () + {} + + + template + void + Matrix::vmult ( + const unsigned int level, + VECTOR& dst, + const VECTOR& src) const + { + matrices[level]->vmult(dst, src); + } + + + template + void + Matrix::vmult_add ( + const unsigned int level, + VECTOR& dst, + const VECTOR& src) const + { + matrices[level]->vmult_add(dst, src); + } + + + template + void + Matrix::Tvmult (const unsigned int level, + VECTOR& dst, + const VECTOR& src) const + { + matrices[level]->Tvmult(dst, src); + } + + + template + void + Matrix::Tvmult_add (const unsigned int level, + VECTOR& dst, + const VECTOR& src) const + { + matrices[level]->Tvmult_add(dst, src); + } + + + template + unsigned int + Matrix::memory_consumption () const + { + return sizeof(*this) + matrices->memory_consumption(); + } +} + +/*----------------------------------------------------------------------*/ + template MGMatrix::MGMatrix (MGLevelObject* p) : diff --git a/deal.II/deal.II/include/multigrid/mg_smoother.h b/deal.II/deal.II/include/multigrid/mg_smoother.h index 39e496d89f..1241e64697 100644 --- a/deal.II/deal.II/include/multigrid/mg_smoother.h +++ b/deal.II/deal.II/include/multigrid/mg_smoother.h @@ -193,6 +193,118 @@ class MGSmootherIdentity : public MGSmootherBase }; +namespace MG +{ +/** + * Smoother using relaxation classes. + * + * A relaxation class is an object that has two member functions, + * @code + * void step(VECTOR& x, const VECTOR& d) const; + * void Tstep(VECTOR& x, const VECTOR& d) const; + * @endcode + * performing one step of the smoothing scheme. + * + * This class performs smoothing on each level. The operation 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 coarser 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. + * + * @p transpose uses the transposed smoothing operation using + * Tstep instead of the regular step of the + * relaxation scheme. + * + * 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 single block. + * + * @author Guido Kanschat, + * @date 2003, 2009, 2010 + */ +template +class SmootherRelaxation : public MGSmoother +{ + public: + /** + * Constructor. Sets memory and + * smoothing parameters. + */ + SmootherRelaxation(const unsigned int steps = 1, + const bool variable = false, + const bool symmetric = false, + const bool transpose = false); + + /** + * Initialize for matrices. This + * function initializes the + * smoothing operator with the + * same smoother for each level. + * + * @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 MGLevelObject& matrices, + const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData()); + + template + void initialize (const MGLevelObject& matrices, + const MGLevelObject& additional_data); + /** + * Initialize for matrix + * blocks. This function + * initializes the smoothing + * operator with the same + * smoother for each level. + * + * @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 MGLevelObject >& matrices, +// const typename RELAX::AdditionalData & additional_data = typename RELAX::AdditionalData()); + + /** + * Empty all vectors. + */ + void clear (); + + /** + * The actual smoothing method. + */ + virtual void smooth (const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const; + + /** + * Memory used by this object. + */ + unsigned int memory_consumption () const; + private: + /** + * Object containing relaxation + * methods. + */ + MGLevelObject smoothers; +}; +} + /** * Smoother using relaxation classes. * @@ -633,6 +745,110 @@ MGSmoother::set_transpose (const bool flag) transpose = flag; } +//----------------------------------------------------------------------// + +namespace MG +{ + template + inline + SmootherRelaxation::SmootherRelaxation( + const unsigned int steps, + const bool variable, + const bool symmetric, + const bool transpose) + : MGSmoother(steps, variable, symmetric, transpose) + {} + + + template + inline void + SmootherRelaxation::clear () + { + smoothers.clear(); + } + + + template + template + inline void + SmootherRelaxation::initialize ( + const MGLevelObject& m, + const typename RELAX::AdditionalData& data) + { + const unsigned int min = m.get_minlevel(); + const unsigned int max = m.get_maxlevel(); + + smoothers.resize(min, max); + + for (unsigned int i=min;i<=max;++i) + smoothers[i].initialize(m[i], data); + } + + + template + template + inline void + SmootherRelaxation::initialize ( + const MGLevelObject& m, + const MGLevelObject& data) + { + const unsigned int min = m.get_minlevel(); + const unsigned int max = m.get_maxlevel(); + + Assert (data.get_minlevel() == min, + ExcDimensionMismatch(data.get_minlevel(), min)); + Assert (data.get_maxlevel() == max, + ExcDimensionMismatch(data.get_maxlevel(), max)); + + smoothers.resize(min, max); + + for (unsigned int i=min;i<=max;++i) + smoothers[i].initialize(m[i], data[i]); + } + + template + inline void + SmootherRelaxation::smooth( + const unsigned int level, + VECTOR& u, + const VECTOR& rhs) const + { + unsigned int maxlevel = smoothers.get_maxlevel(); + unsigned int steps2 = this->steps; + + if (this->variable) + steps2 *= (1<<(maxlevel-level)); + + bool T = this->transpose; + if (this->symmetric && (steps2 % 2 == 0)) + T = false; + if (this->debug > 0) + deallog << 'S' << level << ' '; + + for (unsigned int i=0; isymmetric) + T = !T; + } + } + + + template + inline unsigned int + SmootherRelaxation:: + memory_consumption () const + { + return sizeof(*this) + + smoothers.memory_consumption() + + this->mem->memory_consumption(); + } +} + + //----------------------------------------------------------------------// template @@ -822,7 +1038,7 @@ MGSmootherRelaxation::smooth( VECTOR& u, const VECTOR& rhs) const { - unsigned int maxlevel = matrices.get_maxlevel(); + unsigned int maxlevel = smoothers.get_maxlevel(); unsigned int steps2 = this->steps; if (this->variable)