dune-grid  2.2.0
geostorage.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH
00002 #define DUNE_ALUGRIDGEOMETRYSTORAGE_HH
00003 
00004 // Dune includes
00005 #include <dune/common/misc.hh>
00006 #include <dune/grid/common/grid.hh>
00007 #include <dune/grid/common/gridfactory.hh>
00008 
00009 #if HAVE_ALUGRID
00010 
00011 #include <dune/grid/alugrid/common/declaration.hh>
00012 #include <dune/grid/alugrid/3d/alu3dinclude.hh>
00013 #include <dune/grid/alugrid/2d/alu2dinclude.hh>
00014 
00015 namespace Dune
00016 {
00017   template<int dim, int dimw>
00018   class ALUCubeGrid;
00019 
00020   template<int dim, int dimw>
00021   class ALUSimplexGrid;
00022 
00023   template<int dim, int dimw>
00024   class ALUConformGrid;
00025 
00026   template< class GridImp, class GeometryImpl, int nChild >
00027   class ALULocalGeometryStorage 
00028   {
00029     // array with pointers to the geometries
00030     Dune::array< GeometryImpl *, nChild > geoms_;
00031 
00032     // count local geometry creation
00033     int count_;
00034 
00035     // type of grid impl 
00036     typedef typename GridImp :: ctype ctype; 
00037     enum{ dimension       = GridImp :: dimension };
00038     enum{ dimensionworld  = GridImp :: dimensionworld };
00039 
00040     template <int dummy, int dim, int dimworld, int > 
00041     struct CreateGeometries;
00042 
00043     template <int dummy, int dimworld> 
00044     struct CreateGeometries<dummy, 2, dimworld, ALU2DSPACE triangle >
00045     {
00046       template <class Storage> 
00047       static void createGeometries(Storage& storage, 
00048                                    const GeometryType& type,
00049                                    const bool nonConform )
00050       {
00051         if( nonConform ) 
00052         {
00053           typedef ALUGrid< 2, dimworld, simplex, nonconforming, No_Comm > Grid;
00054           storage.template createGeometries< Grid > (type);
00055         }
00056         else 
00057         {
00058           typedef ALUGrid< 2, dimworld, simplex, conforming, No_Comm > Grid;
00059           storage.template createGeometries< Grid > (type);
00060         }
00061       }
00062     };
00063 
00064     template <int dummy> 
00065     struct CreateGeometries<dummy, 3, 3, ALU3DSPACE tetra >
00066     {
00067       template <class Storage> 
00068       static void createGeometries(Storage& storage, 
00069                                    const GeometryType& type,
00070                                    const bool nonConform )
00071       {
00072         assert ( nonConform ) ;
00073         {
00074           typedef ALUGrid< 3, 3, simplex, nonconforming, No_Comm > Grid;
00075           storage.template createGeometries< Grid > (type);
00076         }
00077         /*
00078          // TODO, implement this for refinement of all edges (conforming)
00079         else 
00080         {
00081           typedef ALUGrid< 3, 3, simplex, conforming, No_Comm > Grid;
00082           storage.template createGeometries< Grid > (type);
00083         }
00084         */
00085       }
00086     };
00087 
00088     template <int dummy, int dimworld> 
00089     struct CreateGeometries<dummy, 2, dimworld, ALU2DSPACE quadrilateral >
00090     {
00091       template <class Storage> 
00092       static void createGeometries(Storage& storage, 
00093                                    const GeometryType& type,
00094                                    const bool nonConform )
00095       {
00096         assert ( nonConform ) ;
00097         {
00098           typedef ALUGrid< 2, dimworld, cube, nonconforming, No_Comm > Grid;
00099           storage.template createGeometries< Grid > (type);
00100         }
00101       }
00102     };
00103 
00104     template <int dummy> 
00105     struct CreateGeometries<dummy, 3, 3, ALU3DSPACE hexa >
00106     {
00107       template <class Storage> 
00108       static void createGeometries(Storage& storage, 
00109                                    const GeometryType& type,
00110                                    const bool nonConform )
00111       {
00112         assert( nonConform );
00113         {
00114           typedef ALUGrid< 3, 3, cube, nonconforming, No_Comm > Grid;
00115           storage.template createGeometries< Grid > (type);
00116         }
00117       }
00118     };
00119 
00120   public:
00121     // create empty storage
00122     ALULocalGeometryStorage ( const GeometryType type, const bool nonConform )
00123     : count_( 0 )
00124     {
00125       geoms_.fill( (GeometryImpl *) 0 );
00126       // the idea is to create a grid containing the reference element,
00127       // refine once and the store the father - child relations 
00128       CreateGeometries<0, dimension, dimensionworld, GridImp :: elementType >
00129         ::createGeometries(*this, type, nonConform);
00130     }
00131 
00132     // desctructor deleteing geometries
00133     ~ALULocalGeometryStorage () 
00134     {
00135       for(size_t i=0; i<geoms_.size(); ++i)
00136         if(geoms_[i]) delete geoms_[i];
00137     }
00138 
00139     // check if geometry has been created
00140     bool geomCreated(int child) const { return geoms_[child] != 0; }
00141 
00142     // return reference to local geometry
00143     const GeometryImpl & operator [] (int child) const
00144     {
00145       assert( geomCreated(child) );
00146       return *(geoms_[child]);
00147     }
00148 
00149     template < class Grid >
00150     void createGeometries(const GeometryType& type) 
00151     {
00152       // create factory without verbosity 
00153       GridFactory< Grid > factory( false );
00154 
00155       const Dune::GenericReferenceElement< ctype, dimension > &refElem
00156         = Dune::GenericReferenceElements< ctype, dimension >::general( type );
00157 
00158       // insert vertices 
00159       FieldVector<ctype, dimensionworld> pos( 0 );
00160       const int vxSize = refElem.size(dimension);  
00161       for(int i=0; i<vxSize; ++i)
00162       {
00163         FieldVector<ctype, dimension> position = refElem.position(i, dimension );
00164         // copy position 
00165         for(int d = 0; d<dimension; ++d )
00166           pos[ d ] = position[ d ];
00167 
00168         factory.insertVertex( pos );
00169       }
00170 
00171       std::vector< unsigned int > vertices( vxSize );
00172       // create grid with reference element
00173       for(size_t i=0; i<vertices.size(); ++i) vertices[ i ] = i;
00174       factory.insertElement(type, vertices);
00175 
00176       // save original sbuf
00177       std::streambuf* cerr_sbuf = std::cerr.rdbuf();
00178       std::stringstream tempout; 
00179       // redirect 'cerr' to a 'fout' to avoid unnecessary output in constructors 
00180       std::cerr.rdbuf(tempout.rdbuf()); 
00181 
00182       Grid* gridPtr = factory.createGrid();
00183       Grid& grid    = *gridPtr; 
00184 
00185       // restore the original stream buffer
00186       std::cerr.rdbuf(cerr_sbuf); 
00187 
00188       //std::cerr = savecerr;
00189       
00190       // refine once to get children 
00191       const int level = 1;
00192       grid.globalRefine( level );
00193 
00194       {
00195         typedef typename Grid :: template Partition< All_Partition >:: LevelGridView MacroGridView;
00196         MacroGridView macroView = grid.template levelView< All_Partition > ( 0 );
00197         typedef typename MacroGridView :: template Codim< 0 > :: Iterator Iterator;
00198 
00199         Iterator it = macroView.template begin<0> (); 
00200 
00201         if( it == macroView.template end<0>() ) 
00202           DUNE_THROW(InvalidStateException,"Empty Grid, should contain at least 1 element");
00203 
00204         typedef typename Iterator :: Entity EntityType;
00205 
00206         const EntityType& entity = *it;
00207         const typename EntityType :: Geometry& geo = entity.geometry();
00208         typedef typename EntityType :: HierarchicIterator HierarchicIteratorType;
00209         const HierarchicIteratorType end = entity.hend( level );
00210 
00211         int childNum = 0;
00212         for( HierarchicIteratorType child = entity.hbegin( level ); 
00213              child != end; ++child, ++childNum ) 
00214         {
00215           create( geo, child->geometry(), childNum );
00216         }
00217       }
00218 
00219       // delete grid 
00220       delete gridPtr;
00221     }
00222 
00223   protected:  
00224     // create local geometry
00225     template< class Geometry >
00226     void create ( const Geometry &father, 
00227                   const Geometry &son, 
00228                   const int child )
00229     {
00230       assert( !geomCreated( child ) );
00231       assert( (child >= 0) && (child < nChild) );
00232 
00233       assert( (count_ < nChild) );
00234       ++count_;
00235 
00236       geoms_[ child ] = new GeometryImpl();
00237       geoms_[ child ]->buildGeomInFather( father, son );
00238     }
00239   };  
00240 
00241 } // namespace Dune
00242 
00243 #endif // end HAVE_ALUGRID
00244 
00245 #endif // #ifndef DUNE_ALUGRIDGEOMETRYSTORAGE_HH