dune-istl  2.2.0
globalaggregates.hh
Go to the documentation of this file.
1 #ifndef DUNE_GLOBALAGGREGATES_HH
2 #define DUNE_GLOBALAGGREGATES_HH
3 
4 #include "aggregates.hh"
5 #include "pinfo.hh"
6 #include <dune/common/parallel/indexset.hh>
7 
8 namespace Dune
9 {
10  namespace Amg
11  {
12 
13  template<typename T, typename TI>
15  {
16  public:
17  typedef TI ParallelIndexSet;
18 
19  typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
20 
21  typedef typename ParallelIndexSet::GlobalIndex IndexedType;
22 
23  typedef typename ParallelIndexSet::LocalIndex LocalIndex;
24 
25  typedef T Vertex;
26 
28  const GlobalLookupIndexSet<ParallelIndexSet>& indexset)
29  : aggregates_(aggregates), indexset_(indexset)
30  {}
31 
32  inline const GlobalIndex& operator[](std::size_t index)const
33  {
34  const Vertex& aggregate = aggregates_[index];
35  if(aggregate >= AggregatesMap<Vertex>::ISOLATED){
36  assert(aggregate != AggregatesMap<Vertex>::UNAGGREGATED);
37  return isolatedMarker;
38  }else{
39  const Dune::IndexPair<GlobalIndex,LocalIndex >* pair = indexset_.pair(aggregate);
40  assert(pair!=0);
41  return pair->global();
42  }
43  }
44 
45 
46  inline GlobalIndex& get(std::size_t index)
47  {
48  const Vertex& aggregate = aggregates_[index];
49  assert(aggregate < AggregatesMap<Vertex>::ISOLATED);
50  const Dune::IndexPair<GlobalIndex,LocalIndex >* pair = indexset_.pair(aggregate);
51  assert(pair!=0);
52  return const_cast<GlobalIndex&>(pair->global());
53  }
54 
55  class Proxy
56  {
57  public:
58  Proxy(const GlobalLookupIndexSet<ParallelIndexSet>& indexset, Vertex& aggregate)
59  : indexset_(&indexset), aggregate_(&aggregate)
60  {}
61 
62  Proxy& operator=(const GlobalIndex& global)
63  {
64  if(global==isolatedMarker)
65  *aggregate_ = AggregatesMap<Vertex>::ISOLATED;
66  else{
67  //assert(global < AggregatesMap<Vertex>::ISOLATED);
68  *aggregate_ = indexset_->operator[](global).local();
69  }
70  return *this;
71  }
72  private:
73  const GlobalLookupIndexSet<ParallelIndexSet>* indexset_;
74  Vertex* aggregate_;
75  };
76 
77  inline Proxy operator[](std::size_t index)
78  {
79  return Proxy(indexset_, aggregates_[index]);
80  }
81 
82  inline void put(const GlobalIndex& global, size_t i)
83  {
84  aggregates_[i]=indexset_[global].local();
85 
86  }
87 
88  private:
89  AggregatesMap<Vertex>& aggregates_;
90  const GlobalLookupIndexSet<ParallelIndexSet>& indexset_;
91  static const GlobalIndex isolatedMarker;
92  };
93 
94  template<typename T, typename TI>
95  const typename TI::GlobalIndex GlobalAggregatesMap<T,TI>::isolatedMarker =
96  std::numeric_limits<typename TI::GlobalIndex>::max();
97 
98  template<typename T, typename TI>
100  {
101  typedef TI ParallelIndexSet;
102  typedef typename ParallelIndexSet::GlobalIndex GlobalIndex;
103 
104  static const GlobalIndex& gather(const GlobalAggregatesMap<T,TI>& ga, size_t i)
105  {
106  return ga[i];
107  }
108 
109  static void scatter(GlobalAggregatesMap<T,TI>& ga, GlobalIndex global, size_t i)
110  {
111  ga[i]=global;
112  }
113  };
114 
115  template<typename T, typename O, typename I>
117  {
118  };
119 
120 #if HAVE_MPI
121 
122 #endif
123 
124  } // namespace Amg
125 
126 #if HAVE_MPI
127  // forward declaration
128  template<class T1, class T2>
130 #endif
131 
132  namespace Amg
133  {
134 
135 #if HAVE_MPI
136 
145  template<typename T, typename O, typename T1, typename T2>
147  {
148  typedef T Vertex;
149  typedef O OverlapFlags;
153 
154  static void publish(AggregatesMap<Vertex>& aggregates,
155  ParallelInformation& pinfo,
156  const GlobalLookupIndexSet& globalLookup)
157  {
159  GlobalMap gmap(aggregates, globalLookup);
160  pinfo.copyOwnerToAll(gmap,gmap);
161  // communication only needed for ALU
162  // (ghosts with same global id as owners on the same process)
163  if (pinfo.getSolverCategory() ==
164  static_cast<int>(SolverCategory::nonoverlapping))
165  pinfo.copyCopyToAll(gmap,gmap);
166 
167  typedef typename ParallelInformation::RemoteIndices::const_iterator Lists;
168  Lists lists = pinfo.remoteIndices().find(pinfo.communicator().rank());
169  if(lists!=pinfo.remoteIndices().end()){
170 
171  // For periodic boundary conditions we must renumber
172  // the aggregates of vertices in the overlap whose owners are
173  // on the same process
174  Vertex maxAggregate =0;
175  typedef typename AggregatesMap<Vertex>::const_iterator Iter;
176  for(Iter i=aggregates.begin(), end=aggregates.end(); i!=end; ++i)
177  maxAggregate = std::max(maxAggregate, *i);
178 
179  // Compute new mapping of aggregates in the overlap that we also own
180  std::map<Vertex,Vertex> newMapping;
181 
182  // insert all elements into map
183  typedef typename ParallelInformation::RemoteIndices::RemoteIndexList
184  ::const_iterator RIter;
185  for(RIter ri=lists->second.first->begin(), rend = lists->second.first->end();
186  ri!=rend; ++ri)
187  if(O::contains(ri->localIndexPair().local().attribute()))
188  newMapping.insert(std::make_pair(aggregates[ri->localIndexPair().local()],
189  maxAggregate));
190  // renumber
191  typedef typename std::map<Vertex,Vertex>::iterator MIter;
192  for(MIter mi=newMapping.begin(), mend=newMapping.end();
193  mi != mend; ++mi)
194  mi->second=++maxAggregate;
195 
196 
197  for(RIter ri=lists->second.first->begin(), rend = lists->second.first->end();
198  ri!=rend; ++ri)
199  if(O::contains(ri->localIndexPair().local().attribute()))
200  aggregates[ri->localIndexPair().local()] =
201  newMapping[aggregates[ri->localIndexPair().local()]];
202  }
203  }
204  };
205 #endif
206 
207  template<typename T, typename O>
209  {
210  typedef T Vertex;
213 
214  static void publish(AggregatesMap<Vertex>& aggregates,
215  ParallelInformation& pinfo,
216  const GlobalLookupIndexSet& globalLookup)
217  {}
218  };
219 
220  } // end Amg namespace
221 
222 
223 #if HAVE_MPI
224  template<typename T, typename TI>
225  struct CommPolicy<Amg::GlobalAggregatesMap<T,TI> >
226  {
229  typedef SizeOne IndexedTypeFlag;
230  static int getSize(const Type&, int)
231  {
232  return 1;
233  }
234  };
235 #endif
236 
237 } // end Dune namespace
238 
239 #endif