dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=4 sw=2 sts=2: 00003 #ifndef DUNE_SCHWARZ_HH 00004 #define DUNE_SCHWARZ_HH 00005 00006 #include <iostream> // for input/output to shell 00007 #include <fstream> // for input/output to files 00008 #include <vector> // STL vector class 00009 #include <sstream> 00010 00011 #include <cmath> // Yes, we do some math here 00012 00013 #include <dune/common/timer.hh> 00014 00015 #include"io.hh" 00016 #include"bvector.hh" 00017 #include"vbvector.hh" 00018 #include"bcrsmatrix.hh" 00019 #include"io.hh" 00020 #include"gsetc.hh" 00021 #include"ilu.hh" 00022 #include"operators.hh" 00023 #include"solvers.hh" 00024 #include"preconditioners.hh" 00025 #include"scalarproducts.hh" 00026 #include"owneroverlapcopy.hh" 00027 00028 namespace Dune { 00029 00058 template<class M, class X, class Y, class C> 00059 class OverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y> 00060 { 00061 public: 00063 typedef M matrix_type; 00065 typedef X domain_type; 00067 typedef Y range_type; 00069 typedef typename X::field_type field_type; 00071 typedef C communication_type; 00072 00073 enum { 00075 category=SolverCategory::overlapping 00076 }; 00077 00085 OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com) 00086 : _A_(A), communication(com) 00087 {} 00088 00090 virtual void apply (const X& x, Y& y) const 00091 { 00092 y = 0; 00093 _A_.umv(x,y); // result is consistent on interior+border 00094 communication.project(y); // we want this here to avoid it before the preconditioner 00095 // since there d is const! 00096 } 00097 00099 virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const 00100 { 00101 _A_.usmv(alpha,x,y); // result is consistent on interior+border 00102 communication.project(y); // we want this here to avoid it before the preconditioner 00103 // since there d is const! 00104 } 00105 00107 virtual const matrix_type& getmat () const 00108 { 00109 return _A_; 00110 } 00111 00112 private: 00113 const matrix_type& _A_; 00114 const communication_type& communication; 00115 }; 00116 00128 template<class X, class C> 00129 class OverlappingSchwarzScalarProduct : public ScalarProduct<X> 00130 { 00131 public: 00133 typedef X domain_type; 00135 typedef typename X::field_type field_type; 00137 typedef C communication_type; 00138 00140 enum {category=SolverCategory::overlapping}; 00141 00146 OverlappingSchwarzScalarProduct (const communication_type& com) 00147 : communication(com) 00148 {} 00149 00154 virtual field_type dot (const X& x, const X& y) 00155 { 00156 field_type result; 00157 communication.dot(x,y,result); 00158 return result; 00159 } 00160 00164 virtual double norm (const X& x) 00165 { 00166 return communication.norm(x); 00167 } 00168 00169 private: 00170 const communication_type& communication; 00171 }; 00172 00173 template<class X, class C> 00174 struct ScalarProductChooser<X,C,SolverCategory::overlapping> 00175 { 00177 typedef OverlappingSchwarzScalarProduct<X,C> ScalarProduct; 00179 typedef C communication_type; 00180 00181 enum{ 00183 solverCategory=SolverCategory::overlapping 00184 }; 00185 00186 static ScalarProduct* construct(const communication_type& comm) 00187 { 00188 return new ScalarProduct(comm); 00189 } 00190 }; 00191 00198 00199 template<class M, class X, class Y, class C> 00200 class ParSSOR : public Preconditioner<X,Y> { 00201 public: 00203 typedef M matrix_type; 00205 typedef X domain_type; 00207 typedef Y range_type; 00209 typedef typename X::field_type field_type; 00211 typedef C communication_type; 00212 00213 // define the category 00214 enum { 00216 category=SolverCategory::overlapping}; 00217 00227 ParSSOR (const matrix_type& A, int n, field_type w, const communication_type& c) 00228 : _A_(A), _n(n), _w(w), communication(c) 00229 { } 00230 00236 virtual void pre (X& x, Y& b) 00237 { 00238 communication.copyOwnerToAll(x,x); // make dirichlet values consistent 00239 } 00240 00246 virtual void apply (X& v, const Y& d) 00247 { 00248 for (int i=0; i<_n; i++){ 00249 bsorf(_A_,v,d,_w); 00250 bsorb(_A_,v,d,_w); 00251 } 00252 communication.copyOwnerToAll(v,v); 00253 } 00254 00260 virtual void post (X& x) {} 00261 00262 private: 00264 const matrix_type& _A_; 00266 int _n; 00268 field_type _w; 00270 const communication_type& communication; 00271 }; 00272 00273 namespace Amg 00274 { 00275 template<class T> class ConstructionTraits; 00276 } 00277 00286 template<class X, class Y, class C, class T=Preconditioner<X,Y> > 00287 class BlockPreconditioner : public Preconditioner<X,Y> { 00288 friend class Amg::ConstructionTraits<BlockPreconditioner<X,Y,C,T> >; 00289 public: 00291 typedef X domain_type; 00293 typedef Y range_type; 00295 typedef typename X::field_type field_type; 00297 typedef C communication_type; 00298 00299 // define the category 00300 enum { 00302 category=SolverCategory::overlapping}; 00303 00311 BlockPreconditioner (T& p, const communication_type& c) 00312 : preconditioner(p), communication(c) 00313 { } 00314 00320 virtual void pre (X& x, Y& b) 00321 { 00322 communication.copyOwnerToAll(x,x); // make dirichlet values consistent 00323 preconditioner.pre(x,b); 00324 } 00325 00331 virtual void apply (X& v, const Y& d) 00332 { 00333 preconditioner.apply(v,d); 00334 communication.copyOwnerToAll(v,v); 00335 } 00336 00337 template<bool forward> 00338 void apply (X& v, const Y& d) 00339 { 00340 preconditioner.template apply<forward>(v,d); 00341 communication.copyOwnerToAll(v,v); 00342 } 00343 00349 virtual void post (X& x) 00350 { 00351 preconditioner.post(x); 00352 } 00353 00354 private: 00356 T& preconditioner; 00357 00359 const communication_type& communication; 00360 }; 00361 00364 } // end namespace 00365 00366 #endif