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