dune-istl  2.2.0
globalaggregates.hh
Go to the documentation of this file.
00001 #ifndef DUNE_GLOBALAGGREGATES_HH
00002 #define DUNE_GLOBALAGGREGATES_HH
00003 
00004 #include "aggregates.hh"
00005 #include "pinfo.hh"
00006 #include <dune/common/parallel/indexset.hh>
00007 
00008 namespace Dune
00009 {
00010   namespace Amg
00011   {
00012     
00013     template<typename T, typename TI>
00014     struct GlobalAggregatesMap
00015     {
00016     public:
00017       typedef TI ParallelIndexSet;
00018       
00019       typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00020       
00021       typedef typename ParallelIndexSet::GlobalIndex IndexedType;
00022   
00023       typedef typename ParallelIndexSet::LocalIndex LocalIndex;
00024       
00025       typedef T Vertex;
00026       
00027       GlobalAggregatesMap(AggregatesMap<Vertex>& aggregates,
00028                           const GlobalLookupIndexSet<ParallelIndexSet>& indexset)
00029         : aggregates_(aggregates), indexset_(indexset)
00030       {}
00031   
00032       inline const GlobalIndex& operator[](std::size_t index)const
00033       {
00034         const Vertex& aggregate = aggregates_[index];
00035         if(aggregate >= AggregatesMap<Vertex>::ISOLATED){
00036           assert(aggregate != AggregatesMap<Vertex>::UNAGGREGATED);
00037           return isolatedMarker;
00038         }else{
00039           const Dune::IndexPair<GlobalIndex,LocalIndex >* pair = indexset_.pair(aggregate);
00040           assert(pair!=0);
00041           return pair->global();
00042         }
00043       }
00044   
00045       
00046       inline GlobalIndex& get(std::size_t index)
00047       {
00048         const Vertex& aggregate = aggregates_[index];
00049         assert(aggregate < AggregatesMap<Vertex>::ISOLATED);
00050         const Dune::IndexPair<GlobalIndex,LocalIndex >* pair = indexset_.pair(aggregate);
00051         assert(pair!=0);
00052         return const_cast<GlobalIndex&>(pair->global());
00053       }
00054 
00055       class Proxy
00056       {
00057       public:
00058         Proxy(const GlobalLookupIndexSet<ParallelIndexSet>& indexset, Vertex& aggregate)
00059           : indexset_(&indexset), aggregate_(&aggregate)
00060         {}
00061         
00062         Proxy& operator=(const GlobalIndex& global)
00063         {
00064           if(global==isolatedMarker)
00065             *aggregate_ = AggregatesMap<Vertex>::ISOLATED;
00066           else{
00067             //assert(global < AggregatesMap<Vertex>::ISOLATED);
00068             *aggregate_ = indexset_->operator[](global).local();
00069           }
00070           return *this;
00071         }
00072       private:
00073         const GlobalLookupIndexSet<ParallelIndexSet>* indexset_;
00074         Vertex* aggregate_;
00075       };
00076       
00077       inline Proxy operator[](std::size_t index)
00078       {
00079         return Proxy(indexset_, aggregates_[index]);
00080       }
00081 
00082       inline void put(const GlobalIndex& global, size_t i)
00083       {
00084         aggregates_[i]=indexset_[global].local();
00085     
00086       }
00087 
00088     private:
00089       AggregatesMap<Vertex>& aggregates_;
00090       const GlobalLookupIndexSet<ParallelIndexSet>& indexset_;
00091       static const GlobalIndex isolatedMarker;
00092     };
00093 
00094     template<typename T, typename TI>
00095     const typename TI::GlobalIndex GlobalAggregatesMap<T,TI>::isolatedMarker = 
00096       std::numeric_limits<typename TI::GlobalIndex>::max();
00097     
00098     template<typename T, typename TI>
00099     struct AggregatesGatherScatter
00100     {
00101       typedef TI ParallelIndexSet;
00102       typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
00103       
00104       static const GlobalIndex& gather(const GlobalAggregatesMap<T,TI>& ga, size_t i)
00105       {
00106         return ga[i];
00107       }
00108   
00109       static void scatter(GlobalAggregatesMap<T,TI>& ga, GlobalIndex global, size_t i)
00110       {
00111         ga[i]=global;
00112       }
00113     };
00114 
00115     template<typename T, typename O, typename I>
00116     struct AggregatesPublisher
00117     {
00118     };
00119 
00120 #if HAVE_MPI    
00121     
00122 #endif
00123 
00124   } // namespace Amg
00125 
00126 #if HAVE_MPI  
00127     // forward declaration
00128     template<class T1, class T2>
00129     class OwnerOverlapCopyCommunication;
00130 #endif
00131 
00132   namespace Amg
00133   {
00134 
00135 #if HAVE_MPI    
00136 
00145     template<typename T, typename O, typename T1, typename T2>
00146     struct AggregatesPublisher<T,O,OwnerOverlapCopyCommunication<T1,T2> >
00147     {
00148       typedef T Vertex;
00149       typedef O OverlapFlags;
00150       typedef OwnerOverlapCopyCommunication<T1,T2> ParallelInformation;
00151       typedef typename ParallelInformation::GlobalLookupIndexSet GlobalLookupIndexSet;
00152       typedef typename ParallelInformation::ParallelIndexSet IndexSet;
00153       
00154       static void publish(AggregatesMap<Vertex>& aggregates, 
00155                           ParallelInformation& pinfo,
00156                           const GlobalLookupIndexSet& globalLookup)
00157       {
00158         typedef Dune::Amg::GlobalAggregatesMap<Vertex,IndexSet> GlobalMap;
00159         GlobalMap gmap(aggregates, globalLookup);
00160         pinfo.copyOwnerToAll(gmap,gmap);
00161         // communication only needed for ALU 
00162         // (ghosts with same global id as owners on the same process)
00163         if (pinfo.getSolverCategory() == 
00164             static_cast<int>(SolverCategory::nonoverlapping))
00165           pinfo.copyCopyToAll(gmap,gmap);
00166         
00167         typedef typename ParallelInformation::RemoteIndices::const_iterator Lists;
00168         Lists lists = pinfo.remoteIndices().find(pinfo.communicator().rank());
00169         if(lists!=pinfo.remoteIndices().end()){
00170 
00171           // For periodic boundary conditions we must renumber
00172           // the aggregates of vertices in the overlap whose owners are 
00173           // on the same process
00174           Vertex maxAggregate =0;
00175           typedef typename AggregatesMap<Vertex>::const_iterator Iter;
00176           for(Iter i=aggregates.begin(), end=aggregates.end(); i!=end; ++i)
00177             maxAggregate = std::max(maxAggregate, *i);
00178           
00179           // Compute new mapping of aggregates in the overlap that we also own
00180           std::map<Vertex,Vertex> newMapping;
00181           
00182           // insert all elements into map
00183           typedef typename ParallelInformation::RemoteIndices::RemoteIndexList
00184             ::const_iterator RIter;
00185           for(RIter ri=lists->second.first->begin(), rend = lists->second.first->end();
00186               ri!=rend; ++ri)
00187             if(O::contains(ri->localIndexPair().local().attribute()))
00188               newMapping.insert(std::make_pair(aggregates[ri->localIndexPair().local()],
00189                                                maxAggregate));
00190           // renumber
00191           typedef typename std::map<Vertex,Vertex>::iterator MIter;
00192           for(MIter mi=newMapping.begin(), mend=newMapping.end();
00193               mi != mend; ++mi)
00194             mi->second=++maxAggregate;
00195           
00196           
00197         for(RIter ri=lists->second.first->begin(), rend = lists->second.first->end();
00198               ri!=rend; ++ri)
00199           if(O::contains(ri->localIndexPair().local().attribute()))
00200               aggregates[ri->localIndexPair().local()] = 
00201                 newMapping[aggregates[ri->localIndexPair().local()]];
00202           }
00203       }
00204     };
00205 #endif
00206 
00207     template<typename T, typename O>
00208     struct AggregatesPublisher<T,O,SequentialInformation>
00209     {
00210       typedef T Vertex;
00211       typedef SequentialInformation ParallelInformation;
00212       typedef typename ParallelInformation::GlobalLookupIndexSet GlobalLookupIndexSet;
00213       
00214       static void publish(AggregatesMap<Vertex>& aggregates, 
00215                           ParallelInformation& pinfo,
00216                           const GlobalLookupIndexSet& globalLookup)
00217       {} 
00218     };
00219     
00220   } // end Amg namespace
00221 
00222 
00223 #if HAVE_MPI 
00224   template<typename T, typename TI>
00225   struct CommPolicy<Amg::GlobalAggregatesMap<T,TI> >
00226   {
00227     typedef Amg::AggregatesMap<T> Type;
00228     typedef typename Amg::GlobalAggregatesMap<T,TI>::IndexedType IndexedType;
00229     typedef SizeOne IndexedTypeFlag;
00230     static int getSize(const Type&, int)
00231     {
00232       return 1;
00233     }
00234   };
00235 #endif
00236 
00237 } // end Dune namespace
00238 
00239 #endif