dune-grid  2.2.0
misc.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALBERTA_MISC_HH
00002 #define DUNE_ALBERTA_MISC_HH
00003 
00004 #include <cassert>
00005 
00006 #include <dune/common/exceptions.hh>
00007 #include <dune/common/typetraits.hh>
00008 #include <dune/common/forloop.hh>
00009 
00010 #include <dune/grid/albertagrid/albertaheader.hh>
00011 
00012 #if HAVE_ALBERTA
00013 
00014 // should the coordinates be cached in a vector (required for ALBERTA 2.0)?
00015 #ifndef DUNE_ALBERTA_CACHE_COORDINATES
00016 #define DUNE_ALBERTA_CACHE_COORDINATES 1
00017 #endif
00018 
00019 // set to 1 to use generic geometries in AlbertaGrid
00020 #ifndef DUNE_ALBERTA_USE_GENERICGEOMETRY
00021 #define DUNE_ALBERTA_USE_GENERICGEOMETRY 0
00022 #endif
00023 
00024 namespace Dune
00025 {
00026 
00027   // Exceptions
00028   // ----------
00029 
00030   class AlbertaError
00031   : public Exception
00032   {};
00033 
00034   class AlbertaIOError
00035   : public IOError
00036   {};
00037 
00038 
00039 
00040   namespace Alberta
00041   {
00042 
00043     // Import Types
00044     // ------------
00045 
00046     static const int dimWorld = DIM_OF_WORLD;
00047 
00048     typedef ALBERTA REAL Real;
00049     typedef ALBERTA REAL_B LocalVector; // in barycentric coordinates
00050     typedef ALBERTA REAL_D GlobalVector;
00051     typedef ALBERTA REAL_DD GlobalMatrix;
00052 
00053 #if DUNE_ALBERTA_VERSION >= 0x300
00054     typedef ALBERTA AFF_TRAFO AffineTransformation;
00055 #else
00056     struct AffineTransformation
00057     {
00058       GlobalMatrix M;
00059       GlobalVector t;
00060     };
00061 #endif
00062 
00063     typedef ALBERTA MESH Mesh;
00064     typedef ALBERTA EL Element;
00065 
00066     static const int meshRefined = MESH_REFINED;
00067     static const int meshCoarsened = MESH_COARSENED;
00068 
00069     static const int InteriorBoundary = INTERIOR;
00070     static const int DirichletBoundary = DIRICHLET;
00071 #if DUNE_ALBERTA_VERSION >= 0x300
00072     typedef ALBERTA BNDRY_TYPE BoundaryId;
00073 #else
00074     typedef S_CHAR BoundaryId;
00075 #endif
00076 
00077     typedef U_CHAR ElementType;
00078 
00079     typedef ALBERTA FE_SPACE DofSpace;
00080 
00081 
00082 
00083     // Memory Manipulation Functions
00084     // -----------------------------
00085 
00086     template< class Data >
00087     inline Data *memAlloc ( size_t size )
00088     {
00089       return MEM_ALLOC( size, Data );
00090     }
00091 
00092     template< class Data >
00093     inline Data *memCAlloc ( size_t size )
00094     {
00095       return MEM_CALLOC( size, Data );
00096     }
00097 
00098     template< class Data >
00099     inline Data *memReAlloc ( Data *ptr, size_t oldSize, size_t newSize )
00100     {
00101       return MEM_REALLOC( ptr, oldSize, newSize, Data );
00102     }
00103 
00104     template< class Data >
00105     inline void memFree ( Data *ptr, size_t size )
00106     {
00107       return MEM_FREE( ptr, size, Data );
00108     }
00109 
00110 
00111 
00112     // GlobalSpace
00113     // -----------
00114 
00115     class GlobalSpace
00116     {
00117       typedef GlobalSpace This;
00118 
00119     public:
00120       typedef GlobalMatrix Matrix;
00121       typedef GlobalVector Vector;
00122 
00123     private:
00124       Matrix identityMatrix_;
00125       Vector nullVector_;
00126 
00127       GlobalSpace ()
00128       {
00129         for( int i = 0; i < dimWorld; ++i )
00130         {
00131           for( int j = 0; j < dimWorld; ++j )
00132             identityMatrix_[ i ][ j ] = Real( 0 );
00133           identityMatrix_[ i ][ i ] = Real( 1 );
00134           nullVector_[ i ] = Real( 0 );
00135         }
00136       }
00137 
00138       static This &instance ()
00139       {
00140         static This theInstance;
00141         return theInstance;
00142       }
00143 
00144     public:
00145       static const Matrix &identityMatrix ()
00146       {
00147         return instance().identityMatrix_;
00148       }
00149 
00150       static const Vector &nullVector ()
00151       {
00152         return instance().nullVector_;
00153       }
00154     };
00155 
00156 
00157 
00158     // NumSubEntities
00159     // --------------
00160 
00161     template< int dim, int codim >
00162     struct NumSubEntities;
00163 
00164     template< int dim >
00165     struct NumSubEntities< dim, 0 >
00166     {
00167       static const int value = 1;
00168     };
00169 
00170     template< int dim >
00171     struct NumSubEntities< dim, dim >
00172     {
00173       static const int value = dim+1;
00174     };
00175 
00176     template<>
00177     struct NumSubEntities< 0, 0 >
00178     {
00179       static const int value = 1;
00180     };
00181 
00182     template<>
00183     struct NumSubEntities< 2, 1 >
00184     {
00185       static const int value = 3;
00186     };
00187 
00188     template<>
00189     struct NumSubEntities< 3, 1 >
00190     {
00191       static const int value = 4;
00192     };
00193 
00194     template<>
00195     struct NumSubEntities< 3, 2 >
00196     {
00197       static const int value = 6;
00198     };
00199 
00200 
00201 
00202     // CodimType
00203     // ---------
00204 
00205     template< int dim, int codim >
00206     struct CodimType;
00207 
00208     template< int dim >
00209     struct CodimType< dim, 0 >
00210     {
00211       static const int value = CENTER;
00212     };
00213 
00214     template< int dim >
00215     struct CodimType< dim, dim >
00216     {
00217       static const int value = VERTEX;
00218     };
00219 
00220     template<>
00221     struct CodimType< 2, 1 >
00222     {
00223       static const int value = EDGE;
00224     };
00225 
00226     template<>
00227     struct CodimType< 3, 1 >
00228     {
00229       static const int value = FACE;
00230     };
00231 
00232     template<>
00233     struct CodimType< 3, 2 >
00234     {
00235       static const int value = EDGE;
00236     };
00237 
00238 
00239 
00240     // FillFlags
00241     // ---------
00242 
00243     template< int dim >
00244     struct FillFlags
00245     {
00246       typedef ALBERTA FLAGS Flags;
00247 
00248       static const Flags nothing = FILL_NOTHING;
00249 
00250       static const Flags coords = FILL_COORDS;
00251 
00252       static const Flags neighbor = FILL_NEIGH;
00253 
00254       static const Flags orientation = (dim == 3 ? FILL_ORIENTATION : FILL_NOTHING);
00255 
00256       static const Flags projection = FILL_PROJECTION;
00257 
00258 #if DUNE_ALBERTA_VERSION >= 0x300
00259       static const Flags elementType = FILL_NOTHING;
00260 #else
00261       static const Flags elementType = (dim == 3 ? FILL_EL_TYPE : FILL_NOTHING);
00262 #endif
00263 
00264 #if DUNE_ALBERTA_VERSION >= 0x300
00265       static const Flags boundaryId = FILL_MACRO_WALLS;
00266 #else
00267       static const Flags boundaryId = FILL_BOUND;
00268 #endif
00269 
00270 #if DUNE_ALBERTA_VERSION >= 0x300
00271       static const Flags nonPeriodic = FILL_NON_PERIODIC;
00272 #else
00273       static const Flags nonPeriodic = FILL_NOTHING;
00274 #endif
00275 
00276       static const Flags all = coords | neighbor | boundaryId | nonPeriodic
00277                                | orientation | projection | elementType;
00278 
00279 #if DUNE_ALBERTA_VERSION >= 0x300
00280       static const Flags standardWithCoords = all & ~nonPeriodic & ~projection;
00281 #else
00282       static const Flags standardWithCoords = all;
00283 #endif
00284 
00285 #if DUNE_ALBERTA_CACHE_COORDINATES
00286       static const Flags standard = standardWithCoords & ~coords;
00287 #else
00288       static const Flags standard = standardWithCoords;
00289 #endif
00290     };
00291 
00292 
00293 
00294     // RefinementEdge
00295     // --------------
00296 
00297     template< int dim >
00298     struct RefinementEdge
00299     {
00300       static const int value = 0;
00301     };
00302 
00303     template<>
00304     struct RefinementEdge< 2 >
00305     {
00306       static const int value = 2;
00307     };
00308 
00309 
00310 
00311     // Dune2AlbertaNumbering
00312     // ---------------------
00313 
00314     template< int dim, int codim >
00315     struct Dune2AlbertaNumbering
00316     {
00317       static int apply ( const int i )
00318       {
00319         assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
00320         return i;
00321       }
00322     };
00323 
00324     template<>
00325     struct Dune2AlbertaNumbering< 3, 2 >
00326     {
00327       static const int numSubEntities = NumSubEntities< 3, 2 >::value;
00328 
00329       static int apply ( const int i )
00330       {
00331         assert( (i >= 0) && (i < numSubEntities) );
00332         static int dune2alberta[ numSubEntities ] = { 0, 3, 1, 2, 4, 5 };
00333         return dune2alberta[ i ];
00334       }
00335     };
00336 
00337 
00338 
00339     // Generic2AlbertaNumbering
00340     // ------------------------
00341 
00342     template< int dim, int codim >
00343     struct Generic2AlbertaNumbering
00344     {
00345       static int apply ( const int i )
00346       {
00347         assert( (i >= 0) && (i < NumSubEntities< dim, codim >::value) );
00348         return i;
00349       }
00350     };
00351 
00352     template< int dim >
00353     struct Generic2AlbertaNumbering< dim, 1 >
00354     {
00355       static int apply ( const int i )
00356       {
00357         assert( (i >= 0) && (i < NumSubEntities< dim, 1 >::value) );
00358         return dim - i;
00359       }
00360     };
00361 
00362     template<>
00363     struct Generic2AlbertaNumbering< 1, 1 >
00364     {
00365       static int apply ( const int i )
00366       {
00367         assert( (i >= 0) && (i < NumSubEntities< 1, 1 >::value) );
00368         return i;
00369       }
00370     };
00371 
00372     template<>
00373     struct Generic2AlbertaNumbering< 3, 2 >
00374     {
00375       static const int numSubEntities = NumSubEntities< 3, 2 >::value;
00376 
00377       static int apply ( const int i )
00378       {
00379         assert( (i >= 0) && (i < numSubEntities) );
00380         static int generic2alberta[ numSubEntities ] = { 0, 1, 3, 2, 4, 5 };
00381         return generic2alberta[ i ];
00382       }
00383     };
00384 
00385 
00386 
00387     // NumberingMap
00388     // ------------
00389 
00390     template< int dim, template< int, int > class Numbering = Generic2AlbertaNumbering >
00391     class NumberingMap
00392     {
00393       typedef NumberingMap< dim, Numbering > This;
00394 
00395       template< int codim >
00396       struct Initialize;
00397 
00398       int *dune2alberta_[ dim+1 ];
00399       int *alberta2dune_[ dim+1 ];
00400       int numSubEntities_[ dim+1 ];
00401 
00402       NumberingMap ( const This & );
00403       This &operator= ( const This & );
00404 
00405     public:
00406       NumberingMap ()
00407       {
00408         ForLoop< Initialize, 0, dim >::apply( *this );
00409       }
00410 
00411       ~NumberingMap ()
00412       {
00413         for( int codim = 0; codim <= dim; ++codim )
00414         {
00415           delete[]( dune2alberta_[ codim ] );
00416           delete[]( alberta2dune_[ codim ] );
00417         }
00418       }
00419 
00420       int dune2alberta ( int codim, int i ) const
00421       {
00422         assert( (codim >= 0) && (codim <= dim) );
00423         assert( (i >= 0) && (i < numSubEntities( codim )) );
00424         return dune2alberta_[ codim ][ i ];
00425       }
00426 
00427       int alberta2dune ( int codim, int i ) const
00428       {
00429         assert( (codim >= 0) && (codim <= dim) );
00430         assert( (i >= 0) && (i < numSubEntities( codim )) );
00431         return alberta2dune_[ codim ][ i ];
00432       }
00433 
00434       int numSubEntities ( int codim ) const
00435       {
00436         assert( (codim >= 0) && (codim <= dim) );
00437         return numSubEntities_[ codim ];
00438       }
00439     };
00440 
00441 
00442 
00443     // NumberingMap::Initialize
00444     // ------------------------
00445 
00446     template< int dim, template< int, int > class Numbering >
00447     template< int codim >
00448     struct NumberingMap< dim, Numbering >::Initialize
00449     {
00450       static const int numSubEntities = NumSubEntities< dim, codim >::value;
00451 
00452       static void apply ( NumberingMap< dim, Numbering > &map )
00453       {
00454         map.numSubEntities_[ codim ] = numSubEntities;
00455         map.dune2alberta_[ codim ] = new int[ numSubEntities ];
00456         map.alberta2dune_[ codim ] = new int[ numSubEntities ];
00457 
00458         for( int i = 0; i < numSubEntities; ++i )
00459         {
00460           const int j = Numbering< dim, codim >::apply( i );
00461           map.dune2alberta_[ codim ][ i ] = j;
00462           map.alberta2dune_[ codim ][ j ] = i;
00463         }
00464       }
00465     };
00466 
00467 
00468 
00469     // MapVertices
00470     // -----------
00471 
00472     template< int dim, int codim >
00473     struct MapVertices;
00474 
00475     template< int dim >
00476     struct MapVertices< dim, 0 >
00477     {
00478       static int apply ( int subEntity, int vertex )
00479       {
00480         assert( subEntity == 0 );
00481         assert( (vertex >= 0) && (vertex <= NumSubEntities< dim, dim >::value) );
00482         return vertex;
00483       }
00484     };
00485 
00486     template<>
00487     struct MapVertices< 2, 1 >
00488     {
00489       static int apply ( int subEntity, int vertex )
00490       {
00491         assert( (subEntity >= 0) && (subEntity < 3) );
00492         assert( (vertex >= 0) && (vertex < 2) );
00493         //static const int map[ 3 ][ 2 ] = { {1,2}, {2,0}, {0,1} };
00494         static const int map[ 3 ][ 2 ] = { {1,2}, {0,2}, {0,1} };
00495         return map[ subEntity ][ vertex ];
00496       }
00497     };
00498 
00499     template<>
00500     struct MapVertices< 3, 1 >
00501     {
00502       static int apply ( int subEntity, int vertex )
00503       {
00504         assert( (subEntity >= 0) && (subEntity < 4) );
00505         assert( (vertex >= 0) && (vertex < 3) );
00506         //static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,3,2}, {0,1,3}, {0,2,1} };
00507         static const int map[ 4 ][ 3 ] = { {1,2,3}, {0,2,3}, {0,1,3}, {0,1,2} };
00508         return map[ subEntity ][ vertex ];
00509       }
00510     };
00511 
00512     template<>
00513     struct MapVertices< 3, 2 >
00514     {
00515       static int apply ( int subEntity, int vertex )
00516       {
00517         assert( (subEntity >= 0) && (subEntity < 6) );
00518         assert( (vertex >= 0) && (vertex < 2) );
00519         static const int map[ 6 ][ 2 ] = { {0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3} };
00520         return map[ subEntity ][ vertex ];
00521       }
00522     };
00523 
00524     template< int dim >
00525     struct MapVertices< dim, dim >
00526     {
00527       static int apply ( int subEntity, int vertex )
00528       {
00529         assert( (subEntity >= 0) && (subEntity < NumSubEntities< dim, 1 >::value) );
00530         assert( vertex == 0 );
00531         return subEntity;
00532       }
00533     };
00534 
00535 
00536 
00537     // Twist
00538     // -----
00539 
00540     // ******************************************************************
00541     // Meaning of the twist (same as in ALU)
00542     // -------------------------------------
00543     //
00544     // Consider a fixed ordering of the vertices v_1, ... v_n of a face
00545     // (here, we assume their indices to be increasing). Denote by k the
00546     // local number of a vertex v within the element and by t the twist.
00547     // Then, v = v_j, where j is computed by the following formula:
00548     //
00549     //        / (2n + 1 - k + t) % n, if t < 0
00550     //   j = <
00551     //        \ (k + t) % n,          if t >= 0
00552     //
00553     //  Note: We use the order of the 0-th vertex dof to assign the twist.
00554     //        This is ok for two reasons:
00555     //        - ALBERTA preserves the relative order of the dofs during
00556     //          dof compression.
00557     //        - ALBERTA enforces the first vertex dof admin to be periodic.
00558     // ******************************************************************
00559 
00560     template< int dim, int subdim >
00561     struct Twist
00562     {
00563       static const int numSubEntities = NumSubEntities< dim, dim-subdim >::value;
00564 
00565       static const int minTwist = 0;
00566       static const int maxTwist = 0;
00567 
00568       static int twist ( const Element *element, int subEntity )
00569       {
00570         assert( (subEntity >= 0) && (subEntity < numSubEntities) );
00571         return 0;
00572       }
00573     };
00574 
00575     template< int dim >
00576     struct Twist< dim, 1 >
00577     {
00578       static const int numSubEntities = NumSubEntities< dim, dim-1 >::value;
00579 
00580       static const int minTwist = 0;
00581       static const int maxTwist = 1;
00582 
00583       static int twist ( const Element *element, int subEntity )
00584       {
00585         assert( (subEntity >= 0) && (subEntity < numSubEntities) );
00586         const int numVertices = NumSubEntities< 1, 1 >::value;
00587         int dof[ numVertices ];
00588         for( int i = 0; i < numVertices; ++i )
00589         {
00590           const int j = MapVertices< dim, dim-1 >::apply( subEntity, i );
00591           dof[ i ] = element->dof[ j ][ 0 ];
00592         }
00593         return (dof[ 0 ] < dof[ 1 ] ? 0 : 1);
00594       }
00595     };
00596 
00597 
00598     template<>
00599     struct Twist< 1, 1 >
00600     {
00601       static const int minTwist = 0;
00602       static const int maxTwist = 0;
00603 
00604       static int twist ( const Element *element, int subEntity )
00605       {
00606         assert( subEntity == 0 );
00607         return 0;
00608       }
00609     };
00610 
00611 
00612     template< int dim >
00613     struct Twist< dim, 2 >
00614     {
00615       static const int numSubEntities = NumSubEntities< dim, dim-2 >::value;
00616 
00617       static const int minTwist = -3;
00618       static const int maxTwist = 2;
00619 
00620       static int twist ( const Element *element, int subEntity )
00621       {
00622         assert( (subEntity >= 0) && (subEntity < numSubEntities) );
00623         const int numVertices = NumSubEntities< 2, 2 >::value;
00624         int dof[ numVertices ];
00625         for( int i = 0; i < numVertices; ++i )
00626         {
00627           const int j = MapVertices< dim, dim-2 >::apply( subEntity, i );
00628           dof[ i ] = element->dof[ j ][ 0 ];
00629         }
00630 
00631         const int twist[ 8 ] = { -2, 1, 666, -1, 2, 666, -3, 0 };
00632         const int k = int( dof[ 0 ] < dof[ 1 ] )
00633                       | (int( dof[ 0 ] < dof[ 2 ] ) << 1)
00634                       | (int( dof[ 1 ] < dof[ 2 ] ) << 2);
00635         assert( twist[ k ] != 666 );
00636         return twist[ k ];
00637       }
00638     };
00639 
00640 
00641     template<>
00642     struct Twist< 2, 2 >
00643     {
00644       static const int minTwist = 0;
00645       static const int maxTwist = 0;
00646 
00647       static int twist ( const Element *element, int subEntity )
00648       {
00649         assert( subEntity == 0 );
00650         return 0;
00651       }
00652     };
00653 
00654 
00655 
00656     template< int dim >
00657     inline int applyTwist ( int twist, int i )
00658     {
00659       const int numCorners = NumSubEntities< dim, dim >::value;
00660       return (twist < 0 ? (2*numCorners + 1 - i + twist) : i + twist) % numCorners;
00661     }
00662 
00663     template< int dim >
00664     inline int applyInverseTwist ( int twist, int i )
00665     {
00666       const int numCorners = NumSubEntities< dim, dim >::value;
00667       return (twist < 0 ? (2*numCorners + 1 - i + twist) : numCorners + i - twist) % numCorners;
00668     }
00669 
00670   }
00671 
00672 }
00673 
00674 #endif // #if HAVE_ALBERTA
00675 
00676 #endif // #ifndef DUNE_ALBERTA_MISC_HH