1 #ifndef DUNE_DGFPARSERYASP_HH
2 #define DUNE_DGFPARSERYASP_HH
6 #include "dgfparser.hh"
13 template<
class Gr
idImp,
template<
class >
class IntersectionImp >
52 dwarn <<
"GridParameterBlock: found keyword `overlap' but no value, defaulting to `" <<
_overlap <<
"' !\n";
62 dwarn <<
"YaspGridParameterBlock: Parameter 'overlap' not specified, "
63 <<
"defaulting to '" <<
_overlap <<
"'." << std::endl;
86 typedef FieldVector< double, dimension > Point;
93 generate( input, comm );
99 std::ifstream input( filename.c_str() );
100 generate( input, comm );
105 delete boundaryDomainBlock_;
113 template <
class Intersection>
118 template <
class Intersection>
121 if( boundaryDomainBlock_->isactive() )
123 std::vector< Point > corners;
124 getCorners( intersection.
geometry(), corners );
125 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
135 template<
int codim >
144 return boundaryDomainBlock_->hasParameter();
147 template <
class GG,
template <
class >
class II >
153 std::vector< Point > corners;
154 getCorners( intersection.
geometry(), corners );
155 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
165 template<
class Entity >
174 template<
class Geometry >
175 static void getCorners (
const Geometry &geometry, std::vector< Point > &corners )
177 corners.resize( geometry.
corners() );
178 for(
int i = 0; i < geometry.
corners(); ++i )
182 corners[ i ][ j ] = corner[ j ];
187 dgf::BoundaryDomBlock *boundaryDomainBlock_;
188 std::vector<double> emptyParam;
192 inline void DGFGridFactory< YaspGrid< dim > >
193 ::generate ( std::istream &gridin, MPICommunicatorType comm )
195 dgf::IntervalBlock intervalBlock( gridin );
197 if( !intervalBlock.isactive() )
198 DUNE_THROW( DGFException,
"YaspGrid can only be created from an interval block." );
200 if( intervalBlock.numIntervals() != 1 )
201 DUNE_THROW( DGFException,
"YaspGrid can only handle 1 interval block." );
203 if( intervalBlock.dimw() != dim )
205 DUNE_THROW( DGFException,
206 "Cannot read an interval of dimension " << intervalBlock.dimw()
207 <<
"into a YaspGrid< " << dim <<
" >." );
210 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
212 FieldVector<double,dim> lang;
213 FieldVector<int,dim> anz;
214 for(
int i = 0; i < dim; ++i )
217 if( fabs( interval.p[ 0 ][ i ] ) > 1e-10 )
219 DUNE_THROW( DGFException,
220 "YaspGrid cannot handle grids with non-zero left lower corner." );
223 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
224 anz[ i ] = interval.n[ i ];
227 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
228 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
229 FieldVector< bool, dim > per(
false );
230 const int numTrafos = trafoBlock.numTransformations();
231 for(
int k = 0; k < numTrafos; ++k )
233 const Transformation &trafo = trafoBlock.transformation( k );
235 bool identity =
true;
236 for(
int i = 0; i < dim; ++i )
237 for(
int j = 0; j < dim; ++j )
238 identity &= (fabs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
240 DUNE_THROW( DGFException,
"YaspGrid can only handle shifts as periodic face transformations." );
244 for(
int i = 0; i < dim; ++i )
246 if( fabs( trafo.shift[ i ] ) < 1e-10 )
251 if( (numDirs != 1) || (fabs( fabs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
253 std::cerr <<
"Tranformation '" << trafo
254 <<
"' does not map boundaries on boundaries." << std::endl;
261 dgf::YaspGridParameterBlock grdParam( gridin );
264 grid_ =
new YaspGrid< dim >( comm, lang, anz, per, grdParam.overlap() );
266 grid_ =
new YaspGrid< dim >( lang, anz, per, grdParam.overlap() );
269 boundaryDomainBlock_ =
new dgf::BoundaryDomBlock( gridin, dimension );
280 #endif // #ifndef DUNE_DGFPARSERYASP_HH