dune-grid  2.2.0
meshpointer.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALBERTA_MESHPOINTER_HH
00002 #define DUNE_ALBERTA_MESHPOINTER_HH
00003 
00009 #include <limits>
00010 #include <string>
00011 
00012 #include <dune/grid/albertagrid/misc.hh>
00013 #include <dune/grid/albertagrid/elementinfo.hh>
00014 #include <dune/grid/albertagrid/macrodata.hh>
00015 #include <dune/grid/albertagrid/projection.hh>
00016 
00017 #if HAVE_ALBERTA
00018 
00019 namespace Dune
00020 {
00021 
00022   namespace Alberta
00023   {
00024 
00025     // External Forward Declarations
00026     // -----------------------------
00027 
00028     template< int dim >
00029     class HierarchyDofNumbering;
00030 
00031 
00032 
00033     // MeshPointer
00034     // -----------
00035 
00036     template< int dim >
00037     class MeshPointer
00038     {
00039       typedef Alberta::ElementInfo< dim > ElementInfo;
00040       typedef typename ElementInfo::MacroElement MacroElement;
00041       typedef typename ElementInfo::FillFlags FillFlags;
00042 
00043       class BoundaryProvider;
00044 
00045       template< int dimWorld >
00046       struct Library;
00047 
00048     public:
00049       class MacroIterator;
00050 
00051       MeshPointer ()
00052       : mesh_( 0 )
00053       {}
00054 
00055       explicit MeshPointer ( Mesh *mesh )
00056       : mesh_( mesh )
00057       {}
00058 
00059       operator Mesh * () const
00060       {
00061         return mesh_;
00062       }
00063 
00064       operator bool () const
00065       {
00066         return (bool)mesh_;
00067       }
00068 
00069       MacroIterator begin () const
00070       {
00071         return MacroIterator( *this, false );
00072       }
00073 
00074       MacroIterator end () const
00075       {
00076         return MacroIterator( *this, true );
00077       }
00078 
00079       int numMacroElements () const;
00080       int size ( int codim ) const;
00081 
00082       // create a mesh from a macrodata structure
00083       // params:  macroData - macro data structure
00084       // returns: number of boundary segments
00085       unsigned int create ( const MacroData< dim > &macroData );
00086 
00087       // create a mesh from a macrodata structure, adding projections
00088       // params:  macroData         - macro data structure
00089       //          projectionFactory - factory for the projections
00090       // returns: number of boundary segments
00091       template< class Proj, class Impl >
00092       unsigned int create ( const MacroData< dim > &macroData,
00093                             const ProjectionFactoryInterface< Proj, Impl > &projectionFactory );
00094 
00095       // create a mesh from a file
00096       // params:  filename - file name of an Alberta macro triangulation
00097       //          binary   - read binary?
00098       // returns: number of boundary segments
00099       unsigned int create ( const std::string &filename, bool binary = false );
00100 
00101       // read back a mesh from a file
00102       // params:  filename - file name of an Alberta save file
00103       //          time     - variable to receive the time stored in the file
00104       // returns: number of boundary segments
00105       //
00106       // notes: - projections are not preserved
00107       //        - we assume that projections are added in the same order they
00108       //          inserted in when the grid was created (otherwise the boundary
00109       //          indices change)
00110       unsigned int read ( const std::string &filename, Real &time );
00111 
00112       bool write ( const std::string &filename, Real time ) const;
00113 
00114       void release ();
00115 
00116       template< class Functor >
00117       void hierarchicTraverse ( Functor &functor,
00118                                 typename FillFlags::Flags fillFlags = FillFlags::standard ) const;
00119 
00120       template< class Functor >
00121       void leafTraverse ( Functor &functor,
00122                           typename FillFlags::Flags fillFlags = FillFlags::standard ) const;
00123 
00124       bool coarsen ( typename FillFlags::Flags fillFlags = FillFlags::nothing );
00125 
00126       bool refine ( typename FillFlags::Flags fillFlags = FillFlags::nothing );
00127 
00128     private:
00129       static ALBERTA NODE_PROJECTION *
00130       initNodeProjection ( Mesh *mesh, ALBERTA MACRO_EL *macroElement, int n );
00131       template< class ProjectionProvider >
00132       static ALBERTA NODE_PROJECTION *
00133       initNodeProjection ( Mesh *mesh, ALBERTA MACRO_EL *macroElement, int n );
00134 
00135       Mesh *mesh_;
00136     };
00137 
00138 
00139 
00140     // MeshPointer::Library
00141     // --------------------
00142 
00143     template< int dim >
00144     template< int dimWorld >
00145     struct MeshPointer< dim >::Library
00146     {
00147       typedef Alberta::MeshPointer< dim > MeshPointer;
00148 
00149       static unsigned int boundaryCount;
00150       static const void *projectionFactory;
00151 
00152       static void
00153       create ( MeshPointer &ptr, const MacroData< dim > &macroData,
00154                ALBERTA NODE_PROJECTION *(*initNodeProjection)( Mesh *, ALBERTA MACRO_EL *, int ) );
00155       static void release ( MeshPointer &ptr );
00156     };
00157 
00158 
00159 
00160     // MeshPointer::MacroIterator
00161     // --------------------------
00162 
00163     template< int dim >
00164     class MeshPointer< dim >::MacroIterator
00165     {
00166       typedef MacroIterator This;
00167 
00168       friend class MeshPointer< dim >;
00169 
00170     public:
00171       typedef Alberta::MeshPointer< dim > MeshPointer;
00172       typedef Alberta::ElementInfo< dim > ElementInfo;
00173 
00174     private:
00175       explicit MacroIterator ( const MeshPointer &mesh, bool end = false )
00176       : mesh_( mesh ),
00177         index_( end ? mesh.numMacroElements() : 0 )
00178       {}
00179 
00180     public:
00181       bool done () const
00182       {
00183         return (index_ >= mesh().numMacroElements());
00184       }
00185 
00186       bool equals ( const MacroIterator &other ) const
00187       {
00188         return (index_ == other.index_);
00189       }
00190 
00191       void increment ()
00192       {
00193         assert( !done() );
00194         ++index_;
00195       }
00196 
00197       const MacroElement &macroElement () const
00198       {
00199         assert( !done() );
00200         return static_cast< const MacroElement & >( mesh().mesh_->macro_els[ index_ ] );
00201       }
00202 
00203       const MeshPointer &mesh () const
00204       {
00205         return mesh_;
00206       }
00207 
00208       This &operator++ ()
00209       {
00210         increment();
00211         return *this;
00212       }
00213 
00214       ElementInfo operator* () const
00215       {
00216         return elementInfo();
00217       }
00218 
00219       bool operator== ( const MacroIterator &other ) const
00220       {
00221         return equals( other );
00222       }
00223 
00224       bool operator!= ( const MacroIterator &other ) const
00225       {
00226         return !equals( other );
00227       }
00228 
00229       ElementInfo
00230       elementInfo ( typename FillFlags::Flags fillFlags = FillFlags::standard ) const
00231       {
00232         if( done() )
00233           return ElementInfo();
00234         else
00235           return ElementInfo( mesh(), macroElement(), fillFlags );
00236       }
00237 
00238     private:
00239       MeshPointer mesh_;
00240       int index_;
00241     };
00242     
00243     
00244     
00245     // Implementation of MeshPointer
00246     // -----------------------------
00247 
00248     template< int dim >
00249     inline int MeshPointer< dim >::numMacroElements () const
00250     {
00251       return (mesh_ ? mesh_->n_macro_el : 0);
00252     }
00253 
00254 
00255     template<>
00256     inline int MeshPointer< 1 >::size( int codim ) const
00257     {
00258       assert( (codim >= 0) && (codim <= 1) );
00259       return (codim == 0 ? mesh_->n_elements : mesh_->n_vertices);
00260     }
00261 
00262     template<>
00263     inline int MeshPointer< 2 >::size( int codim ) const
00264     {
00265       assert( (codim >= 0) && (codim <= 2) );
00266       if( codim == 0 )
00267         return mesh_->n_elements;
00268       else if( codim == 2 )
00269         return mesh_->n_vertices;
00270       else
00271         return mesh_->n_edges;
00272     }
00273 
00274     template<>
00275     inline int MeshPointer< 3 >::size( int codim ) const
00276     {
00277       assert( (codim >= 0) && (codim <= 3) );
00278       if( codim == 0 )
00279         return mesh_->n_elements;
00280       else if( codim == 3 )
00281         return mesh_->n_vertices;
00282       else if( codim == 1 )
00283         return mesh_->n_faces;
00284       else
00285         return mesh_->n_edges;
00286     }
00287 
00288 
00289     template< int dim >
00290     inline unsigned int MeshPointer< dim >
00291       ::create ( const MacroData< dim > &macroData )
00292     {
00293       release();
00294 
00295       Library< dimWorld >::boundaryCount = 0;
00296       Library< dimWorld >::create( *this, macroData, &initNodeProjection );
00297       return Library< dimWorld >::boundaryCount;
00298     }
00299 
00300 
00301     template< int dim >
00302     template< class Proj, class Impl >
00303     inline unsigned int MeshPointer< dim >
00304       ::create ( const MacroData< dim > &macroData,
00305                  const ProjectionFactoryInterface< Proj, Impl > &projectionFactory )
00306     {
00307       typedef ProjectionFactoryInterface< Proj, Impl > ProjectionFactory;
00308 
00309       release();
00310 
00311       Library< dimWorld >::boundaryCount = 0;
00312       Library< dimWorld >::projectionFactory = &projectionFactory;
00313       Library< dimWorld >::create( *this, macroData, &initNodeProjection< ProjectionFactory > );
00314       Library< dimWorld >::projectionFactory = 0;
00315       return Library< dimWorld >::boundaryCount;
00316     }
00317 
00318 
00319 
00320 
00321     template< int dim >
00322     inline unsigned int MeshPointer< dim >
00323       ::create ( const std::string &filename, bool binary )
00324     {
00325       MacroData< dim > macroData;
00326       macroData.read( filename, binary );
00327       const unsigned int boundaryCount = create( macroData );
00328       macroData.release();
00329       return boundaryCount;
00330     }
00331 
00332 
00333     template< int dim >
00334     inline unsigned int MeshPointer< dim >::read ( const std::string &filename, Real &time )
00335     {
00336       release();
00337 
00338       Library< dimWorld >::boundaryCount = 0;
00339 #if DUNE_ALBERTA_VERSION >= 0x300
00340       mesh_ = ALBERTA read_mesh_xdr( filename.c_str(), &time, NULL, NULL );
00341 #else
00342       mesh_ = ALBERTA read_mesh_xdr( filename.c_str(), &time, NULL );
00343 #endif
00344       return Library< dimWorld >::boundaryCount;
00345     }
00346 
00347 
00348     template< int dim >
00349     inline bool MeshPointer< dim >::write ( const std::string &filename, Real time ) const
00350     {
00351       int success = ALBERTA write_mesh_xdr( mesh_, filename.c_str(), time );
00352       return (success == 0);
00353     }
00354 
00355 
00356     template< int dim >
00357     inline void MeshPointer< dim >::release ()
00358     {
00359       Library< dimWorld >::release( *this );
00360     }
00361 
00362 
00363     template< int dim >
00364     template< class Functor >
00365     inline void MeshPointer< dim >
00366       ::hierarchicTraverse ( Functor &functor,
00367                             typename FillFlags::Flags fillFlags ) const
00368     {
00369       const MacroIterator eit = end();
00370       for( MacroIterator it = begin(); it != eit; ++it )
00371       {
00372         const ElementInfo info = it.elementInfo( fillFlags );
00373         info.hierarchicTraverse( functor );
00374       }
00375     }
00376 
00377 
00378     template< int dim >
00379     template< class Functor >
00380     inline void MeshPointer< dim >
00381       ::leafTraverse ( Functor &functor,
00382                        typename FillFlags::Flags fillFlags ) const
00383     {
00384       const MacroIterator eit = end();
00385       for( MacroIterator it = begin(); it != eit; ++it )
00386       {
00387         const ElementInfo info = it.elementInfo( fillFlags );
00388         info.leafTraverse( functor );
00389       }
00390     }
00391 
00392 
00393 #if DUNE_ALBERTA_VERSION >= 0x300
00394     template< int dim >
00395     inline bool MeshPointer< dim >::coarsen ( typename FillFlags::Flags fillFlags )
00396     {
00397       const bool coarsened = (ALBERTA coarsen( mesh_, fillFlags ) == meshCoarsened);
00398       if( coarsened )
00399         ALBERTA dof_compress( mesh_ );
00400       return coarsened;
00401     }
00402 #endif // #if DUNE_ALBERTA_VERSION >= 0x300
00403 
00404 #if DUNE_ALBERTA_VERSION <= 0x200
00405     template< int dim >
00406     inline bool MeshPointer< dim >::coarsen ( typename FillFlags::Flags fillFlags )
00407     {
00408       assert( fillFlags == FillFlags::nothing );
00409       const bool coarsened = (ALBERTA coarsen( mesh_ ) == meshCoarsened);
00410       if( coarsened )
00411         ALBERTA dof_compress( mesh_ );
00412       return coarsened;
00413     }
00414 #endif // #if DUNE_ALBERTA_VERSION <= 0x200
00415 
00416 
00417 #if DUNE_ALBERTA_VERSION >= 0x300
00418     template< int dim >
00419     inline bool MeshPointer< dim >::refine ( typename FillFlags::Flags fillFlags )
00420     {
00421       return (ALBERTA refine( mesh_, fillFlags ) == meshRefined);
00422     }
00423 #endif // #if DUNE_ALBERTA_VERSION >= 0x300
00424 
00425 #if DUNE_ALBERTA_VERSION <= 0x200
00426     template< int dim >
00427     inline bool MeshPointer< dim >::refine ( typename FillFlags::Flags fillFlags )
00428     {
00429       assert( fillFlags == FillFlags::nothing );
00430       return (ALBERTA refine( mesh_ ) == meshRefined);
00431     }
00432 #endif // #if DUNE_ALBERTA_VERSION <= 0x200
00433 
00434 
00435     template< int dim >
00436     inline ALBERTA NODE_PROJECTION *
00437     MeshPointer< dim >::initNodeProjection ( Mesh *mesh, ALBERTA MACRO_EL *macroEl, int n )
00438     {
00439       const MacroElement &macroElement = static_cast< const MacroElement & >( *macroEl );
00440       if( (n > 0) && macroElement.isBoundary( n-1 ) )
00441         return new BasicNodeProjection( Library< dimWorld >::boundaryCount++ );
00442       else
00443         return 0;
00444     }
00445 
00446 
00447     template< int dim >
00448     template< class ProjectionFactory >
00449     inline ALBERTA NODE_PROJECTION *
00450     MeshPointer< dim >::initNodeProjection ( Mesh *mesh, ALBERTA MACRO_EL *macroEl, int n )
00451     {
00452       typedef typename ProjectionFactory::Projection Projection;
00453 
00454       const MacroElement &macroElement = static_cast< const MacroElement & >( *macroEl );
00455 
00456       MeshPointer< dim > meshPointer( mesh );
00457       ElementInfo elementInfo( meshPointer, macroElement, FillFlags::standard );
00458       const ProjectionFactory &projectionFactory = *static_cast< const ProjectionFactory * >( Library< dimWorld >::projectionFactory );
00459       if( (n > 0) && macroElement.isBoundary( n-1 ) )
00460       {
00461         const unsigned int boundaryIndex = Library< dimWorld >::boundaryCount++;
00462         if( projectionFactory.hasProjection( elementInfo, n-1 ) )
00463         {
00464           Projection projection = projectionFactory.projection( elementInfo, n-1 );
00465           return new NodeProjection< dim, Projection >( boundaryIndex, projection );
00466         }
00467         else
00468           return new BasicNodeProjection( boundaryIndex );
00469       }
00470       else if( (dim < dimWorld) && (n == 0) )
00471       {
00472         const unsigned int boundaryIndex = std::numeric_limits< unsigned int >::max();
00473         if( projectionFactory.hasProjection( elementInfo ) )
00474         {
00475           Projection projection = projectionFactory.projection( elementInfo );
00476           return new NodeProjection< dim, Projection >( boundaryIndex, projection );
00477         }
00478         else
00479           return 0;
00480       }
00481       else
00482         return 0;
00483     }
00484 
00485   } // namespace Alberta
00486 
00487 } // namespace Dune
00488 
00489 #endif // #if HAVE_ALBERTA
00490 
00491 #endif // #ifndef DUNE_ALBERTA_MESHPOINTER_HH