dune-istl
2.2.0
|
00001 #ifndef DUNE_AMGSMOOTHER_HH 00002 #define DUNE_AMGSMOOTHER_HH 00003 00004 #include<dune/istl/paamg/construction.hh> 00005 #include<dune/istl/paamg/aggregates.hh> 00006 #include<dune/istl/preconditioners.hh> 00007 #include<dune/istl/schwarz.hh> 00008 #include<dune/istl/novlpschwarz.hh> 00009 #include<dune/common/propertymap.hh> 00010 00011 namespace Dune 00012 { 00013 namespace Amg 00014 { 00015 00031 template<class T> 00032 struct DefaultSmootherArgs 00033 { 00037 typedef T RelaxationFactor; 00038 00042 int iterations; 00046 RelaxationFactor relaxationFactor; 00047 00051 DefaultSmootherArgs() 00052 : iterations(1), relaxationFactor(1.0) 00053 {} 00054 }; 00055 00059 template<class T> 00060 struct SmootherTraits 00061 { 00062 typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments; 00063 00064 }; 00065 00066 template<class X, class Y, class C, class T> 00067 struct SmootherTraits<BlockPreconditioner<X,Y,C,T> > 00068 { 00069 typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments; 00070 00071 }; 00072 00073 template<class C, class T> 00074 struct SmootherTraits<NonoverlappingBlockPreconditioner<C,T> > 00075 { 00076 typedef DefaultSmootherArgs<typename T::matrix_type::field_type> Arguments; 00077 00078 }; 00079 00083 template<class T> 00084 class DefaultConstructionArgs 00085 { 00086 typedef typename T::matrix_type Matrix; 00087 00088 typedef typename SmootherTraits<T>::Arguments SmootherArgs; 00089 00090 typedef Dune::Amg::AggregatesMap<typename MatrixGraph<Matrix>::VertexDescriptor> AggregatesMap; 00091 00092 public: 00093 virtual ~DefaultConstructionArgs() 00094 {} 00095 00096 void setMatrix(const Matrix& matrix) 00097 { 00098 matrix_=&matrix; 00099 } 00100 virtual void setMatrix(const Matrix& matrix, const AggregatesMap& amap) 00101 { 00102 setMatrix(matrix); 00103 } 00104 00105 00106 const Matrix& getMatrix() const 00107 { 00108 return *matrix_; 00109 } 00110 00111 void setArgs(const SmootherArgs& args) 00112 { 00113 args_=&args; 00114 } 00115 00116 template<class T1> 00117 void setComm(T1& comm) 00118 {} 00119 00120 const SmootherArgs getArgs() const 00121 { 00122 return *args_; 00123 } 00124 00125 protected: 00126 const Matrix* matrix_; 00127 private: 00128 const SmootherArgs* args_; 00129 }; 00130 00131 template<class T> 00132 struct ConstructionArgs 00133 : public DefaultConstructionArgs<T> 00134 {}; 00135 00136 template<class T, class C=SequentialInformation> 00137 class DefaultParallelConstructionArgs 00138 : public ConstructionArgs<T> 00139 { 00140 public: 00141 virtual ~DefaultParallelConstructionArgs() 00142 {} 00143 00144 void setComm(const C& comm) 00145 { 00146 comm_ = &comm; 00147 } 00148 00149 const C& getComm() const 00150 { 00151 return *comm_; 00152 } 00153 private: 00154 const C* comm_; 00155 }; 00156 00157 00158 template<class T> 00159 class ConstructionTraits; 00160 00164 template<class M, class X, class Y, int l> 00165 struct ConstructionTraits<SeqSSOR<M,X,Y,l> > 00166 { 00167 typedef DefaultConstructionArgs<SeqSSOR<M,X,Y,l> > Arguments; 00168 00169 static inline SeqSSOR<M,X,Y,l>* construct(Arguments& args) 00170 { 00171 return new SeqSSOR<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations, 00172 args.getArgs().relaxationFactor); 00173 } 00174 00175 static inline void deconstruct(SeqSSOR<M,X,Y,l>* ssor) 00176 { 00177 delete ssor; 00178 } 00179 00180 }; 00181 00182 00186 template<class M, class X, class Y, int l> 00187 struct ConstructionTraits<SeqSOR<M,X,Y,l> > 00188 { 00189 typedef DefaultConstructionArgs<SeqSOR<M,X,Y,l> > Arguments; 00190 00191 static inline SeqSOR<M,X,Y,l>* construct(Arguments& args) 00192 { 00193 return new SeqSOR<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations, 00194 args.getArgs().relaxationFactor); 00195 } 00196 00197 static inline void deconstruct(SeqSOR<M,X,Y,l>* sor) 00198 { 00199 delete sor; 00200 } 00201 00202 }; 00206 template<class M, class X, class Y, int l> 00207 struct ConstructionTraits<SeqJac<M,X,Y,l> > 00208 { 00209 typedef DefaultConstructionArgs<SeqJac<M,X,Y,l> > Arguments; 00210 00211 static inline SeqJac<M,X,Y,l>* construct(Arguments& args) 00212 { 00213 return new SeqJac<M,X,Y,l>(args.getMatrix(), args.getArgs().iterations, 00214 args.getArgs().relaxationFactor); 00215 } 00216 00217 static void deconstruct(SeqJac<M,X,Y,l>* jac) 00218 { 00219 delete jac; 00220 } 00221 00222 }; 00223 00224 00228 template<class M, class X, class Y> 00229 struct ConstructionTraits<SeqILU0<M,X,Y> > 00230 { 00231 typedef DefaultConstructionArgs<SeqILU0<M,X,Y> > Arguments; 00232 00233 static inline SeqILU0<M,X,Y>* construct(Arguments& args) 00234 { 00235 return new SeqILU0<M,X,Y>(args.getMatrix(), 00236 args.getArgs().relaxationFactor); 00237 } 00238 00239 static void deconstruct(SeqILU0<M,X,Y>* ilu) 00240 { 00241 delete ilu; 00242 } 00243 00244 }; 00245 00246 template<class M, class X, class Y> 00247 class ConstructionArgs<SeqILUn<M,X,Y> > 00248 : public DefaultConstructionArgs<SeqILUn<M,X,Y> > 00249 { 00250 public: 00251 ConstructionArgs(int n=1) 00252 : n_(n) 00253 {} 00254 00255 void setN(int n) 00256 { 00257 n_ = n; 00258 } 00259 int getN() 00260 { 00261 return n_; 00262 } 00263 00264 private: 00265 int n_; 00266 }; 00267 00268 00272 template<class M, class X, class Y> 00273 struct ConstructionTraits<SeqILUn<M,X,Y> > 00274 { 00275 typedef ConstructionArgs<SeqILUn<M,X,Y> > Arguments; 00276 00277 static inline SeqILUn<M,X,Y>* construct(Arguments& args) 00278 { 00279 return new SeqILUn<M,X,Y>(args.getMatrix(), args.getN(), 00280 args.getArgs().relaxationFactor); 00281 } 00282 00283 static void deconstruct(SeqILUn<M,X,Y>* ilu) 00284 { 00285 delete ilu; 00286 } 00287 00288 }; 00289 00293 template<class M, class X, class Y, class C> 00294 struct ConstructionTraits<ParSSOR<M,X,Y,C> > 00295 { 00296 typedef DefaultParallelConstructionArgs<M,C> Arguments; 00297 00298 static inline ParSSOR<M,X,Y,C>* construct(Arguments& args) 00299 { 00300 return new ParSSOR<M,X,Y,C>(args.getMatrix(), args.getArgs().iterations, 00301 args.getArgs().relaxationFactor, 00302 args.getComm()); 00303 } 00304 static inline void deconstruct(ParSSOR<M,X,Y,C>* ssor) 00305 { 00306 delete ssor; 00307 } 00308 }; 00309 00310 template<class X, class Y, class C, class T> 00311 struct ConstructionTraits<BlockPreconditioner<X,Y,C,T> > 00312 { 00313 typedef DefaultParallelConstructionArgs<T,C> Arguments; 00314 typedef ConstructionTraits<T> SeqConstructionTraits; 00315 static inline BlockPreconditioner<X,Y,C,T>* construct(Arguments& args) 00316 { 00317 return new BlockPreconditioner<X,Y,C,T>(*SeqConstructionTraits::construct(args), 00318 args.getComm()); 00319 } 00320 00321 static inline void deconstruct(BlockPreconditioner<X,Y,C,T>* bp) 00322 { 00323 SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner)); 00324 delete bp; 00325 } 00326 00327 }; 00328 00329 template<class C, class T> 00330 struct ConstructionTraits<NonoverlappingBlockPreconditioner<C,T> > 00331 { 00332 typedef DefaultParallelConstructionArgs<T,C> Arguments; 00333 typedef ConstructionTraits<T> SeqConstructionTraits; 00334 static inline NonoverlappingBlockPreconditioner<C,T>* construct(Arguments& args) 00335 { 00336 return new NonoverlappingBlockPreconditioner<C,T>(*SeqConstructionTraits::construct(args), 00337 args.getComm()); 00338 } 00339 00340 static inline void deconstruct(NonoverlappingBlockPreconditioner<C,T>* bp) 00341 { 00342 SeqConstructionTraits::deconstruct(static_cast<T*>(&bp->preconditioner)); 00343 delete bp; 00344 } 00345 00346 }; 00347 00358 template<class T> 00359 struct SmootherApplier 00360 { 00361 typedef T Smoother; 00362 typedef typename Smoother::range_type Range; 00363 typedef typename Smoother::domain_type Domain; 00364 00372 static void preSmooth(Smoother& smoother, Domain& v, const Range& d) 00373 { 00374 smoother.apply(v,d); 00375 } 00376 00384 static void postSmooth(Smoother& smoother, Domain& v, const Range& d) 00385 { 00386 smoother.apply(v,d); 00387 } 00388 }; 00389 00390 template<class M, class X, class Y, int l> 00391 struct SmootherApplier<SeqSOR<M,X,Y,l> > 00392 { 00393 typedef SeqSOR<M,X,Y,l> Smoother; 00394 typedef typename Smoother::range_type Range; 00395 typedef typename Smoother::domain_type Domain; 00396 00397 static void preSmooth(Smoother& smoother, Domain& v, Range& d) 00398 { 00399 smoother.template apply<true>(v,d); 00400 } 00401 00402 00403 static void postSmooth(Smoother& smoother, Domain& v, Range& d) 00404 { 00405 smoother.template apply<false>(v,d); 00406 } 00407 }; 00408 00409 template<class M, class X, class Y, class C, int l> 00410 struct SmootherApplier<BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > > 00411 { 00412 typedef BlockPreconditioner<X,Y,C,SeqSOR<M,X,Y,l> > Smoother; 00413 typedef typename Smoother::range_type Range; 00414 typedef typename Smoother::domain_type Domain; 00415 00416 static void preSmooth(Smoother& smoother, Domain& v, Range& d) 00417 { 00418 smoother.template apply<true>(v,d); 00419 } 00420 00421 00422 static void postSmooth(Smoother& smoother, Domain& v, Range& d) 00423 { 00424 smoother.template apply<false>(v,d); 00425 } 00426 }; 00427 00428 template<class M, class X, class Y, class C, int l> 00429 struct SmootherApplier<NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > > 00430 { 00431 typedef NonoverlappingBlockPreconditioner<C,SeqSOR<M,X,Y,l> > Smoother; 00432 typedef typename Smoother::range_type Range; 00433 typedef typename Smoother::domain_type Domain; 00434 00435 static void preSmooth(Smoother& smoother, Domain& v, Range& d) 00436 { 00437 smoother.template apply<true>(v,d); 00438 } 00439 00440 00441 static void postSmooth(Smoother& smoother, Domain& v, Range& d) 00442 { 00443 smoother.template apply<false>(v,d); 00444 } 00445 }; 00446 00447 }// end namespace Amg 00448 00449 // forward declarations 00450 template<class M, class X, class MO, class MS, class A> 00451 class SeqOverlappingSchwarz; 00452 00453 class MultiplicativeSchwarzMode; 00454 00455 namespace Amg 00456 { 00457 template<class M, class X, class MS, class TA> 00458 struct SmootherApplier<SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode, 00459 MS,TA> > 00460 { 00461 typedef SeqOverlappingSchwarz<M,X,MultiplicativeSchwarzMode,MS,TA> Smoother; 00462 typedef typename Smoother::range_type Range; 00463 typedef typename Smoother::domain_type Domain; 00464 00465 static void preSmooth(Smoother& smoother, Domain& v, const Range& d) 00466 { 00467 smoother.template apply<true>(v,d); 00468 } 00469 00470 00471 static void postSmooth(Smoother& smoother, Domain& v, const Range& d) 00472 { 00473 smoother.template apply<false>(v,d); 00474 00475 } 00476 }; 00477 00478 // template<class M, class X, class TM, class TA> 00479 // class SeqOverlappingSchwarz; 00480 00481 template<class T> 00482 struct SeqOverlappingSchwarzSmootherArgs 00483 : public DefaultSmootherArgs<T> 00484 { 00485 enum Overlap {vertex, aggregate, pairwise, none}; 00486 00487 Overlap overlap; 00488 bool onthefly; 00489 00490 SeqOverlappingSchwarzSmootherArgs(Overlap overlap_=vertex, 00491 bool onthefly_=false) 00492 : overlap(overlap_), onthefly(onthefly_) 00493 {} 00494 }; 00495 00496 template<class M, class X, class TM, class TS, class TA> 00497 struct SmootherTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> > 00498 { 00499 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> Arguments; 00500 }; 00501 00502 template<class M, class X, class TM, class TS, class TA> 00503 class ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > 00504 : public DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > 00505 { 00506 typedef DefaultConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Father; 00507 00508 public: 00509 typedef typename MatrixGraph<M>::VertexDescriptor VertexDescriptor; 00510 typedef Dune::Amg::AggregatesMap<VertexDescriptor> AggregatesMap; 00511 typedef typename AggregatesMap::AggregateDescriptor AggregateDescriptor; 00512 typedef typename SeqOverlappingSchwarz<M,X,TM,TS,TA>::subdomain_vector Vector; 00513 typedef typename Vector::value_type Subdomain; 00514 00515 virtual void setMatrix(const M& matrix, const AggregatesMap& amap) 00516 { 00517 Father::setMatrix(matrix); 00518 00519 std::vector<bool> visited(amap.noVertices(), false); 00520 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType; 00521 VisitedMapType visitedMap(visited.begin()); 00522 00523 MatrixGraph<const M> graph(matrix); 00524 00525 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs; 00526 00527 switch(Father::getArgs().overlap){ 00528 case SmootherArgs::vertex: 00529 { 00530 VertexAdder visitor(subdomains, amap); 00531 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00532 } 00533 break; 00534 case SmootherArgs::pairwise: 00535 { 00536 createPairDomains(graph); 00537 } 00538 break; 00539 case SmootherArgs::aggregate: 00540 { 00541 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap); 00542 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00543 } 00544 break; 00545 case SmootherArgs::none: 00546 NoneAdder visitor; 00547 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00548 break; 00549 default: 00550 DUNE_THROW(NotImplemented, "This overlapping scheme is not supported!"); 00551 } 00552 } 00553 void setMatrix(const M& matrix) 00554 { 00555 Father::setMatrix(matrix); 00556 00557 /* Create aggregates map where each aggregate is just one vertex. */ 00558 AggregatesMap amap(matrix.N()); 00559 VertexDescriptor v=0; 00560 for(typename AggregatesMap::iterator iter=amap.begin(); 00561 iter!=amap.end(); ++iter) 00562 *iter=v++; 00563 00564 std::vector<bool> visited(amap.noVertices(), false); 00565 typedef IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap> VisitedMapType; 00566 VisitedMapType visitedMap(visited.begin()); 00567 00568 MatrixGraph<const M> graph(matrix); 00569 00570 typedef SeqOverlappingSchwarzSmootherArgs<typename M::field_type> SmootherArgs; 00571 00572 switch(Father::getArgs().overlap){ 00573 case SmootherArgs::vertex: 00574 { 00575 VertexAdder visitor(subdomains, amap); 00576 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00577 } 00578 break; 00579 case SmootherArgs::aggregate: 00580 { 00581 DUNE_THROW(NotImplemented, "Aggregate overlap is not supported yet"); 00582 /* 00583 AggregateAdder<VisitedMapType> visitor(subdomains, amap, graph, visitedMap); 00584 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00585 */ 00586 } 00587 break; 00588 case SmootherArgs::pairwise: 00589 { 00590 createPairDomains(graph); 00591 } 00592 break; 00593 case SmootherArgs::none: 00594 NoneAdder visitor; 00595 createSubdomains(matrix, graph, amap, visitor, visitedMap); 00596 00597 } 00598 } 00599 00600 const Vector& getSubDomains() 00601 { 00602 return subdomains; 00603 } 00604 00605 private: 00606 struct VertexAdder 00607 { 00608 VertexAdder(Vector& subdomains_, const AggregatesMap& aggregates_) 00609 : subdomains(subdomains_), max(-1), subdomain(-1), aggregates(aggregates_) 00610 {} 00611 template<class T> 00612 void operator()(const T& edge) 00613 { 00614 if(aggregates[edge.target()]!=AggregatesMap::ISOLATED) 00615 subdomains[subdomain].insert(edge.target()); 00616 } 00617 int setAggregate(const AggregateDescriptor& aggregate_) 00618 { 00619 subdomain=aggregate_; 00620 max = std::max(subdomain, aggregate_); 00621 return subdomain; 00622 } 00623 int noSubdomains() const 00624 { 00625 return max+1; 00626 } 00627 private: 00628 Vector& subdomains; 00629 AggregateDescriptor max; 00630 AggregateDescriptor subdomain; 00631 const AggregatesMap& aggregates; 00632 }; 00633 struct NoneAdder 00634 { 00635 template<class T> 00636 void operator()(const T& edge) 00637 {} 00638 int setAggregate(const AggregateDescriptor& aggregate_) 00639 { 00640 return -1; 00641 } 00642 int noSubdomains() const 00643 { 00644 return -1; 00645 } 00646 }; 00647 00648 template<class VM> 00649 struct AggregateAdder 00650 { 00651 AggregateAdder(Vector& subdomains_, const AggregatesMap& aggregates_, 00652 const MatrixGraph<const M>& graph_, VM& visitedMap_) 00653 : subdomains(subdomains_), subdomain(-1), aggregates(aggregates_), 00654 adder(subdomains_, aggregates_), graph(graph_), visitedMap(visitedMap_) 00655 {} 00656 template<class T> 00657 void operator()(const T& edge) 00658 { 00659 subdomains[subdomain].insert(edge.target()); 00660 // If we (the neighbouring vertex of the aggregate) 00661 // are not isolated, add the aggregate we belong to 00662 // to the same subdomain using the OneOverlapAdder 00663 if(aggregates[edge.target()]!=AggregatesMap::ISOLATED){ 00664 assert(aggregates[edge.target()]!=aggregate); 00665 typename AggregatesMap::VertexList vlist; 00666 aggregates.template breadthFirstSearch<true,false>(edge.target(), aggregate, 00667 graph, vlist, adder, adder, 00668 visitedMap); 00669 } 00670 } 00671 00672 int setAggregate(const AggregateDescriptor& aggregate_) 00673 { 00674 adder.setAggregate(aggregate_); 00675 aggregate=aggregate_; 00676 return ++subdomain; 00677 } 00678 int noSubdomains() const 00679 { 00680 return subdomain+1; 00681 } 00682 00683 private: 00684 AggregateDescriptor aggregate; 00685 Vector& subdomains; 00686 int subdomain; 00687 const AggregatesMap& aggregates; 00688 VertexAdder adder; 00689 const MatrixGraph<const M>& graph; 00690 VM& visitedMap; 00691 }; 00692 00693 void createPairDomains(const MatrixGraph<const M>& graph) 00694 { 00695 typedef typename MatrixGraph<const M>::ConstVertexIterator VIter; 00696 typedef typename MatrixGraph<const M>::ConstEdgeIterator EIter; 00697 typedef typename M::size_type size_type; 00698 00699 std::set<std::pair<size_type,size_type> > pairs; 00700 int total=0; 00701 for(VIter v=graph.begin(), ve=graph.end(); ve != v; ++v) 00702 for(EIter e = v.begin(), ee=v.end(); ee!=e; ++e) 00703 { 00704 ++total; 00705 if(e.source()<e.target()) 00706 pairs.insert(std::make_pair(e.source(),e.target())); 00707 else 00708 pairs.insert(std::make_pair(e.target(),e.source())); 00709 } 00710 00711 00712 subdomains.resize(pairs.size()); 00713 Dune::dinfo <<std::endl<< "Created "<<pairs.size()<<" ("<<total<<") pair domains"<<std::endl<<std::endl; 00714 typedef typename std::set<std::pair<size_type,size_type> >::const_iterator SIter; 00715 typename Vector::iterator subdomain=subdomains.begin(); 00716 00717 for(SIter s=pairs.begin(), se =pairs.end(); se!=s; ++s) 00718 { 00719 subdomain->insert(s->first); 00720 subdomain->insert(s->second); 00721 ++subdomain; 00722 } 00723 std::size_t minsize=10000; 00724 std::size_t maxsize=0; 00725 int sum=0; 00726 for(typename Vector::size_type i=0; i < subdomains.size(); ++i){ 00727 sum+=subdomains[i].size(); 00728 minsize=std::min(minsize, subdomains[i].size()); 00729 maxsize=std::max(maxsize, subdomains[i].size()); 00730 } 00731 Dune::dinfo<<"Subdomain size: min="<<minsize<<" max="<<maxsize<<" avg="<<(sum/subdomains.size()) 00732 <<" no="<<subdomains.size()<<std::endl; 00733 } 00734 00735 template<class Visitor> 00736 void createSubdomains(const M& matrix, const MatrixGraph<const M>& graph, 00737 const AggregatesMap& amap, Visitor& overlapVisitor, 00738 IteratorPropertyMap<std::vector<bool>::iterator,IdentityMap>& visitedMap ) 00739 { 00740 // count number ag aggregates. We asume that the 00741 // aggregates are numbered consecutively from 0 exept 00742 // for the isolated ones. All isolated vertices form 00743 // one aggregate, here. 00744 int isolated=0; 00745 AggregateDescriptor maxAggregate=0; 00746 00747 for(std::size_t i=0; i < amap.noVertices(); ++i) 00748 if(amap[i]==AggregatesMap::ISOLATED) 00749 isolated++; 00750 else 00751 maxAggregate = std::max(maxAggregate, amap[i]); 00752 00753 subdomains.resize(maxAggregate+1+isolated); 00754 00755 // reset the subdomains 00756 for(typename Vector::size_type i=0; i < subdomains.size(); ++i) 00757 subdomains[i].clear(); 00758 00759 // Create the subdomains from the aggregates mapping. 00760 // For each aggregate we mark all entries and the 00761 // neighbouring vertices as belonging to the same subdomain 00762 VertexAdder aggregateVisitor(subdomains, amap); 00763 00764 for(VertexDescriptor i=0; i < amap.noVertices(); ++i) 00765 if(!get(visitedMap, i)){ 00766 AggregateDescriptor aggregate=amap[i]; 00767 00768 if(amap[i]==AggregatesMap::ISOLATED){ 00769 // isolated vertex gets its own aggregate 00770 subdomains.push_back(Subdomain()); 00771 aggregate=subdomains.size()-1; 00772 } 00773 overlapVisitor.setAggregate(aggregate); 00774 aggregateVisitor.setAggregate(aggregate); 00775 subdomains[aggregate].insert(i); 00776 typename AggregatesMap::VertexList vlist; 00777 amap.template breadthFirstSearch<false,false>(i, aggregate, graph, vlist, aggregateVisitor, 00778 overlapVisitor, visitedMap); 00779 } 00780 00781 std::size_t minsize=10000; 00782 std::size_t maxsize=0; 00783 int sum=0; 00784 for(typename Vector::size_type i=0; i < subdomains.size(); ++i){ 00785 sum+=subdomains[i].size(); 00786 minsize=std::min(minsize, subdomains[i].size()); 00787 maxsize=std::max(maxsize, subdomains[i].size()); 00788 } 00789 Dune::dinfo<<"Subdomain size: min="<<minsize<<" max="<<maxsize<<" avg="<<(sum/subdomains.size()) 00790 <<" no="<<subdomains.size()<<" isolated="<<isolated<<std::endl; 00791 00792 00793 00794 } 00795 Vector subdomains; 00796 }; 00797 00798 00799 template<class M, class X, class TM, class TS, class TA> 00800 struct ConstructionTraits<SeqOverlappingSchwarz<M,X,TM,TS,TA> > 00801 { 00802 typedef ConstructionArgs<SeqOverlappingSchwarz<M,X,TM,TS,TA> > Arguments; 00803 00804 static inline SeqOverlappingSchwarz<M,X,TM,TS,TA>* construct(Arguments& args) 00805 { 00806 return new SeqOverlappingSchwarz<M,X,TM,TS,TA>(args.getMatrix(), 00807 args.getSubDomains(), 00808 args.getArgs().relaxationFactor, 00809 args.getArgs().onthefly); 00810 } 00811 00812 static void deconstruct(SeqOverlappingSchwarz<M,X,TM,TS,TA>* schwarz) 00813 { 00814 delete schwarz; 00815 } 00816 }; 00817 00818 00819 } // namespace Amg 00820 } // namespace Dune 00821 00822 00823 00824 #endif