dune-grid
2.2.0
|
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> ¶meter ( 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