include/deal.II/multigrid/mg_block_smoother.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: mg_block_smoother.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2005, 2006, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
00010 //
00011 //---------------------------------------------------------------------------
00012 #ifndef __deal2__mg_block_smoother_h
00013 #define __deal2__mg_block_smoother_h
00014 
00015 
00016 #include <deal.II/base/config.h>
00017 #include <deal.II/base/smartpointer.h>
00018 #include <deal.II/lac/pointer_matrix.h>
00019 #include <deal.II/lac/vector_memory.h>
00020 #include <deal.II/lac/block_vector.h>
00021 #include <deal.II/multigrid/mg_base.h>
00022 #include <deal.II/base/mg_level_object.h>
00023 #include <vector>
00024 
00025 DEAL_II_NAMESPACE_OPEN
00026 
00027 template <int,int> class MGDoFHandler;
00028 
00029 
00030 /*
00031  * MGSmootherBase is defined in mg_base.h
00032  */
00033 
00036 
00045 template <class MATRIX, class RELAX, typename number>
00046 class MGSmootherBlock
00047   : public MGSmootherBase<BlockVector<number> >
00048 {
00049   public:
00054     MGSmootherBlock(VectorMemory<BlockVector<number> >& mem,
00055                     const unsigned int steps = 1,
00056                     const bool variable = false,
00057                     const bool symmetric = false,
00058                     const bool transpose = false,
00059                     const bool reverse = false);
00060 
00083     template <class MGMATRIX, class MGRELAX>
00084     void initialize (const MGMATRIX& matrices,
00085                      const MGRELAX& smoothers);
00086 
00090     void clear ();
00091 
00096     void set_steps (const unsigned int);
00097 
00102     void set_variable (const bool);
00103 
00108     void set_symmetric (const bool);
00109 
00115     void set_transpose (const bool);
00116 
00122     void set_reverse (const bool);
00123 
00134     virtual void smooth (const unsigned int         level,
00135                          BlockVector<number>&       u,
00136                          const BlockVector<number>& rhs) const;
00137  private:
00141     MGLevelObject<PointerMatrix<MATRIX, BlockVector<number> > > matrices;
00142 
00146     MGLevelObject<PointerMatrix<RELAX, BlockVector<number> > > smoothers;
00147 
00151     unsigned int steps;
00152 
00156     bool variable;
00157 
00161     bool symmetric;
00162 
00163                                      /*
00164                                       * Transposed?
00165                                       */
00166     bool transpose;
00167 
00171     bool reverse;
00172 
00176     VectorMemory<BlockVector<number> >& mem;
00177 
00178 };
00179 
00182 //---------------------------------------------------------------------------
00183 
00184 #ifndef DOXYGEN
00185 
00186 template <class MATRIX, class RELAX, typename number>
00187 inline
00188 MGSmootherBlock<MATRIX, RELAX, number>::MGSmootherBlock(
00189   VectorMemory<BlockVector<number> >& mem,
00190   const unsigned int steps,
00191   const bool variable,
00192   const bool symmetric,
00193   const bool transpose,
00194   const bool reverse)
00195                 :
00196                 steps(steps),
00197                 variable(variable),
00198                 symmetric(symmetric),
00199                 transpose(transpose),
00200                 reverse(reverse),
00201                 mem(mem)
00202 {}
00203 
00204 
00205 template <class MATRIX, class RELAX, typename number>
00206 inline void
00207 MGSmootherBlock<MATRIX, RELAX, number>::clear ()
00208 {
00209   unsigned int i=matrices.get_minlevel(),
00210        max_level=matrices.get_maxlevel();
00211   for (; i<=max_level; ++i)
00212     {
00213       smoothers[i] = 0;
00214       matrices[i] = 0;
00215     }
00216 }
00217 
00218 
00219 template <class MATRIX, class RELAX, typename number>
00220 template <class MGMATRIX, class MGRELAX>
00221 inline void
00222 MGSmootherBlock<MATRIX, RELAX, number>::initialize (
00223   const MGMATRIX& m,
00224   const MGRELAX& s)
00225 {
00226   const unsigned int min = m.get_minlevel();
00227   const unsigned int max = m.get_maxlevel();
00228 
00229   matrices.resize(min, max);
00230   smoothers.resize(min, max);
00231 
00232   for (unsigned int i=min;i<=max;++i)
00233     {
00234       matrices[i] = &m[i];
00235       smoothers[i] = &s[i];
00236     }
00237 }
00238 
00239 template <class MATRIX, class RELAX, typename number>
00240 inline void
00241 MGSmootherBlock<MATRIX, RELAX, number>::
00242 set_steps (const unsigned int s)
00243 {
00244   steps = s;
00245 }
00246 
00247 
00248 template <class MATRIX, class RELAX, typename number>
00249 inline void
00250 MGSmootherBlock<MATRIX, RELAX, number>::
00251 set_variable (const bool flag)
00252 {
00253   variable = flag;
00254 }
00255 
00256 
00257 template <class MATRIX, class RELAX, typename number>
00258 inline void
00259 MGSmootherBlock<MATRIX, RELAX, number>::
00260 set_symmetric (const bool flag)
00261 {
00262   symmetric = flag;
00263 }
00264 
00265 
00266 template <class MATRIX, class RELAX, typename number>
00267 inline void
00268 MGSmootherBlock<MATRIX, RELAX, number>::
00269 set_transpose (const bool flag)
00270 {
00271   transpose = flag;
00272 }
00273 
00274 
00275 template <class MATRIX, class RELAX, typename number>
00276 inline void
00277 MGSmootherBlock<MATRIX, RELAX, number>::
00278 set_reverse (const bool flag)
00279 {
00280   reverse = flag;
00281 }
00282 
00283 
00284 template <class MATRIX, class RELAX, typename number>
00285 inline void
00286 MGSmootherBlock<MATRIX, RELAX, number>::smooth(
00287   const unsigned int level,
00288   BlockVector<number>& u,
00289   const BlockVector<number>& rhs) const
00290 {
00291   deallog.push("Smooth");
00292 
00293   unsigned int maxlevel = matrices.get_maxlevel();
00294   unsigned int steps2 = steps;
00295 
00296   if (variable)
00297     steps2 *= (1<<(maxlevel-level));
00298 
00299   BlockVector<number>* r = mem.alloc();
00300   BlockVector<number>* d = mem.alloc();
00301   r->reinit(u);
00302   d->reinit(u);
00303 
00304   bool T = transpose;
00305   if (symmetric && (steps2 % 2 == 0))
00306     T = false;
00307 //  cerr << 'S' << level;
00308 //  cerr << '(' << matrices[level]->m() << ',' << matrices[level]->n() << ')';
00309 
00310   for (unsigned int i=0; i<steps2; ++i)
00311     {
00312       if (T)
00313         {
00314 //        cerr << 'T';
00315           matrices[level].vmult(*r,u);
00316           r->sadd(-1.,1.,rhs);
00317           smoothers[level].Tvmult(*d, *r);
00318         } else {
00319 //        cerr << 'N';
00320           matrices[level].vmult(*r,u);
00321           r->sadd(-1.,1.,rhs);
00322           smoothers[level].vmult(*d, *r);
00323         }
00324 //      cerr << '{' << r->l2_norm() << '}';
00325       u += *d;
00326       if (symmetric)
00327         T = !T;
00328  }
00329 
00330   mem.free(r);
00331   mem.free(d);
00332   deallog.pop();
00333 }
00334 
00335 #endif // DOXYGEN
00336 
00337 DEAL_II_NAMESPACE_CLOSE
00338 
00339 #endif
00340 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Tue May 22 2012 12:06:16 by doxygen 1.7.3