Dune Core Modules (2.4.1)

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
54namespace 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
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:73
void scatter(MessageBufferImp &buff, const EntityType &e, size_t n)
unpack data from message buffer to user.
Definition: datahandleif.hh:133
bool contains(int dim, int codim) const
returns true if data for given valid codim should be communicated
Definition: datahandleif.hh:88
void gather(MessageBufferImp &buff, const EntityType &e) const
pack data from user to message buffer
Definition: datahandleif.hh:121
Wrapper class for entities.
Definition: entity.hh:62
PartitionType partitionType() const
Partition type of this entity.
Definition: entity.hh:128
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:59
GridFamily::Traits::CollectiveCommunication CollectiveCommunication
A type that is a model of Dune::CollectiveCommunication. It provides a portable way for collective co...
Definition: grid.hh:545
Describes the parallel communication interface class for MessageBuffers and DataHandles.
const IndexSet & indexSet() const
obtain the index set
Definition: gridview.hh:175
Traits::Grid Grid
type of the grid
Definition: gridview.hh:77
const Grid & grid() const
obtain a const reference to the underlying hierarchic grid
Definition: gridview.hh:164
const CollectiveCommunication & comm() const
obtain collective communication object
Definition: gridview.hh:249
Traits::IndexSet IndexSet
type of the index set
Definition: gridview.hh:80
int size(int codim) const
obtain number of entities in a given codimension
Definition: gridview.hh:181
void communicate(CommDataHandleIF< DataHandleImp, DataType > &data, InterfaceType iftype, CommunicationDirection dir) const
Communicate data on this view.
Definition: gridview.hh:268
@ dimension
The dimension of the grid.
Definition: gridview.hh:130
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: alignment.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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)