dune-grid  2.2.0
dgfs.hh
Go to the documentation of this file.
00001 #ifndef DUNE_DGFS_HH
00002 #define DUNE_DGFS_HH
00003 
00004 #include <dune/grid/common/intersection.hh>
00005 #include <dune/grid/sgrid.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   template< int dim, int dimworld, class ctype >
00018   struct DGFGridFactory< SGrid< dim, dimworld, ctype > >
00019   {
00020     typedef SGrid< dim, dimworld, ctype > Grid;
00021 
00022     const static int dimension = Grid::dimension;
00023 
00024     typedef MPIHelper::MPICommunicator MPICommunicatorType;
00025 
00026   private:
00027     typedef FieldVector< double, dimension > Point;
00028 
00029     typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
00030 
00031   public:
00032     explicit DGFGridFactory ( std::istream &input,
00033                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
00034     {
00035       generate( input, comm );
00036     }
00037 
00038     explicit DGFGridFactory ( const std::string &filename,
00039                               MPICommunicatorType comm = MPIHelper::getCommunicator() )
00040     {
00041       std::ifstream input( filename.c_str() );
00042       generate( input, comm );
00043     }
00044 
00045     Grid *grid() const
00046     {
00047       return grid_;
00048     }
00049 
00050     template< class Intersection >
00051     bool wasInserted ( const Intersection &intersection ) const
00052     {
00053       return false;
00054     }
00055 
00056     template< class Intersection >
00057     int boundaryId ( const Intersection &intersection ) const
00058     {
00059       if( boundaryDomainBlock_->isactive() )
00060       {
00061         std::vector< Point > corners;
00062         getCorners( intersection.geometry(), corners );
00063         const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
00064         if( data )
00065           return data->id();
00066         else
00067           return intersection.indexInInside();
00068       }
00069       else
00070         return intersection.indexInInside();
00071     }
00072 
00073     template< int codim >
00074     int numParameters () const
00075     {
00076       return 0;
00077     }
00078 
00079     // return true if boundary parameters found
00080     bool haveBoundaryParameters () const
00081     {
00082       return boundaryDomainBlock_->isactive();
00083     }
00084 
00085     template < class GG, template < class > class II >
00086     const typename DGFBoundaryParameter::type &
00087       boundaryParameter ( const Intersection< GG, II > & intersection ) const
00088     {
00089       if( haveBoundaryParameters() )
00090       {
00091         std::vector< Point > corners;
00092         getCorners( intersection.geometry(), corners );
00093         const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
00094         if( data )
00095           return data->parameter();
00096         else
00097           return DGFBoundaryParameter::defaultValue();
00098       }
00099       else
00100         return DGFBoundaryParameter::defaultValue();
00101     }  
00102 
00103 
00104     template< class Entity >
00105     std::vector< double > &parameter ( const Entity &entity )
00106     {
00107       return emptyParam;
00108     }
00109 
00110   private:
00111     void generate( std::istream &gridin, MPICommunicatorType comm );
00112 
00113     template< class Geometry >
00114     static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
00115     {
00116       corners.resize( geometry.corners() );
00117       for( int i = 0; i < geometry.corners(); ++i )
00118       {
00119         const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
00120         for( int j = 0; j < dimension; ++j )
00121           corners[ i ][ j ] = corner[ j ];
00122       }
00123     }
00124 
00125     Grid *grid_;
00126     dgf::BoundaryDomBlock *boundaryDomainBlock_;
00127     std::vector< double > emptyParam;
00128   };
00129 
00130 
00131 
00132   template< int dim, int dimworld, class ctype >
00133   inline void DGFGridFactory< SGrid< dim, dimworld, ctype > >
00134     ::generate ( std::istream &gridin, MPICommunicatorType comm )
00135   {
00136     dgf::IntervalBlock intervalBlock( gridin );
00137 
00138     if( !intervalBlock.isactive() )
00139       DUNE_THROW( DGFException, "SGrid can only be created from an interval block." );
00140 
00141     if( intervalBlock.numIntervals() != 1 )
00142       DUNE_THROW( DGFException, "SGrid can only handle 1 interval block." );
00143 
00144     if( intervalBlock.dimw() != dim )
00145     {
00146       DUNE_THROW( DGFException,
00147                   "Cannot read an interval of dimension " << intervalBlock.dimw()
00148                   << "into a SGrid< " << dim << ", " << dimworld << " >." );
00149     }
00150     
00151     const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
00152 
00153     FieldVector< double, dimension > lower, upper;
00154     FieldVector< int, dimension > anz;
00155     for( int i = 0; i < dimension; ++i )
00156     {
00157       lower[ i ] = interval.p[ 0 ][ i ];
00158       upper[ i ] = interval.p[ 1 ][ i ];
00159       anz[ i ] = interval.n[ i ];
00160     }
00161 
00162     grid_ = new Grid( anz, lower, upper );
00163 
00164     boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
00165   }
00166 
00167 
00168 
00169   template< int dim, int dimworld, class ctype >
00170   struct DGFGridInfo< SGrid< dim, dimworld, ctype > >
00171   {
00172     static int refineStepsForHalf ()
00173     {
00174       return 1;
00175     }
00176 
00177     static double refineWeight ()
00178     {
00179       return 1.0 / double( 1 << dim );
00180     }
00181   };
00182 
00183 }
00184 #endif // #ifndef DUNE_DGFS_HH