← Back to branch summary

~ubuntu-branches/ubuntu/utopic/dune-grid/utopic-proposed

« back to all changes in this revision

Viewing changes to dune/grid/uggrid/uggridindexsets.hh

  • Committer: Package Import Robot
  • Author(s): Ansgar Burchardt
  • Date: 2012-04-06 11:31:20 UTC
  • Revision ID: package-import@ubuntu.com-20120406113120-x0z4e0qqgfhmaj2a
Tags: upstream-2.2~svn7982
Import upstream version 2.2~svn7982

Show diffs side-by-side

added added

removed removed

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