dune-grid  2.2.0
dgfyasp.hh
Go to the documentation of this file.
00001 #ifndef DUNE_DGFPARSERYASP_HH
00002 #define DUNE_DGFPARSERYASP_HH
00003 
00004 #include <dune/grid/common/intersection.hh>
00005 #include <dune/grid/yaspgrid.hh>
00006 #include "dgfparser.hh"
00007 
00008 namespace Dune
00009 {
00010   // forward declaration
00011   // -------------------
00012 
00013   template< class GridImp, template< class > class IntersectionImp >
00014   class Intersection;
00015 
00016 
00017   namespace dgf
00018   {
00019 
00033     class YaspGridParameterBlock
00034     : public GridParameterBlock
00035     {
00036     protected:
00037       int _overlap;     // overlap for YaspGrid
00038 
00039     public:
00041       YaspGridParameterBlock( std::istream &in )
00042         : GridParameterBlock( in ),
00043           _overlap( 0 )  // default value
00044       {
00045         // check overlap
00046         if( findtoken( "overlap" ) )
00047         {
00048           int x;
00049           if( getnextentry(x) ) _overlap = x;
00050           else 
00051           {
00052             dwarn << "GridParameterBlock: found keyword `overlap' but no value, defaulting to `" <<  _overlap  <<"' !\n";
00053           }
00054 
00055           if (_overlap < 0) 
00056           {
00057             DUNE_THROW(DGFException,"Negative overlap specified!");
00058           }
00059         }
00060         else
00061         {
00062           dwarn << "YaspGridParameterBlock: Parameter 'overlap' not specified, "
00063                 << "defaulting to '" << _overlap << "'." << std::endl;
00064         }
00065         
00066       }
00067 
00069       int overlap () const
00070       {
00071         return _overlap;
00072       }
00073 
00074     };
00075 
00076   }
00077 
00078   template <int dim>
00079   struct DGFGridFactory< YaspGrid<dim> >
00080   {
00081     typedef YaspGrid<dim> Grid;
00082     const static int dimension = Grid::dimension;
00083     typedef MPIHelper::MPICommunicator MPICommunicatorType;
00084 
00085   private:
00086     typedef FieldVector< double, dimension > Point;
00087     typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
00088 
00089   public:
00090     explicit DGFGridFactory ( std::istream &input,
00091                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
00092     {
00093       generate( input, comm );
00094     }
00095 
00096     explicit DGFGridFactory ( const std::string &filename,
00097                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
00098     {
00099       std::ifstream input( filename.c_str() );
00100       generate( input, comm );
00101     }
00102 
00103     ~DGFGridFactory ()
00104     {
00105       delete boundaryDomainBlock_;
00106     }
00107 
00108     Grid *grid() const
00109     {
00110       return grid_;
00111     }
00112 
00113     template <class Intersection>
00114     bool wasInserted(const Intersection &intersection) const
00115     {
00116       return false;
00117     }
00118     template <class Intersection>
00119     int boundaryId(const Intersection &intersection) const
00120     {
00121       if( boundaryDomainBlock_->isactive() )
00122       {
00123         std::vector< Point > corners;
00124         getCorners( intersection.geometry(), corners );
00125         const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
00126         if( data )
00127           return data->id();
00128         else
00129           return intersection.indexInInside();
00130       }
00131       else
00132         return intersection.indexInInside();
00133     }
00134 
00135     template< int codim >
00136     int numParameters () const
00137     {
00138       return 0;
00139     }
00140 
00141     // return true if boundary parameters found
00142     bool haveBoundaryParameters () const
00143     {
00144       return boundaryDomainBlock_->hasParameter(); 
00145     }
00146 
00147     template < class GG, template < class > class II >
00148     const typename DGFBoundaryParameter::type &
00149       boundaryParameter ( const Intersection< GG, II > & intersection ) const
00150     {
00151       if( haveBoundaryParameters() )
00152       {
00153         std::vector< Point > corners;
00154         getCorners( intersection.geometry(), corners );
00155         const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
00156         if( data )
00157           return data->parameter();
00158         else
00159           return DGFBoundaryParameter::defaultValue();
00160       }
00161       else
00162         return DGFBoundaryParameter::defaultValue();
00163     }
00164 
00165     template< class Entity >
00166     std::vector<double> &parameter ( const Entity &entity )
00167     {
00168       return emptyParam;
00169     }
00170 
00171   private:
00172     void generate( std::istream &gridin, MPICommunicatorType comm );
00173 
00174     template< class Geometry >
00175     static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
00176     {
00177       corners.resize( geometry.corners() );
00178       for( int i = 0; i < geometry.corners(); ++i )
00179       {
00180         const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
00181         for( int j = 0; j < dimension; ++j )
00182           corners[ i ][ j ] = corner[ j ];
00183       }
00184     }
00185 
00186     Grid *grid_;
00187     dgf::BoundaryDomBlock *boundaryDomainBlock_;
00188     std::vector<double> emptyParam;
00189   };
00190 
00191   template< int dim >
00192   inline void DGFGridFactory< YaspGrid< dim > >
00193     ::generate ( std::istream &gridin, MPICommunicatorType comm )
00194   {
00195     dgf::IntervalBlock intervalBlock( gridin );
00196 
00197     if( !intervalBlock.isactive() )
00198       DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
00199 
00200     if( intervalBlock.numIntervals() != 1 )
00201       DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
00202 
00203     if( intervalBlock.dimw() != dim )
00204     {
00205       DUNE_THROW( DGFException,
00206                   "Cannot read an interval of dimension " << intervalBlock.dimw()
00207                   << "into a YaspGrid< " << dim << " >." );
00208     }
00209     
00210     const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
00211 
00212     FieldVector<double,dim> lang;
00213     FieldVector<int,dim>    anz;
00214     for( int i = 0; i < dim; ++i )
00215     {
00216       // check that start point is 0.0
00217       if( fabs( interval.p[ 0 ][ i ] ) > 1e-10 )
00218       {
00219         DUNE_THROW( DGFException,
00220                     "YaspGrid cannot handle grids with non-zero left lower corner." );
00221       }
00222 
00223       lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
00224       anz[ i ]  = interval.n[ i ];
00225     }
00226 
00227     typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
00228     dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
00229     FieldVector< bool, dim > per( false );
00230     const int numTrafos = trafoBlock.numTransformations();
00231     for( int k = 0; k < numTrafos; ++k )
00232     {
00233       const Transformation &trafo = trafoBlock.transformation( k );
00234 
00235       bool identity = true;
00236       for( int i = 0; i < dim; ++i )
00237         for( int j = 0; j < dim; ++j )
00238           identity &= (fabs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
00239       if( !identity )
00240         DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
00241 
00242       int numDirs = 0;
00243       int dir = -1;
00244       for( int i = 0; i < dim; ++i )
00245       {
00246         if( fabs( trafo.shift[ i ] ) < 1e-10 )
00247           continue;
00248         dir = i;
00249         ++numDirs;
00250       }
00251       if( (numDirs != 1) || (fabs( fabs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
00252       {
00253         std::cerr << "Tranformation '" << trafo
00254                   << "' does not map boundaries on boundaries." << std::endl;
00255       }
00256       else
00257         per[ dir ] = true;
00258     }
00259 
00260     // get grid parameters 
00261     dgf::YaspGridParameterBlock grdParam( gridin );
00262 
00263 #if HAVE_MPI
00264     grid_ = new YaspGrid< dim >( comm, lang, anz, per, grdParam.overlap() );
00265 #else
00266     grid_ = new YaspGrid< dim >( lang, anz, per, grdParam.overlap() );
00267 #endif
00268 
00269     boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
00270   }
00271 
00272   template <int dim>
00273   struct DGFGridInfo< YaspGrid<dim> > {
00274       static int refineStepsForHalf() {return 1;}
00275       static double refineWeight() {return pow(0.5,dim);}
00276   };
00277 
00278 
00279 }
00280 #endif // #ifndef DUNE_DGFPARSERYASP_HH