DUNE-FEM (unstable)

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