dune-istl
2.2.0
|
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