5#ifndef DUNE_ALBERTA_MISC_HH
6#define DUNE_ALBERTA_MISC_HH
12#include <dune/common/hybridutilities.hh>
15#include <dune/grid/albertagrid/albertaheader.hh>
20#ifndef DUNE_ALBERTA_CACHE_COORDINATES
21#define DUNE_ALBERTA_CACHE_COORDINATES 1
46 static const int dimWorld = DIM_OF_WORLD;
48 typedef ALBERTA REAL Real;
49 typedef ALBERTA REAL_B LocalVector;
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;
56 static const int meshRefined = MESH_REFINED;
57 static const int meshCoarsened = MESH_COARSENED;
59 static const int InteriorBoundary = INTERIOR;
60 static const int DirichletBoundary = DIRICHLET;
61 typedef ALBERTA BNDRY_TYPE BoundaryId;
63 typedef U_CHAR ElementType;
65 typedef ALBERTA FE_SPACE DofSpace;
72 template<
class Data >
73 inline Data *memAlloc (
size_t size )
75 return MEM_ALLOC(
size, Data );
78 template<
class Data >
79 inline Data *memCAlloc (
size_t size )
81 return MEM_CALLOC(
size, Data );
84 template<
class Data >
85 inline Data *memReAlloc ( Data *ptr,
size_t oldSize,
size_t newSize )
87 return MEM_REALLOC( ptr, oldSize, newSize, Data );
90 template<
class Data >
91 inline void memFree ( Data *ptr,
size_t size )
93 return MEM_FREE( ptr,
size, Data );
103 typedef GlobalSpace This;
106 typedef GlobalMatrix Matrix;
107 typedef GlobalVector Vector;
110 Matrix identityMatrix_;
115 for(
int i = 0; i < dimWorld; ++i )
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 );
124 static This &instance ()
126 static This theInstance;
131 static const Matrix &identityMatrix ()
133 return instance().identityMatrix_;
136 static const Vector &nullVector ()
138 return instance().nullVector_;
147 template<
int dim,
int codim >
148 struct NumSubEntities;
151 struct NumSubEntities< dim, 0 >
153 static const int value = 1;
157 struct NumSubEntities< dim, dim >
159 static const int value = dim+1;
163 struct NumSubEntities< 0, 0 >
165 static const int value = 1;
169 struct NumSubEntities< 2, 1 >
171 static const int value = 3;
175 struct NumSubEntities< 3, 1 >
177 static const int value = 4;
181 struct NumSubEntities< 3, 2 >
183 static const int value = 6;
191 template<
int dim,
int codim >
195 struct CodimType< dim, 0 >
197 static const int value = CENTER;
201 struct CodimType< dim, dim >
203 static const int value = VERTEX;
207 struct CodimType< 2, 1 >
209 static const int value = EDGE;
213 struct CodimType< 3, 1 >
215 static const int value = FACE;
219 struct CodimType< 3, 2 >
221 static const int value = EDGE;
232 typedef ALBERTA FLAGS Flags;
234 static const Flags nothing = FILL_NOTHING;
236 static const Flags coords = FILL_COORDS;
238 static const Flags neighbor = FILL_NEIGH;
240 static const Flags orientation = (dim == 3 ? FILL_ORIENTATION : FILL_NOTHING);
242 static const Flags projection = FILL_PROJECTION;
244 static const Flags elementType = FILL_NOTHING;
246 static const Flags boundaryId = FILL_MACRO_WALLS;
248 static const Flags nonPeriodic = FILL_NON_PERIODIC;
250 static const Flags
all = coords | neighbor | boundaryId | nonPeriodic
251 | orientation | projection | elementType;
253 static const Flags standardWithCoords =
all & ~nonPeriodic & ~projection;
255#if DUNE_ALBERTA_CACHE_COORDINATES
256 static const Flags standard = standardWithCoords & ~coords;
258 static const Flags standard = standardWithCoords;
268 struct RefinementEdge
270 static const int value = 0;
274 struct RefinementEdge< 2 >
276 static const int value = 2;
284 template<
int dim,
int codim >
285 struct Dune2AlbertaNumbering
287 static int apply (
const int i )
289 assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
295 struct Dune2AlbertaNumbering< 3, 2 >
297 static const int numSubEntities = NumSubEntities< 3, 2 >::value;
299 static int apply (
const int i )
301 assert( (i >= 0) && (i < numSubEntities) );
302 static int dune2alberta[ numSubEntities ] = { 0, 3, 1, 2, 4, 5 };
303 return dune2alberta[ i ];
312 template<
int dim,
int codim >
313 struct Generic2AlbertaNumbering
315 static int apply (
const int i )
317 assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
323 struct Generic2AlbertaNumbering< dim, 1 >
325 static int apply (
const int i )
327 assert( (i >= 0) && (i < NumSubEntities< dim, 1 >::value) );
333 struct Generic2AlbertaNumbering< 1, 1 >
335 static int apply (
const int i )
337 assert( (i >= 0) && (i < NumSubEntities< 1, 1 >::value) );
343 struct Generic2AlbertaNumbering< 3, 2 >
345 static const int numSubEntities = NumSubEntities< 3, 2 >::value;
347 static int apply (
const int i )
349 assert( (i >= 0) && (i < numSubEntities) );
350 static int generic2alberta[ numSubEntities ] = { 0, 1, 3, 2, 4, 5 };
351 return generic2alberta[ i ];
360 template<
int dim,
template<
int,
int >
class Numbering = Generic2AlbertaNumbering >
363 typedef NumberingMap< dim, Numbering > This;
365 template<
int codim >
368 int *dune2alberta_[ dim+1 ];
369 int *alberta2dune_[ dim+1 ];
370 int numSubEntities_[ dim+1 ];
372 NumberingMap (
const This & );
373 This &operator= (
const This & );
378 Hybrid::forEach( std::make_index_sequence< dim+1 >{}, [ & ](
auto i ){ Initialize< i >::apply( *
this ); } );
383 for(
int codim = 0; codim <= dim; ++codim )
385 delete[]( dune2alberta_[ codim ] );
386 delete[]( alberta2dune_[ codim ] );
390 int dune2alberta (
int codim,
int i )
const
392 assert( (codim >= 0) && (codim <= dim) );
393 assert( (i >= 0) && (i < numSubEntities( codim )) );
394 return dune2alberta_[ codim ][ i ];
397 int alberta2dune (
int codim,
int i )
const
399 assert( (codim >= 0) && (codim <= dim) );
400 assert( (i >= 0) && (i < numSubEntities( codim )) );
401 return alberta2dune_[ codim ][ i ];
404 int numSubEntities (
int codim )
const
406 assert( (codim >= 0) && (codim <= dim) );
407 return numSubEntities_[ codim ];
416 template<
int dim,
template<
int,
int >
class Numbering >
417 template<
int codim >
418 struct NumberingMap< dim, Numbering >::Initialize
420 static const int numSubEntities = NumSubEntities< dim, codim >::value;
422 static void apply ( NumberingMap< dim, Numbering > &map )
424 map.numSubEntities_[ codim ] = numSubEntities;
425 map.dune2alberta_[ codim ] =
new int[ numSubEntities ];
426 map.alberta2dune_[ codim ] =
new int[ numSubEntities ];
428 for(
int i = 0; i < numSubEntities; ++i )
430 const int j = Numbering< dim, codim >::apply( i );
431 map.dune2alberta_[ codim ][ i ] = j;
432 map.alberta2dune_[ codim ][ j ] = i;
442 template<
int dim,
int codim >
446 struct MapVertices< dim, 0 >
448 static int apply (
int subEntity,
int vertex )
450 assert( subEntity == 0 );
451 assert( (
vertex >= 0) && (
vertex <= NumSubEntities< dim, dim >::value) );
457 struct MapVertices< 2, 1 >
459 static int apply (
int subEntity,
int vertex )
461 assert( (subEntity >= 0) && (subEntity < 3) );
464 static const int map[ 3 ][ 2 ] = { {1,2}, {0,2}, {0,1} };
465 return map[ subEntity ][
vertex ];
470 struct MapVertices< 3, 1 >
472 static int apply (
int subEntity,
int vertex )
474 assert( (subEntity >= 0) && (subEntity < 4) );
477 static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,2,3}, {0,1,3}, {0,1,2} };
478 return map[ subEntity ][
vertex ];
483 struct MapVertices< 3, 2 >
485 static int apply (
int subEntity,
int vertex )
487 assert( (subEntity >= 0) && (subEntity < 6) );
489 static const int map[ 6 ][ 2 ] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
490 return map[ subEntity ][
vertex ];
495 struct MapVertices< dim, dim >
497 static int apply (
int subEntity,
int vertex )
499 assert( (subEntity >= 0) && (subEntity < NumSubEntities< dim, 1 >::value) );
530 template<
int dim,
int subdim >
533 static const int numSubEntities = NumSubEntities< dim, dim-subdim >::value;
535 static const int minTwist = 0;
536 static const int maxTwist = 0;
538 static int twist ( [[maybe_unused]]
const Element *element,
539 [[maybe_unused]]
int subEntity )
541 assert( (subEntity >= 0) && (subEntity < numSubEntities) );
547 struct Twist< dim, 1 >
549 static const int numSubEntities = NumSubEntities< dim, dim-1 >::value;
551 static const int minTwist = 0;
552 static const int maxTwist = 1;
554 static int twist (
const Element *element,
int subEntity )
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 )
561 const int j = MapVertices< dim, dim-1 >::apply( subEntity, i );
562 dof[ i ] = element->dof[ j ][ 0 ];
564 return (dof[ 0 ] < dof[ 1 ] ? 0 : 1);
572 static const int minTwist = 0;
573 static const int maxTwist = 0;
575 static int twist ( [[maybe_unused]]
const Element *element,
576 [[maybe_unused]]
int subEntity )
578 assert( subEntity == 0 );
585 struct Twist< dim, 2 >
587 static const int numSubEntities = NumSubEntities< dim, dim-2 >::value;
589 static const int minTwist = -3;
590 static const int maxTwist = 2;
592 static int twist (
const Element *element,
int subEntity )
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 )
599 const int j = MapVertices< dim, dim-2 >::apply( subEntity, i );
600 dof[ i ] = element->dof[ j ][ 0 ];
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 );
616 static const int minTwist = 0;
617 static const int maxTwist = 0;
619 static int twist ( [[maybe_unused]]
const Element *element,
620 [[maybe_unused]]
int subEntity )
622 assert( subEntity == 0 );
630 inline int applyTwist (
int twist,
int i )
632 const int numCorners = NumSubEntities< dim, dim >::value;
633 return (twist < 0 ? (2*numCorners + 1 - i + twist) : i + twist) % numCorners;
637 inline int applyInverseTwist (
int twist,
int i )
639 const int numCorners = NumSubEntities< dim, dim >::value;
640 return (twist < 0 ? (2*numCorners + 1 - i + twist) : numCorners + i - twist) % numCorners;
A few common exception classes.
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr All all
PartitionSet for all partitions.
Definition: partitionset.hh:295
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Traits for type conversions and type information.