dune-grid  2.2.1
misc.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALBERTA_MISC_HH
2 #define DUNE_ALBERTA_MISC_HH
3 
4 #include <cassert>
5 
6 #include <dune/common/exceptions.hh>
7 #include <dune/common/typetraits.hh>
8 #include <dune/common/forloop.hh>
9 
11 
12 #if HAVE_ALBERTA
13 
14 // should the coordinates be cached in a vector (required for ALBERTA 2.0)?
15 #ifndef DUNE_ALBERTA_CACHE_COORDINATES
16 #define DUNE_ALBERTA_CACHE_COORDINATES 1
17 #endif
18 
19 // set to 1 to use generic geometries in AlbertaGrid
20 #ifndef DUNE_ALBERTA_USE_GENERICGEOMETRY
21 #define DUNE_ALBERTA_USE_GENERICGEOMETRY 0
22 #endif
23 
24 namespace Dune
25 {
26 
27  // Exceptions
28  // ----------
29 
31  : public Exception
32  {};
33 
35  : public IOError
36  {};
37 
38 
39 
40  namespace Alberta
41  {
42 
43  // Import Types
44  // ------------
45 
46  static const int dimWorld = DIM_OF_WORLD;
47 
48  typedef ALBERTA REAL Real;
49  typedef ALBERTA REAL_B LocalVector; // in barycentric coordinates
50  typedef ALBERTA REAL_D GlobalVector;
51  typedef ALBERTA REAL_DD GlobalMatrix;
52 
53 #if DUNE_ALBERTA_VERSION >= 0x300
54  typedef ALBERTA AFF_TRAFO AffineTransformation;
55 #else
57  {
60  };
61 #endif
62 
63  typedef ALBERTA MESH Mesh;
64  typedef ALBERTA EL Element;
65 
66  static const int meshRefined = MESH_REFINED;
67  static const int meshCoarsened = MESH_COARSENED;
68 
69  static const int InteriorBoundary = INTERIOR;
70  static const int DirichletBoundary = DIRICHLET;
71 #if DUNE_ALBERTA_VERSION >= 0x300
72  typedef ALBERTA BNDRY_TYPE BoundaryId;
73 #else
74  typedef S_CHAR BoundaryId;
75 #endif
76 
77  typedef U_CHAR ElementType;
78 
79  typedef ALBERTA FE_SPACE DofSpace;
80 
81 
82 
83  // Memory Manipulation Functions
84  // -----------------------------
85 
86  template< class Data >
87  inline Data *memAlloc ( size_t size )
88  {
89  return MEM_ALLOC( size, Data );
90  }
91 
92  template< class Data >
93  inline Data *memCAlloc ( size_t size )
94  {
95  return MEM_CALLOC( size, Data );
96  }
97 
98  template< class Data >
99  inline Data *memReAlloc ( Data *ptr, size_t oldSize, size_t newSize )
100  {
101  return MEM_REALLOC( ptr, oldSize, newSize, Data );
102  }
103 
104  template< class Data >
105  inline void memFree ( Data *ptr, size_t size )
106  {
107  return MEM_FREE( ptr, size, Data );
108  }
109 
110 
111 
112  // GlobalSpace
113  // -----------
114 
116  {
117  typedef GlobalSpace This;
118 
119  public:
122 
123  private:
124  Matrix identityMatrix_;
125  Vector nullVector_;
126 
127  GlobalSpace ()
128  {
129  for( int i = 0; i < dimWorld; ++i )
130  {
131  for( int j = 0; j < dimWorld; ++j )
132  identityMatrix_[ i ][ j ] = Real( 0 );
133  identityMatrix_[ i ][ i ] = Real( 1 );
134  nullVector_[ i ] = Real( 0 );
135  }
136  }
137 
138  static This &instance ()
139  {
140  static This theInstance;
141  return theInstance;
142  }
143 
144  public:
145  static const Matrix &identityMatrix ()
146  {
147  return instance().identityMatrix_;
148  }
149 
150  static const Vector &nullVector ()
151  {
152  return instance().nullVector_;
153  }
154  };
155 
156 
157 
158  // NumSubEntities
159  // --------------
160 
161  template< int dim, int codim >
163 
164  template< int dim >
165  struct NumSubEntities< dim, 0 >
166  {
167  static const int value = 1;
168  };
169 
170  template< int dim >
171  struct NumSubEntities< dim, dim >
172  {
173  static const int value = dim+1;
174  };
175 
176  template<>
177  struct NumSubEntities< 0, 0 >
178  {
179  static const int value = 1;
180  };
181 
182  template<>
183  struct NumSubEntities< 2, 1 >
184  {
185  static const int value = 3;
186  };
187 
188  template<>
189  struct NumSubEntities< 3, 1 >
190  {
191  static const int value = 4;
192  };
193 
194  template<>
195  struct NumSubEntities< 3, 2 >
196  {
197  static const int value = 6;
198  };
199 
200 
201 
202  // CodimType
203  // ---------
204 
205  template< int dim, int codim >
206  struct CodimType;
207 
208  template< int dim >
209  struct CodimType< dim, 0 >
210  {
211  static const int value = CENTER;
212  };
213 
214  template< int dim >
215  struct CodimType< dim, dim >
216  {
217  static const int value = VERTEX;
218  };
219 
220  template<>
221  struct CodimType< 2, 1 >
222  {
223  static const int value = EDGE;
224  };
225 
226  template<>
227  struct CodimType< 3, 1 >
228  {
229  static const int value = FACE;
230  };
231 
232  template<>
233  struct CodimType< 3, 2 >
234  {
235  static const int value = EDGE;
236  };
237 
238 
239 
240  // FillFlags
241  // ---------
242 
243  template< int dim >
244  struct FillFlags
245  {
246  typedef ALBERTA FLAGS Flags;
247 
248  static const Flags nothing = FILL_NOTHING;
249 
250  static const Flags coords = FILL_COORDS;
251 
252  static const Flags neighbor = FILL_NEIGH;
253 
254  static const Flags orientation = (dim == 3 ? FILL_ORIENTATION : FILL_NOTHING);
255 
256  static const Flags projection = FILL_PROJECTION;
257 
258 #if DUNE_ALBERTA_VERSION >= 0x300
259  static const Flags elementType = FILL_NOTHING;
260 #else
261  static const Flags elementType = (dim == 3 ? FILL_EL_TYPE : FILL_NOTHING);
262 #endif
263 
264 #if DUNE_ALBERTA_VERSION >= 0x300
265  static const Flags boundaryId = FILL_MACRO_WALLS;
266 #else
267  static const Flags boundaryId = FILL_BOUND;
268 #endif
269 
270 #if DUNE_ALBERTA_VERSION >= 0x300
271  static const Flags nonPeriodic = FILL_NON_PERIODIC;
272 #else
273  static const Flags nonPeriodic = FILL_NOTHING;
274 #endif
275 
278 
279 #if DUNE_ALBERTA_VERSION >= 0x300
280  static const Flags standardWithCoords = all & ~nonPeriodic & ~projection;
281 #else
282  static const Flags standardWithCoords = all;
283 #endif
284 
285 #if DUNE_ALBERTA_CACHE_COORDINATES
287 #else
288  static const Flags standard = standardWithCoords;
289 #endif
290  };
291 
292 
293 
294  // RefinementEdge
295  // --------------
296 
297  template< int dim >
299  {
300  static const int value = 0;
301  };
302 
303  template<>
304  struct RefinementEdge< 2 >
305  {
306  static const int value = 2;
307  };
308 
309 
310 
311  // Dune2AlbertaNumbering
312  // ---------------------
313 
314  template< int dim, int codim >
316  {
317  static int apply ( const int i )
318  {
319  assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
320  return i;
321  }
322  };
323 
324  template<>
325  struct Dune2AlbertaNumbering< 3, 2 >
326  {
327  static const int numSubEntities = NumSubEntities< 3, 2 >::value;
328 
329  static int apply ( const int i )
330  {
331  assert( (i >= 0) && (i < numSubEntities) );
332  static int dune2alberta[ numSubEntities ] = { 0, 3, 1, 2, 4, 5 };
333  return dune2alberta[ i ];
334  }
335  };
336 
337 
338 
339  // Generic2AlbertaNumbering
340  // ------------------------
341 
342  template< int dim, int codim >
344  {
345  static int apply ( const int i )
346  {
347  assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
348  return i;
349  }
350  };
351 
352  template< int dim >
353  struct Generic2AlbertaNumbering< dim, 1 >
354  {
355  static int apply ( const int i )
356  {
357  assert( (i >= 0) && (i < NumSubEntities< dim, 1 >::value) );
358  return dim - i;
359  }
360  };
361 
362  template<>
364  {
365  static int apply ( const int i )
366  {
367  assert( (i >= 0) && (i < NumSubEntities< 1, 1 >::value) );
368  return i;
369  }
370  };
371 
372  template<>
374  {
375  static const int numSubEntities = NumSubEntities< 3, 2 >::value;
376 
377  static int apply ( const int i )
378  {
379  assert( (i >= 0) && (i < numSubEntities) );
380  static int generic2alberta[ numSubEntities ] = { 0, 1, 3, 2, 4, 5 };
381  return generic2alberta[ i ];
382  }
383  };
384 
385 
386 
387  // NumberingMap
388  // ------------
389 
390  template< int dim, template< int, int > class Numbering = Generic2AlbertaNumbering >
392  {
394 
395  template< int codim >
396  struct Initialize;
397 
398  int *dune2alberta_[ dim+1 ];
399  int *alberta2dune_[ dim+1 ];
400  int numSubEntities_[ dim+1 ];
401 
402  NumberingMap ( const This & );
403  This &operator= ( const This & );
404 
405  public:
407  {
408  ForLoop< Initialize, 0, dim >::apply( *this );
409  }
410 
412  {
413  for( int codim = 0; codim <= dim; ++codim )
414  {
415  delete[]( dune2alberta_[ codim ] );
416  delete[]( alberta2dune_[ codim ] );
417  }
418  }
419 
420  int dune2alberta ( int codim, int i ) const
421  {
422  assert( (codim >= 0) && (codim <= dim) );
423  assert( (i >= 0) && (i < numSubEntities( codim )) );
424  return dune2alberta_[ codim ][ i ];
425  }
426 
427  int alberta2dune ( int codim, int i ) const
428  {
429  assert( (codim >= 0) && (codim <= dim) );
430  assert( (i >= 0) && (i < numSubEntities( codim )) );
431  return alberta2dune_[ codim ][ i ];
432  }
433 
434  int numSubEntities ( int codim ) const
435  {
436  assert( (codim >= 0) && (codim <= dim) );
437  return numSubEntities_[ codim ];
438  }
439  };
440 
441 
442 
443  // NumberingMap::Initialize
444  // ------------------------
445 
446  template< int dim, template< int, int > class Numbering >
447  template< int codim >
448  struct NumberingMap< dim, Numbering >::Initialize
449  {
450  static const int numSubEntities = NumSubEntities< dim, codim >::value;
451 
452  static void apply ( NumberingMap< dim, Numbering > &map )
453  {
454  map.numSubEntities_[ codim ] = numSubEntities;
455  map.dune2alberta_[ codim ] = new int[ numSubEntities ];
456  map.alberta2dune_[ codim ] = new int[ numSubEntities ];
457 
458  for( int i = 0; i < numSubEntities; ++i )
459  {
460  const int j = Numbering< dim, codim >::apply( i );
461  map.dune2alberta_[ codim ][ i ] = j;
462  map.alberta2dune_[ codim ][ j ] = i;
463  }
464  }
465  };
466 
467 
468 
469  // MapVertices
470  // -----------
471 
472  template< int dim, int codim >
473  struct MapVertices;
474 
475  template< int dim >
476  struct MapVertices< dim, 0 >
477  {
478  static int apply ( int subEntity, int vertex )
479  {
480  assert( subEntity == 0 );
481  assert( (vertex >= 0) && (vertex <= NumSubEntities< dim, dim >::value) );
482  return vertex;
483  }
484  };
485 
486  template<>
487  struct MapVertices< 2, 1 >
488  {
489  static int apply ( int subEntity, int vertex )
490  {
491  assert( (subEntity >= 0) && (subEntity < 3) );
492  assert( (vertex >= 0) && (vertex < 2) );
493  //static const int map[ 3 ][ 2 ] = { {1,2}, {2,0}, {0,1} };
494  static const int map[ 3 ][ 2 ] = { {1,2}, {0,2}, {0,1} };
495  return map[ subEntity ][ vertex ];
496  }
497  };
498 
499  template<>
500  struct MapVertices< 3, 1 >
501  {
502  static int apply ( int subEntity, int vertex )
503  {
504  assert( (subEntity >= 0) && (subEntity < 4) );
505  assert( (vertex >= 0) && (vertex < 3) );
506  //static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,3,2}, {0,1,3}, {0,2,1} };
507  static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,2,3}, {0,1,3}, {0,1,2} };
508  return map[ subEntity ][ vertex ];
509  }
510  };
511 
512  template<>
513  struct MapVertices< 3, 2 >
514  {
515  static int apply ( int subEntity, int vertex )
516  {
517  assert( (subEntity >= 0) && (subEntity < 6) );
518  assert( (vertex >= 0) && (vertex < 2) );
519  static const int map[ 6 ][ 2 ] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
520  return map[ subEntity ][ vertex ];
521  }
522  };
523 
524  template< int dim >
525  struct MapVertices< dim, dim >
526  {
527  static int apply ( int subEntity, int vertex )
528  {
529  assert( (subEntity >= 0) && (subEntity < NumSubEntities< dim, 1 >::value) );
530  assert( vertex == 0 );
531  return subEntity;
532  }
533  };
534 
535 
536 
537  // Twist
538  // -----
539 
540  // ******************************************************************
541  // Meaning of the twist (same as in ALU)
542  // -------------------------------------
543  //
544  // Consider a fixed ordering of the vertices v_1, ... v_n of a face
545  // (here, we assume their indices to be increasing). Denote by k the
546  // local number of a vertex v within the element and by t the twist.
547  // Then, v = v_j, where j is computed by the following formula:
548  //
549  // / (2n + 1 - k + t) % n, if t < 0
550  // j = <
551  // \ (k + t) % n, if t >= 0
552  //
553  // Note: We use the order of the 0-th vertex dof to assign the twist.
554  // This is ok for two reasons:
555  // - ALBERTA preserves the relative order of the dofs during
556  // dof compression.
557  // - ALBERTA enforces the first vertex dof admin to be periodic.
558  // ******************************************************************
559 
560  template< int dim, int subdim >
561  struct Twist
562  {
563  static const int numSubEntities = NumSubEntities< dim, dim-subdim >::value;
564 
565  static const int minTwist = 0;
566  static const int maxTwist = 0;
567 
568  static int twist ( const Element *element, int subEntity )
569  {
570  assert( (subEntity >= 0) && (subEntity < numSubEntities) );
571  return 0;
572  }
573  };
574 
575  template< int dim >
576  struct Twist< dim, 1 >
577  {
578  static const int numSubEntities = NumSubEntities< dim, dim-1 >::value;
579 
580  static const int minTwist = 0;
581  static const int maxTwist = 1;
582 
583  static int twist ( const Element *element, int subEntity )
584  {
585  assert( (subEntity >= 0) && (subEntity < numSubEntities) );
586  const int numVertices = NumSubEntities< 1, 1 >::value;
587  int dof[ numVertices ];
588  for( int i = 0; i < numVertices; ++i )
589  {
590  const int j = MapVertices< dim, dim-1 >::apply( subEntity, i );
591  dof[ i ] = element->dof[ j ][ 0 ];
592  }
593  return (dof[ 0 ] < dof[ 1 ] ? 0 : 1);
594  }
595  };
596 
597 
598  template<>
599  struct Twist< 1, 1 >
600  {
601  static const int minTwist = 0;
602  static const int maxTwist = 0;
603 
604  static int twist ( const Element *element, int subEntity )
605  {
606  assert( subEntity == 0 );
607  return 0;
608  }
609  };
610 
611 
612  template< int dim >
613  struct Twist< dim, 2 >
614  {
615  static const int numSubEntities = NumSubEntities< dim, dim-2 >::value;
616 
617  static const int minTwist = -3;
618  static const int maxTwist = 2;
619 
620  static int twist ( const Element *element, int subEntity )
621  {
622  assert( (subEntity >= 0) && (subEntity < numSubEntities) );
623  const int numVertices = NumSubEntities< 2, 2 >::value;
624  int dof[ numVertices ];
625  for( int i = 0; i < numVertices; ++i )
626  {
627  const int j = MapVertices< dim, dim-2 >::apply( subEntity, i );
628  dof[ i ] = element->dof[ j ][ 0 ];
629  }
630 
631  const int twist[ 8 ] = { -2, 1, 666, -1, 2, 666, -3, 0 };
632  const int k = int( dof[ 0 ] < dof[ 1 ] )
633  | (int( dof[ 0 ] < dof[ 2 ] ) << 1)
634  | (int( dof[ 1 ] < dof[ 2 ] ) << 2);
635  assert( twist[ k ] != 666 );
636  return twist[ k ];
637  }
638  };
639 
640 
641  template<>
642  struct Twist< 2, 2 >
643  {
644  static const int minTwist = 0;
645  static const int maxTwist = 0;
646 
647  static int twist ( const Element *element, int subEntity )
648  {
649  assert( subEntity == 0 );
650  return 0;
651  }
652  };
653 
654 
655 
656  template< int dim >
657  inline int applyTwist ( int twist, int i )
658  {
659  const int numCorners = NumSubEntities< dim, dim >::value;
660  return (twist < 0 ? (2*numCorners + 1 - i + twist) : i + twist) % numCorners;
661  }
662 
663  template< int dim >
664  inline int applyInverseTwist ( int twist, int i )
665  {
666  const int numCorners = NumSubEntities< dim, dim >::value;
667  return (twist < 0 ? (2*numCorners + 1 - i + twist) : numCorners + i - twist) % numCorners;
668  }
669 
670  }
671 
672 }
673 
674 #endif // #if HAVE_ALBERTA
675 
676 #endif // #ifndef DUNE_ALBERTA_MISC_HH