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