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