1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
|
#ifndef DUNE_UGGRID_INDEXSETS_HH
#define DUNE_UGGRID_INDEXSETS_HH
/** \file
\brief The index and id sets for the UGGrid class
*/
#include <vector>
#include <set>
#include <dune/grid/common/grid.hh>
namespace Dune {
template<class GridImp>
class UGGridLevelIndexSet : public IndexSet<GridImp,UGGridLevelIndexSet<GridImp>, UG::UINT>
{
enum {dim = GridImp::dimension};
public:
/** \brief Default constructor
Unfortunately we can't force the user to init grid_ and level_, because
UGGridLevelIndexSets are meant to be stored in an array.
\todo I want to make this constructor private, but I can't, because
it is called by UGGrid through a std::vector::resize()
*/
UGGridLevelIndexSet () {}
//! get index of an entity
template<int cd>
unsigned int index (const typename GridImp::Traits::template Codim<cd>::Entity& e) const
{
return UG_NS<dim>::levelIndex(grid_->getRealImplementation(e).target_);
}
//! get index of subEntity of a codim 0 entity
template<int cc>
unsigned int subIndex (const typename GridImp::Traits::template Codim<cc>::Entity& e,
int i,
unsigned int codim) const
{
if (cc==dim)
return UG_NS<dim>::levelIndex(grid_->getRealImplementation(e).target_);
if (codim==dim)
return UG_NS<dim>::levelIndex(UG_NS<dim>::Corner(grid_->getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(i,e.type())));
if (codim==0)
return UG_NS<dim>::levelIndex(grid_->getRealImplementation(e).target_);
if (codim==dim-1) {
int a=GenericReferenceElements<double,dim>::general(e.type()).subEntity(i,dim-1,0,dim);
int b=GenericReferenceElements<double,dim>::general(e.type()).subEntity(i,dim-1,1,dim);
return UG_NS<dim>::levelIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(grid_->getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(a,e.type())),
UG_NS<dim>::Corner(grid_->getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(b,e.type()))));
}
if (codim==1)
return UG_NS<dim>::levelIndex(UG_NS<dim>::SideVector(grid_->getRealImplementation(e).target_,
UGGridRenumberer<dim>::facesDUNEtoUG(i,e.type())));
DUNE_THROW(GridError, "UGGrid<" << dim << "," << dim << ">::subIndex isn't implemented for codim==" << codim );
}
//! get number of entities of given codim, type and on this level
int size (int codim) const {
if (codim==0)
return numSimplices_+numPyramids_+numPrisms_+numCubes_;
if (codim==dim)
return numVertices_;
if (codim==dim-1)
return numEdges_;
if (codim==1)
return numTriFaces_+numQuadFaces_;
DUNE_THROW(NotImplemented, "wrong codim!");
}
//! get number of entities of given codim, type and on this level
int size (GeometryType type) const
{
int codim = GridImp::dimension-type.dim();
if (codim==0) {
if (type.isSimplex())
return numSimplices_;
else if (type.isPyramid())
return numPyramids_;
else if (type.isPrism())
return numPrisms_;
else if (type.isCube())
return numCubes_;
else
return 0;
}
if (codim==dim) {
return numVertices_;
}
if (codim==dim-1) {
return numEdges_;
}
if (codim==1) {
if (type.isSimplex())
return numTriFaces_;
else if (type.isCube())
return numQuadFaces_;
else
return 0;
}
DUNE_THROW(NotImplemented, "Wrong codim!");
}
/** \brief Deliver all geometry types used in this grid */
const std::vector<GeometryType>& geomTypes (int codim) const
{
return myTypes_[codim];
}
/** \brief Return true if e is contained in the index set.
\note If 'entity' is from another grid this method may still return 'true'.
This is acceptable by the Dune grid interface specification.
*/
template <class EntityType>
bool contains (const EntityType& entity) const
{
return entity.level() == level_;
}
/** \brief Update the level indices. This method is called after each grid change */
void update(const GridImp& grid, int level, std::vector<unsigned int>* nodePermutation=0);
const GridImp* grid_;
int level_;
int numSimplices_;
int numPyramids_;
int numPrisms_;
int numCubes_;
int numVertices_;
int numEdges_;
int numTriFaces_;
int numQuadFaces_;
std::vector<GeometryType> myTypes_[dim+1];
};
template<class GridImp>
class UGGridLeafIndexSet: public IndexSet<GridImp,UGGridLeafIndexSet<GridImp>, UG::UINT>
{
public:
/*
We use the remove_const to extract the Type from the mutable class,
because the const class is not instantiated yet.
*/
enum {dim = remove_const<GridImp>::type::dimension};
//! constructor stores reference to a grid and level
UGGridLeafIndexSet (const GridImp& g)
: grid_(g), coarsestLevelWithLeafElements_(0)
{
}
//! get index of an entity
/*
We use the remove_const to extract the Type from the mutable class,
because the const class is not instantiated yet.
*/
template<int cd>
int index (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
{
return UG_NS<dim>::leafIndex(grid_.getRealImplementation(e).target_);
}
//! get index of subEntity of a codim 0 entity
/*
We use the remove_const to extract the Type from the mutable class,
because the const class is not instantiated yet.
*/
template<int cc>
unsigned int subIndex (const typename remove_const<GridImp>::type::Traits::template Codim<cc>::Entity& e,
int i,
unsigned int codim) const
{
if (cc==dim)
return UG_NS<dim>::levelIndex(grid_.getRealImplementation(e).target_);
if (codim==0)
return UG_NS<dim>::leafIndex(grid_.getRealImplementation(e).target_);
const GeometryType type = e.type();
if (codim==dim)
return UG_NS<dim>::leafIndex(UG_NS<dim>::Corner(grid_.getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(i,type)));
if (codim==dim-1) {
int a=GenericReferenceElements<double,dim>::general(type).subEntity(i,dim-1,0,dim);
int b=GenericReferenceElements<double,dim>::general(type).subEntity(i,dim-1,1,dim);
return UG_NS<dim>::leafIndex(UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(grid_.getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(a,type)),
UG_NS<dim>::Corner(grid_.getRealImplementation(e).target_,
UGGridRenumberer<dim>::verticesDUNEtoUG(b,type))));
}
if (codim==1)
return UG_NS<dim>::leafIndex(UG_NS<dim>::SideVector(grid_.getRealImplementation(e).target_,
UGGridRenumberer<dim>::facesDUNEtoUG(i,type)));
DUNE_THROW(GridError, "UGGrid<" << dim << "," << dim << ">::subLeafIndex isn't implemented for codim==" << codim );
}
//! get number of entities of given codim and type
int size (GeometryType type) const
{
if (type.dim()==GridImp::dimension) {
if (type.isSimplex())
return numSimplices_;
else if (type.isPyramid())
return numPyramids_;
else if (type.isPrism())
return numPrisms_;
else if (type.isCube())
return numCubes_;
else
return 0;
}
if (type.dim()==0) {
return numVertices_;
}
if (type.dim()==1) {
return numEdges_;
}
if (type.isTriangle())
return numTriFaces_;
else if (type.isQuadrilateral())
return numQuadFaces_;
return 0;
}
//! get number of entities of given codim
int size (int codim) const
{
int s=0;
const std::vector<GeometryType>& geomTs = geomTypes(codim);
for (unsigned int i=0; i<geomTs.size(); ++i)
s += size(geomTs[i]);
return s;
}
/** deliver all geometry types used in this grid */
const std::vector<GeometryType>& geomTypes (int codim) const
{
return myTypes_[codim];
}
/** \brief Return true if e is contained in the index set.
\note If 'entity' is from another grid this method may still return 'true'.
This is acceptable by the Dune grid interface specification.
*/
template <class EntityType>
bool contains (const EntityType& entity) const
{
return UG_NS<dim>::isLeaf(GridImp::getRealImplementation(entity).target_);
}
/** \brief Update the leaf indices. This method is called after each grid change. */
void update(std::vector<unsigned int>* nodePermutation=0);
const GridImp& grid_;
/** \brief The lowest level that contains leaf elements
This corresponds to UG's fullRefineLevel, which is, unfortunately only
computed if you use some nontrivial UG algebra. Thus we compute it
ourselves, and use it to speed up the leaf iterators.
*/
unsigned int coarsestLevelWithLeafElements_;
int numSimplices_;
int numPyramids_;
int numPrisms_;
int numCubes_;
int numVertices_;
int numEdges_;
int numTriFaces_;
int numQuadFaces_;
std::vector<GeometryType> myTypes_[dim+1];
};
/** \brief Implementation class for the UGGrid Id sets
The UGGridGlobalIdSet and the UGGridLocalIdSet are virtually identical. This
class implements them both at once. You can select the one you want using
the <tt>Local</tt> template parameter.
\tparam Local false for GlobalIdSet, true for LocalIdSet
*/
template <class GridImp, bool Local>
class UGGridIdSet : public IdSet<GridImp,UGGridIdSet<GridImp,Local>,unsigned int>
{
enum {dim = remove_const<GridImp>::type::dimension};
typedef typename std::pair<const typename UG_NS<dim>::Element*, int> Face;
/** \brief Look for copy of a face on the next-lower grid level.
\todo This method is not implemented very efficiently, but I will not put
further work into it as long as the much-awaited face objects are not included
in UG.
*/
static Face getFatherFace(const Face& face) {
// set up result object
Face resultFace;
resultFace.first = UG_NS<dim>::EFather(face.first);
// If there is no father element then we know there is no father face
/** \bug This is not true when doing vertical load balancing. */
if (resultFace.first == NULL)
return resultFace;
// Get all corners of the face
std::set<const typename UG_NS<dim>::Vertex*> corners;
for (int i=0; i<UG_NS<dim>::Corners_Of_Side(face.first, face.second); i++)
corners.insert(UG_NS<dim>::Corner(face.first, UG_NS<dim>::Corner_Of_Side(face.first, face.second, i))->myvertex);
// Loop over all faces of the father element and look for a face that has the same vertices
for (int i=0; i<UG_NS<dim>::Sides_Of_Elem(resultFace.first); i++) {
// Copy father face into set
std::set<const typename UG_NS<dim>::Vertex*> fatherFaceCorners;
for (int j=0; j<UG_NS<dim>::Corners_Of_Side(resultFace.first, i); j++)
fatherFaceCorners.insert(UG_NS<dim>::Corner(resultFace.first, UG_NS<dim>::Corner_Of_Side(resultFace.first, i, j))->myvertex);
// Do the vertex sets match?
if (corners==fatherFaceCorners) {
resultFace.second = i;
return resultFace;
}
}
// No father face found
resultFace.first = NULL;
return resultFace;
}
public:
//! constructor stores reference to a grid
UGGridIdSet (const GridImp& g) : grid_(g) {}
//! get id of an entity
/*
We use the remove_const to extract the Type from the mutable class,
because the const class is not instantiated yet.
*/
template<int cd>
unsigned int id (const typename remove_const<GridImp>::type::Traits::template Codim<cd>::Entity& e) const
{
if (cd==0) {
// If we're asked for the id of an element, and that element is a copy of its father, then
// we return the id of the lowest ancestor that the element is a copy from. That way copies
// of elements have the same id
const typename UG_NS<dim>::Element* ancestor = (typename UG_NS<dim>::Element* const)(grid_.getRealImplementation(e).target_);
/** \todo We should really be using an isCopy() method rather than hasCopy() */
while (UG_NS<dim>::EFather(ancestor) && UG_NS<dim>::hasCopy(UG_NS<dim>::EFather(ancestor)))
ancestor = UG_NS<dim>::EFather(ancestor);
#if defined ModelP
return ancestor->ge.ddd.gid;
#else
return UG_NS<dim>::id(ancestor);
#endif
}
#if defined ModelP
if (cd == dim) {
typename UG_NS<dim>::Node *node =
reinterpret_cast<typename UG_NS<dim>::Node *>(grid_.getRealImplementation(e).target_);
return node->ddd.gid;
}
else {
DUNE_THROW(NotImplemented,
(Local ? "Local" : "Global") <<
" persistent index for entities which are neither nodes nor elements.");
}
#else
return UG_NS<dim>::id(grid_.getRealImplementation(e).target_);
#endif
}
//! get id of subEntity
/*
We use the remove_const to extract the Type from the mutable class,
because the const class is not instantiated yet.
*/
unsigned int subId (const typename remove_const<GridImp>::type::Traits::template Codim<0>::Entity& e,
int i,
unsigned int codim) const
{
if (codim==0)
return id<0>(e);
const typename UG_NS<dim>::Element* target = grid_.getRealImplementation(e).target_;
GeometryType type = e.type();
if (dim-codim==1) {
int a=GenericReferenceElements<double,dim>::general(type).subEntity(i,dim-1,0,dim);
int b=GenericReferenceElements<double,dim>::general(type).subEntity(i,dim-1,1,dim);
const typename UG_NS<dim>::Edge* edge = UG_NS<dim>::GetEdge(UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(a,type)),
UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(b,type)));
// If this edge is the copy of an edge on a lower level we return the id of that lower
// edge, because Dune wants entities which are copies of each other to have the same id.
const typename UG_NS<dim>::Edge* fatherEdge;
fatherEdge = GetFatherEdge(edge);
while (fatherEdge // fatherEdge exists
// ... and it must be a true copy father
&& ( (fatherEdge->links[0].nbnode->myvertex == edge->links[0].nbnode->myvertex
&& fatherEdge->links[1].nbnode->myvertex == edge->links[1].nbnode->myvertex)
||
(fatherEdge->links[0].nbnode->myvertex == edge->links[1].nbnode->myvertex
&& fatherEdge->links[1].nbnode->myvertex == edge->links[0].nbnode->myvertex) ) ) {
edge = fatherEdge;
fatherEdge = GetFatherEdge(edge);
}
#ifdef ModelP
//return (Local) ? edge->id : edge->ddd.gid;
DUNE_THROW(NotImplemented, "!");
#else
return edge->id;
#endif
}
if (codim==1) { // Faces
Face face(target, UGGridRenumberer<dim>::facesDUNEtoUG(i,type));
// If this face is the copy of a face on a lower level we return the id of that lower
// face, because Dune wants entities which are copies of each other to have the same id.
Face fatherFace;
fatherFace = getFatherFace(face);
while (fatherFace.first) {
face = fatherFace;
fatherFace = getFatherFace(face);
}
#ifdef ModelP
return (Local)
? UG_NS<dim>::SideVector(face.first, face.second)->id
: UG_NS<dim>::SideVector(face.first, face.second)->ddd.gid;
#else
return UG_NS<dim>::SideVector(face.first, face.second)->id;
#endif
}
if (codim==dim) {
#ifdef ModelP
return (Local)
? UG_NS<dim>::id(UG_NS<dim>::Corner(target,UGGridRenumberer<dim>::verticesDUNEtoUG(i,type)))
: UG_NS<dim>::Corner(target, UGGridRenumberer<dim>::verticesDUNEtoUG(i,type))->ddd.gid;
#else
return UG_NS<dim>::id(UG_NS<dim>::Corner(target,UGGridRenumberer<dim>::verticesDUNEtoUG(i,type)));
#endif
}
DUNE_THROW(GridError, "UGGrid<" << dim << ">::subId isn't implemented for codim==" << codim );
}
//private:
const GridImp& grid_;
};
} // namespace Dune
#endif
|