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