00001
00002
00003
00004
00005
00006
00007
00008
00009
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
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
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
00308
00309
00310 for (unsigned int i=0; i<steps2; ++i)
00311 {
00312 if (T)
00313 {
00314
00315 matrices[level].vmult(*r,u);
00316 r->sadd(-1.,1.,rhs);
00317 smoothers[level].Tvmult(*d, *r);
00318 } else {
00319
00320 matrices[level].vmult(*r,u);
00321 r->sadd(-1.,1.,rhs);
00322 smoothers[level].vmult(*d, *r);
00323 }
00324
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