Dune Core Modules (2.6.0)

globalindexset.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
33 #ifndef DUNE_GRID_UTILITY_GLOBALINDEXSET_HH
34 #define DUNE_GRID_UTILITY_GLOBALINDEXSET_HH
35 
37 #include <vector>
38 #include <iostream>
39 #include <fstream>
40 #include <memory>
41 #include <map>
42 #include <utility>
43 #include <algorithm>
44 
46 #include <dune/grid/common/gridenums.hh>
48 
50 #if HAVE_MPI
52 #endif
53 
54 namespace Dune
55 {
56 
59  template<class GridView>
61  {
62  public:
64  typedef int Index;
65 
71  template <class Entity, int Codim>
73  {
76  static PartitionType get(const Entity& entity, int codim, int i)
77  {
78  if (codim==Codim)
79  return entity.template subEntity<Codim>(i).partitionType();
80  else
81  return SubPartitionTypeProvider<Entity,Codim-1>::get(entity, codim, i);
82  }
83  };
84 
85  template <class Entity>
86  struct SubPartitionTypeProvider<Entity,0>
87  {
88  static PartitionType get(const Entity& entity, int codim, int i)
89  {
90  return entity.template subEntity<0>(i).partitionType();
91  }
92  };
93 
94  private:
96  typedef typename GridView::Grid Grid;
97 
98  typedef typename GridView::Grid::GlobalIdSet GlobalIdSet;
99  typedef typename GridView::Grid::GlobalIdSet::IdType IdType;
100  typedef typename GridView::Traits::template Codim<0>::Iterator Iterator;
101 
102  typedef typename Grid::CollectiveCommunication CollectiveCommunication;
103 
104  typedef std::map<IdType,Index> MapId2Index;
105  typedef std::map<Index,Index> IndexMap;
106 
107  /*********************************************************************************************/
108  /* calculate unique partitioning for all entities of a given codim in a given GridView, */
109  /* assuming they all have the same geometry, i.e. codim, type */
110  /*********************************************************************************************/
111  class UniqueEntityPartition
112  {
113  private:
114  /* A DataHandle class to calculate the minimum of a std::vector which is accompanied by an index set */
115  template<class IS, class V> // mapper type and vector type
116  class MinimumExchange
117  : public Dune::CommDataHandleIF<MinimumExchange<IS,V>,typename V::value_type>
118  {
119  public:
121  typedef typename V::value_type DataType;
122 
124  bool contains (int dim, unsigned int codim) const
125  {
126  return codim==indexSetCodim_;
127  }
128 
130  bool fixedSize (int dim, int codim) const
131  {
132  return true ;
133  }
134 
138  template<class EntityType>
139  size_t size (EntityType& e) const
140  {
141  return 1 ;
142  }
143 
145  template<class MessageBuffer, class EntityType>
146  void gather (MessageBuffer& buff, const EntityType& e) const
147  {
148  buff.write(v_[indexset_.index(e)]);
149  }
150 
155  template<class MessageBuffer, class EntityType>
156  void scatter (MessageBuffer& buff, const EntityType& e, size_t n)
157  {
158  DataType x;
159  buff.read(x);
160  if (x>=0) // other is -1 means, he does not want it
161  v_[indexset_.index(e)] = std::min(x,v_[indexset_.index(e)]);
162  }
163 
165  MinimumExchange (const IS& indexset, V& v, unsigned int indexSetCodim)
166  : indexset_(indexset),
167  v_(v),
168  indexSetCodim_(indexSetCodim)
169  {}
170 
171  private:
172  const IS& indexset_;
173  V& v_;
174  unsigned int indexSetCodim_;
175  };
176 
177  public:
180  UniqueEntityPartition (const GridView& gridview, unsigned int codim)
181  : assignment_(gridview.size(codim))
182  {
184  typedef typename GridView::IndexSet IndexSet;
185 
186  // assign own rank to entities that I might have
187  for (auto it = gridview.template begin<0>(); it!=gridview.template end<0>(); ++it)
188  for (unsigned int i=0; i<it->subEntities(codim); i++)
189  {
190  // Evil hack: I need to call subEntity, which needs the entity codimension as a static parameter.
191  // However, we only have it as a run-time parameter.
192  PartitionType subPartitionType = SubPartitionTypeProvider<typename GridView::template Codim<0>::Entity, GridView::dimension>::get(*it,codim,i);
193 
194  assignment_[gridview.indexSet().subIndex(*it,i,codim)]
195  = ( subPartitionType==Dune::InteriorEntity or subPartitionType==Dune::BorderEntity )
196  ? gridview.comm().rank() // set to own rank
197  : - 1; // it is a ghost entity, I will not possibly own it.
198  }
199 
201  MinimumExchange<IndexSet,std::vector<Index> > dh(gridview.indexSet(),assignment_,codim);
202 
203  gridview.communicate(dh,Dune::All_All_Interface,Dune::ForwardCommunication);
204  }
205 
207  int owner(size_t i)
208  {
209  return assignment_[i];
210  }
211 
213  size_t numOwners(int rank) const
214  {
215  return std::count(assignment_.begin(), assignment_.end(), rank);
216  }
217 
218  private:
219  std::vector<int> assignment_;
220  };
221 
222  private:
223  /* A DataHandle class to communicate the global index from the
224  * owning to the non-owning entity; the class is based on the MinimumExchange
225  * class in the parallelsolver.hh header file.
226  */
227  class IndexExchange
228  : public Dune::CommDataHandleIF<IndexExchange,Index>
229  {
230  public:
232  bool contains (int dim, unsigned int codim) const
233  {
234  return codim==indexSetCodim_;
235  }
236 
238  bool fixedSize (int dim, int codim) const
239  {
240  return true;
241  }
242 
247  template<class EntityType>
248  size_t size (EntityType& e) const
249  {
250  return 1;
251  }
252 
254  template<class MessageBuffer, class EntityType>
255  void gather (MessageBuffer& buff, const EntityType& e) const
256  {
257  IdType id=globalidset_.id(e);
258 
259  if (indexSetCodim_==0)
260  buff.write(mapid2entity_[id]);
261  else
262  buff.write((*mapid2entity_.find(id)).second);
263  }
264 
269  template<class MessageBuffer, class EntityType>
270  void scatter (MessageBuffer& buff, const EntityType& entity, size_t n)
271  {
272  Index x;
273  buff.read(x);
274 
282  if(x >= 0) {
283  const IdType id = globalidset_.id(entity);
284 
285  if (indexSetCodim_==0)
286  mapid2entity_[id] = x;
287  else
288  {
289  mapid2entity_.erase(id);
290  mapid2entity_.insert(std::make_pair(id,x));
291 
292  const Index lindex = indexSet_.index(entity);
293  localGlobalMap_[lindex] = x;
294  }
295  }
296  }
297 
299  IndexExchange (const GlobalIdSet& globalidset, MapId2Index& mapid2entity,
300  const typename GridView::IndexSet& localIndexSet, IndexMap& localGlobal,
301  unsigned int indexSetCodim)
302  : globalidset_(globalidset),
303  mapid2entity_(mapid2entity),
304  indexSet_(localIndexSet),
305  localGlobalMap_(localGlobal),
306  indexSetCodim_(indexSetCodim)
307  {}
308 
309  private:
310  const GlobalIdSet& globalidset_;
311  MapId2Index& mapid2entity_;
312 
313  const typename GridView::IndexSet& indexSet_;
314  IndexMap& localGlobalMap_;
315  unsigned int indexSetCodim_;
316  };
317 
318  public:
324  GlobalIndexSet(const GridView& gridview, int codim)
325  : gridview_(gridview),
326  codim_(codim)
327  {
328  int rank = gridview.comm().rank();
329  int size = gridview.comm().size();
330 
331  const typename GridView::IndexSet& indexSet = gridview.indexSet();
332 
333  std::unique_ptr<UniqueEntityPartition> uniqueEntityPartition;
334  if (codim_!=0)
335  uniqueEntityPartition = std::unique_ptr<UniqueEntityPartition>(new UniqueEntityPartition(gridview,codim_));
336 
337  int nLocalEntity = (codim_==0)
338  ? std::distance(gridview.template begin<0, Dune::Interior_Partition>(), gridview.template end<0, Dune::Interior_Partition>())
339  : uniqueEntityPartition->numOwners(rank);
340 
341  // Compute the global, non-redundant number of entities, i.e. the number of entities in the set
342  // without double, aka. redundant entities, on the interprocessor boundary via global reduce. */
343  nGlobalEntity_ = gridview.comm().template sum<int>(nLocalEntity);
344 
345  /* communicate the number of locally owned entities to all other processes so that the respective offset
346  * can be calculated on the respective processor; we use the Dune mpi collective communication facility
347  * for this; first, we gather the number of locally owned entities on the root process and, second, we
348  * broadcast the array to all processes where the respective offset can be calculated. */
349 
350  std::vector<int> offset(size);
351  std::fill(offset.begin(), offset.end(), 0);
352 
354  gridview_.comm().template allgather<int>(&nLocalEntity, 1, offset.data());
355 
356  int myoffset = 0;
357  for (int i=1; i<rank+1; i++)
358  myoffset += offset[i-1];
359 
360  /* compute globally unique index over all processes; the idea of the algorithm is as follows: if
361  * an entity is owned by the process, it is assigned an index that is the addition of the offset
362  * specific for this process and a consecutively incremented counter; if the entity is not owned
363  * by the process, it is assigned -1, which signals that this specific entity will get its global
364  * unique index through communication afterwards;
365  *
366  * thus, the calculation of the globally unique index is divided into 2 stages:
367  *
368  * (1) we calculate the global index independently;
369  *
370  * (2) we achieve parallel adjustment by communicating the index
371  * from the owning entity to the non-owning entity.
372  *
373  */
374 
375  // 1st stage of global index calculation: calculate global index for owned entities
376  // initialize map that stores an entity's global index via it's globally unique id as key
377  globalIndex_.clear();
378 
379  const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet();
381  Index globalcontrib = 0;
383  if (codim_==0) // This case is simpler
384  {
385  for (Iterator iter = gridview_.template begin<0>(); iter!=gridview_.template end<0>(); ++iter)
386  {
387  const IdType id = globalIdSet.id(*iter);
390  if (iter->partitionType() == Dune::InteriorEntity)
391  {
392  const Index gindex = myoffset + globalcontrib;
394  globalIndex_[id] = gindex;
395  globalcontrib++;
396  }
397 
399  else
400  {
401  globalIndex_[id] = -1;
402  }
403  }
404  }
405  else // if (codim==0) else
406  {
407  std::vector<bool> firstTime(gridview_.size(codim_));
408  std::fill(firstTime.begin(), firstTime.end(), true);
409 
410  for(Iterator iter = gridview_.template begin<0>();iter!=gridview_.template end<0>(); ++iter)
411  {
412  for (size_t i=0; i<iter->subEntities(codim_); i++)
413  {
414  IdType id=globalIdSet.subId(*iter,i,codim_);
415 
416  Index idx = gridview_.indexSet().subIndex(*iter,i,codim_);
417 
418  if (!firstTime[idx] )
419  continue;
420 
421  firstTime[idx] = false;
422 
423  if (uniqueEntityPartition->owner(idx) == rank)
424  {
425  const Index gindex = myoffset + globalcontrib;
426  globalIndex_.insert(std::make_pair(id,gindex));
428  const Index lindex = idx;
429  localGlobalMap_[lindex] = gindex;
430 
431  globalcontrib++;
432  }
433  else
434  {
435  globalIndex_.insert(std::make_pair(id,-1));
436  }
437  }
438 
439  }
440  }
441 
442  // 2nd stage of global index calculation: communicate global index for non-owned entities
443 
444  // Create the data handle and communicate.
445  IndexExchange dataHandle(globalIdSet,globalIndex_,indexSet,localGlobalMap_,codim_);
447  }
448 
450  template <class Entity>
451  Index index(const Entity& entity) const
452  {
453  if (codim_==0)
454  {
456  const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet();
457  const IdType id = globalIdSet.id(entity);
458  const Index gindex = globalIndex_.find(id)->second;
460  return gindex;
461  }
462  else
463  return localGlobalMap_.find(gridview_.indexSet().index(entity))->second;
464  }
465 
471  template <class Entity>
472  Index subIndex(const Entity& entity, unsigned int i, unsigned int codim) const
473  {
474  if (codim_==0)
475  {
477  const GlobalIdSet& globalIdSet = gridview_.grid().globalIdSet();
478  const IdType id = globalIdSet.subId(entity,i,codim);
479  const Index gindex = globalIndex_.find(id)->second;
481  return gindex;
482  }
483  else
484  return localGlobalMap_.find(gridview_.indexSet().subIndex(entity,i,codim))->second;
485  }
486 
492  unsigned int size(unsigned int codim) const
493  {
494  return (codim_==codim) ? nGlobalEntity_ : 0;
495  }
496 
497  protected:
498  const GridView gridview_;
499 
501  unsigned int codim_;
502 
505 
506  IndexMap localGlobalMap_;
507 
510  MapId2Index globalIndex_;
511  };
512 
513 } // namespace Dune
514 
515 #endif /* DUNE_GRID_UTILITY_GLOBALINDEXSET_HH */
CommDataHandleIF describes the features of a data handle for communication in parallel runs using the...
Definition: datahandleif.hh:76
void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
unpack data from message buffer to user.
Definition: datahandleif.hh:167
bool contains(int dim, int codim) const
returns true if data for given valid codim should be communicated
Definition: datahandleif.hh:91
void gather(MessageBufferImp &buff, const EntityType &e) const
pack data from user to message buffer
Definition: datahandleif.hh:153
Wrapper class for entities.
Definition: entity.hh:64
PartitionType partitionType() const
Partition type of this entity.
Definition: entity.hh:129
Calculate globally unique index over all processes in a Dune grid.
Definition: globalindexset.hh:61
Index subIndex(const Entity &entity, unsigned int i, unsigned int codim) const
Return the global index of a subentity of a given entity.
Definition: globalindexset.hh:472
MapId2Index globalIndex_
Stores global index of entities with entity's globally unique id as key.
Definition: globalindexset.hh:510
int Index
The number type used for global indices
Definition: globalindexset.hh:64
int nGlobalEntity_
Global number of entities, i.e. number of entities without rendundant entities on interprocessor boun...
Definition: globalindexset.hh:504
Index index(const Entity &entity) const
Return the global index of a given entity.
Definition: globalindexset.hh:451
unsigned int codim_
Codimension of the entities that we hold indices for.
Definition: globalindexset.hh:501
unsigned int size(unsigned int codim) const
Return the total number of entities over all processes that we have indices for.
Definition: globalindexset.hh:492
GlobalIndexSet(const GridView &gridview, int codim)
Constructor for a given GridView.
Definition: globalindexset.hh:324
Grid view abstract base class.
Definition: gridview.hh:60
GridFamily::Traits::CollectiveCommunication CollectiveCommunication
A type that is a model of Dune::CollectiveCommunication. It provides a portable way for collective co...
Definition: grid.hh:519
Describes the parallel communication interface class for MessageBuffers and DataHandles.
int size(int codim) const
obtain number of entities in a given codimension
Definition: gridview.hh:178
Traits ::IndexSet IndexSet
type of the index set
Definition: gridview.hh:80
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir) const
Communicate data on this view.
Definition: gridview.hh:265
const Grid & grid() const
obtain a const reference to the underlying hierarchic grid
Definition: gridview.hh:161
const CollectiveCommunication & comm() const
obtain collective communication object
Definition: gridview.hh:246
Traits ::Grid Grid
type of the grid
Definition: gridview.hh:77
const IndexSet & indexSet() const
obtain the index set
Definition: gridview.hh:172
@ dimension
The dimension of the grid.
Definition: gridview.hh:127
PartitionType
Attributes used in the generic overlap model.
Definition: gridenums.hh:28
@ InteriorEntity
all interior entities
Definition: gridenums.hh:29
@ BorderEntity
on boundary between interior and overlap
Definition: gridenums.hh:30
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:169
@ All_All_Interface
send all and receive all entities
Definition: gridenums.hh:89
Helpers for dealing with MPI.
Dune namespace.
Definition: alignedallocator.hh:10
Static tag representing a codimension.
Definition: dimension.hh:22
Helper class to provide access to subentity PartitionTypes with a run-time codimension.
Definition: globalindexset.hh:73
static PartitionType get(const Entity &entity, int codim, int i)
Get PartitionType of the i-th subentity of codimension 'codim' of entity 'entity'.
Definition: globalindexset.hh:76
std::size_t fixedSize
The number of data items per index if it is fixed, 0 otherwise.
Definition: variablesizecommunicator.hh:245
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.80.0 (Feb 21, 23:30, 2024)