dune-grid  2.2.0
macrodata.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALBERTA_MACRODATA_HH
00002 #define DUNE_ALBERTA_MACRODATA_HH
00003 
00009 #include <dune/common/fvector.hh>
00010 #include <dune/common/fmatrix.hh>
00011 
00012 #include <dune/grid/albertagrid/misc.hh>
00013 #include <dune/grid/albertagrid/algebra.hh>
00014 #include <dune/grid/albertagrid/albertaheader.hh>
00015 
00016 #if HAVE_ALBERTA
00017 
00018 namespace Dune
00019 {
00020 
00021   namespace Alberta
00022   {
00023 
00024     template< int dim >
00025     class MacroData
00026     {
00027       typedef MacroData< dim > This;
00028 
00029       typedef ALBERTA MACRO_DATA Data;
00030 
00031       static const int dimension = dim;
00032       static const int numVertices = NumSubEntities< dimension, dimension >::value;
00033       static const int numEdges = NumSubEntities< dimension, dimension-1 >::value;
00034 
00035       static const int initialSize = 4096;
00036 
00037     public:
00038       template< int >
00039       struct Library;
00040 
00041       template< int > friend struct InstantiateMacroDataLibrary;
00042 
00043     public:
00044       typedef int ElementId[ numVertices ];
00045 
00046       static const int supportPeriodicity = (DUNE_ALBERTA_VERSION >= 0x300);
00047 
00048       MacroData ()
00049       : data_( NULL ),
00050         vertexCount_( -1 ),
00051         elementCount_( -1 )
00052       {}
00053 
00054       operator Data * () const
00055       {
00056         return data_;
00057       }
00058 
00059       int vertexCount () const
00060       {
00061         return (vertexCount_ < 0 ? data_->n_total_vertices : vertexCount_);
00062       }
00063 
00064       int elementCount () const
00065       {
00066         return (elementCount_ < 0 ? data_->n_macro_elements : elementCount_);
00067       }
00068 
00069       ElementId &element ( int i ) const;
00070       GlobalVector &vertex ( int i ) const;
00071       int &neighbor ( int element, int i ) const;
00072       BoundaryId &boundaryId ( int element, int i ) const;
00073 
00078       void create ();
00079 
00088       void finalize ();
00089 
00098       void markLongestEdge ();
00099 
00108       void setOrientation ( const Real orientation );
00109 
00120       bool checkNeighbors () const;
00121 
00123       void release ()
00124       {
00125         if( data_ != NULL )
00126         {
00127           ALBERTA free_macro_data( data_ );
00128           data_ = NULL;
00129         }
00130         vertexCount_ = elementCount_ = -1;
00131       }
00132 
00138       int insertElement ( const ElementId &id );
00139 
00145       int insertVertex ( const GlobalVector &coords )
00146       {
00147         assert( vertexCount_ >= 0 );
00148         if( vertexCount_ >= data_->n_total_vertices )
00149           resizeVertices( 2*vertexCount_ );
00150         copy( coords, vertex( vertexCount_ ) );
00151         return vertexCount_++;
00152       }
00153 
00159       int insertVertex ( const FieldVector< Real, dimWorld > &coords )
00160       {
00161         assert( vertexCount_ >= 0 );
00162         if( vertexCount_ >= data_->n_total_vertices )
00163           resizeVertices( 2*vertexCount_ );
00164         copy( coords, vertex( vertexCount_ ) );
00165         return vertexCount_++;
00166       }
00167 
00168       void insertWallTrafo ( const GlobalMatrix &m, const GlobalVector &t );
00169       void insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
00170                              const FieldVector< Real, dimWorld > &shift );
00171 
00172       void checkCycles ();
00173 
00174       void read ( const std::string &filename, bool binary = false );
00175 
00176       bool write ( const std::string &filename, bool binary = false ) const
00177       {
00178         if( binary )
00179           return ALBERTA write_macro_data_xdr( data_, filename.c_str() );
00180         else
00181           return ALBERTA write_macro_data( data_, filename.c_str() );
00182       }
00183 
00184     private:
00185       template< class Vector >
00186       void copy ( const Vector &x, GlobalVector &y )
00187       {
00188         for( int i = 0; i < dimWorld; ++i )
00189           y[ i ] = x[ i ];
00190       }
00191 
00192       void resizeElements ( const int newSize );
00193 
00194       void resizeVertices ( const int newSize )
00195       {
00196         const int oldSize = data_->n_total_vertices;
00197         data_->n_total_vertices = newSize;
00198         data_->coords = memReAlloc< GlobalVector >( data_->coords, oldSize, newSize );
00199         assert( (data_->coords != NULL) || (newSize == 0) );
00200       }
00201 
00202     private:
00203       Data *data_;
00204       int vertexCount_;
00205       int elementCount_;
00206     };
00207 
00208 
00209 
00210     // MacroData::Library
00211     // ------------------
00212 
00213     template< int dim >
00214     template< int >
00215     struct MacroData< dim >::Library
00216     {
00217       typedef Alberta::MacroData< dim > MacroData;
00218 
00219       static bool checkNeighbors ( const MacroData &macroData );
00220       static void markLongestEdge ( MacroData &macroData );
00221       static void setOrientation ( MacroData &macroData, const Real orientation );
00222 
00223     private:
00224       static Real edgeLength ( const MacroData &macroData, const ElementId &e, int edge );
00225       static int longestEdge ( const MacroData &macroData, const ElementId &e );
00226 
00227       template< class Type >
00228       static void rotate ( Type *array, int i, int shift );
00229 
00230       static void rotate ( MacroData &macroData, int i, int shift );
00231       static void swap ( MacroData &macroData, int el, int v1, int v2 );
00232     };
00233 
00234 
00235 
00236     // Implementation of MacroData
00237     // ---------------------------
00238 
00239     template< int dim >
00240     inline typename MacroData< dim >::ElementId &
00241     MacroData< dim >::element ( int i ) const
00242     {
00243       assert( (i >= 0) && (i < data_->n_macro_elements) );
00244       const int offset = i * numVertices;
00245       return *reinterpret_cast< ElementId * >( data_->mel_vertices + offset );
00246     }
00247 
00248 
00249     template< int dim >
00250     inline GlobalVector &MacroData< dim >::vertex ( int i ) const
00251     {
00252       assert( (i >= 0) && (i < data_->n_total_vertices) );
00253       return data_->coords[ i ];
00254     }
00255 
00256 
00257     template< int dim >
00258     inline int &MacroData< dim >::neighbor ( int element, int i ) const
00259     {
00260       assert( (element >= 0) && (element < data_->n_macro_elements) );
00261       assert( (i >= 0) && (i < numVertices) );
00262       return data_->neigh[ element*numVertices + i ];
00263     }
00264 
00265 
00266     template< int dim >
00267     inline BoundaryId &MacroData< dim >::boundaryId ( int element, int i ) const
00268     {
00269       assert( (element >= 0) && (element < data_->n_macro_elements) );
00270       assert( (i >= 0) && (i < numVertices) );
00271       return data_->boundary[ element*numVertices + i ];
00272     }
00273 
00274 
00275 #if DUNE_ALBERTA_VERSION >= 0x300
00276     template< int dim >
00277     inline void MacroData< dim >::create ()
00278     {
00279       release();
00280       data_ = ALBERTA alloc_macro_data( dim, initialSize, initialSize );
00281       data_->boundary = memAlloc< BoundaryId >( initialSize*numVertices );
00282       if( dim == 3 )
00283         data_->el_type = memAlloc< ElementType >( initialSize );
00284       vertexCount_ = elementCount_ = 0;
00285       elementCount_ = 0;
00286     }
00287 #endif // #if DUNE_ALBERTA_VERSION >= 0x300
00288 
00289 #if DUNE_ALBERTA_VERSION == 0x200
00290     template< int dim >
00291     inline void MacroData< dim >::create ()
00292     {
00293       release();
00294       data_ = ALBERTA alloc_macro_data( dim, initialSize, initialSize, 0 );
00295       data_->boundary = memAlloc< BoundaryId >( initialSize*numVertices );
00296       if( dim == 3 )
00297         data_->el_type = memAlloc< ElementType >( initialSize );
00298       vertexCount_ = elementCount_ = 0;
00299       elementCount_ = 0;
00300     }
00301 #endif // #if DUNE_ALBERTA_VERSION == 0x200
00302 
00303 
00304     template< int dim >
00305     inline void MacroData< dim >::finalize ()
00306     {
00307       if( (vertexCount_ >= 0) && (elementCount_ >= 0) )
00308       {
00309         resizeVertices( vertexCount_ );
00310         resizeElements( elementCount_ );
00311         ALBERTA compute_neigh_fast( data_ );
00312 
00313         // assign default boundary id (if none is assigned)
00314         for( int element = 0; element < elementCount_; ++element )
00315         {
00316           for( int i = 0; i < numVertices; ++i )
00317           {
00318             BoundaryId &id = boundaryId( element, i );
00319             if( neighbor( element, i ) >= 0 )
00320             {
00321               assert( id == InteriorBoundary );
00322               id = InteriorBoundary;
00323             }
00324             else
00325               id = (id == InteriorBoundary ? DirichletBoundary : id);
00326           }
00327         }
00328 
00329         vertexCount_ = elementCount_ = -1;
00330       }
00331       assert( (vertexCount_ < 0) && (elementCount_ < 0) );
00332     }
00333 
00334 
00335     template< int dim >
00336     inline void MacroData< dim >::markLongestEdge ()
00337     {
00338       Library< dimWorld >::markLongestEdge( *this );
00339     }
00340 
00341 
00342     template< int dim >
00343     inline void MacroData< dim >::setOrientation ( const Real orientation )
00344     {
00345       Library< dimWorld >::setOrientation( *this, orientation );
00346     }
00347 
00348 
00349     template< int dim >
00350     inline bool MacroData< dim >::checkNeighbors () const
00351     {
00352       return Library< dimWorld >::checkNeighbors( *this );
00353     }
00354 
00355 
00356     template< int dim >
00357     inline int MacroData< dim >::insertElement ( const ElementId &id )
00358     {
00359       assert( elementCount_ >= 0 );
00360       if( elementCount_ >= data_->n_macro_elements )
00361         resizeElements( 2*elementCount_ );
00362 
00363       ElementId &e = element( elementCount_ );
00364       for( int i = 0; i < numVertices; ++i )
00365       {
00366         e[ i ] = id[ i ];
00367         boundaryId( elementCount_, i ) = InteriorBoundary;
00368       }
00369       if( dim == 3 )
00370         data_->el_type[ elementCount_ ] = 0;
00371 
00372       return elementCount_++;
00373     }
00374 
00375 
00376 #if DUNE_ALBERTA_VERSION >= 0x300
00377     template< int dim >
00378     inline void MacroData< dim >
00379       ::insertWallTrafo ( const GlobalMatrix &matrix, const GlobalVector &shift )
00380     {
00381       int &count = data_->n_wall_trafos;
00382       AffineTransformation *&array = data_->wall_trafos;
00383 
00384       // resize wall trafo array
00385       array = memReAlloc< AffineTransformation >( array, count, count+1 );
00386       assert( data_->wall_trafos != NULL );
00387 
00388       // copy matrix and shift
00389       for( int i = 0; i < dimWorld; ++i )
00390         copy( matrix[ i ], array[ count ].M[ i ] );
00391       copy( shift, array[ count ].t );
00392       ++count;
00393     }
00394 
00395     template< int dim >
00396     inline void MacroData< dim >
00397       ::insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
00398                           const FieldVector< Real, dimWorld > &shift )
00399     {
00400       int &count = data_->n_wall_trafos;
00401       AffineTransformation *&array = data_->wall_trafos;
00402 
00403       // resize wall trafo array
00404       array = memReAlloc< AffineTransformation >( array, count, count+1 );
00405       assert( data_->wall_trafos != NULL );
00406 
00407       // copy matrix and shift
00408       for( int i = 0; i < dimWorld; ++i )
00409         copy( matrix[ i ], array[ count ].M[ i ] );
00410       copy( shift, array[ count ].t );
00411       ++count;
00412     }
00413 #endif // #if DUNE_ALBERTA_VERSION >= 0x300
00414 
00415 #if DUNE_ALBERTA_VERSION <= 0x200
00416     template< int dim >
00417     inline void MacroData< dim >
00418       ::insertWallTrafo ( const GlobalMatrix &m, const GlobalVector &t )
00419     {
00420       DUNE_THROW( AlbertaError,
00421                   "Periodic grids are only supported in ALBERTA 3.0 or higher." );
00422     }
00423 
00424     template< int dim >
00425     inline void MacroData< dim >
00426       ::insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
00427                           const FieldVector< Real, dimWorld > &shift )
00428     {
00429       DUNE_THROW( AlbertaError,
00430                   "Periodic grids are only supported in ALBERTA 3.0 or higher." );
00431     }
00432 #endif // #if DUNE_ALBERTA_VERSION <= 0x200
00433 
00434 
00435     template< int dim >
00436     inline void MacroData< dim >::checkCycles ()
00437     {
00438       // ensure that the macro data has been finalized
00439       finalize();
00440       ALBERTA macro_test( data_, NULL );
00441     }
00442 
00443 
00444     template< int dim >
00445     inline void MacroData< dim >::read ( const std::string &filename, bool binary )
00446     {
00447       release();
00448       if( binary )
00449         data_ = ALBERTA read_macro_xdr( filename.c_str() );
00450       else
00451         data_ = ALBERTA read_macro( filename.c_str() );
00452     }
00453 
00454 
00455     template< int dim >
00456     inline void MacroData< dim >::resizeElements ( const int newSize )
00457     {
00458       const int oldSize = data_->n_macro_elements;
00459       data_->n_macro_elements = newSize;
00460       data_->mel_vertices = memReAlloc( data_->mel_vertices, oldSize*numVertices, newSize*numVertices );
00461       data_->boundary = memReAlloc( data_->boundary, oldSize*numVertices, newSize*numVertices );
00462       if( dim == 3 )
00463         data_->el_type = memReAlloc( data_->el_type, oldSize, newSize );
00464       assert( (newSize == 0) || (data_->mel_vertices != NULL) );
00465     }
00466 
00467   }
00468 
00469 }
00470 
00471 #endif // #if HAVE_ALBERTA
00472 
00473 #endif