3#ifndef DUNE_DGFPARSERYASP_HH
4#define DUNE_DGFPARSERYASP_HH
6#include <dune/grid/common/intersection.hh>
16 template<
class Gr
idImp,
class IntersectionImp >
49 if( findtoken(
"overlap" ) )
52 if( getnextentry(x) ) _overlap = x;
55 dwarn <<
"GridParameterBlock: found keyword `overlap' but no value, defaulting to `" << _overlap <<
"' !\n";
65 dwarn <<
"YaspGridParameterBlock: Parameter 'overlap' not specified, "
66 <<
"defaulting to '" << _overlap <<
"'." << std::endl;
84 template <
typename ctype,
int dim>
93 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
96 explicit DGFGridFactory ( std::istream &input,
99 generate( input, comm );
102 explicit DGFGridFactory (
const std::string &filename,
105 std::ifstream input( filename.c_str() );
108 generate( input, comm );
113 delete boundaryDomainBlock_;
121 template <
class Intersection>
122 bool wasInserted(
const Intersection &intersection)
const
127 template <
class Intersection>
130 if( boundaryDomainBlock_->isactive() )
132 std::vector< Point > corners;
133 getCorners( intersection.
geometry(), corners );
134 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
144 template<
int codim >
145 int numParameters ()
const
151 bool haveBoundaryParameters ()
const
153 return boundaryDomainBlock_->hasParameter();
156 template<
class GG,
class II >
160 if( haveBoundaryParameters() )
162 std::vector< Point > corners;
163 getCorners( intersection.
geometry(), corners );
164 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
166 return data->parameter();
174 template<
class Entity >
175 std::vector<double> ¶meter (
const Entity & )
181 void generate( std::istream &gridin, MPICommunicatorType comm );
183 template<
class Geometry >
184 static void getCorners (
const Geometry &geometry, std::vector< Point > &corners )
186 corners.resize( geometry.
corners() );
187 for(
int i = 0; i < geometry.
corners(); ++i )
190 for(
int j = 0; j < dimension; ++j )
191 corners[ i ][ j ] = corner[ j ];
196 dgf::BoundaryDomBlock *boundaryDomainBlock_;
197 std::vector<double> emptyParam;
201 template<
typename ctype,
int dim >
202 inline void DGFGridFactory< YaspGrid< dim , EquidistantCoordinates<ctype, dim> > >
203 ::generate ( std::istream &gridin, MPICommunicatorType comm )
206 dgf::IntervalBlock intervalBlock( gridin );
208 if( !intervalBlock.isactive() )
211 if( intervalBlock.numIntervals() != 1 )
214 if( intervalBlock.dimw() != dim )
217 "Cannot read an interval of dimension " << intervalBlock.dimw()
218 <<
" into a YaspGrid< " << dim <<
" >." );
221 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
223 FieldVector<ctype, dim> lang;
224 std::array<int,dim> anz;
225 for(
int i = 0; i < dim; ++i )
228 if( abs( interval.p[ 0 ][ i ] ) > 1e-10 )
231 "YaspGrid cannot handle grids with non-zero left lower corner." );
234 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
235 anz[ i ] = interval.n[ i ];
238 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
239 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
240 std::bitset< dim > per;
241 const int numTrafos = trafoBlock.numTransformations();
242 for(
int k = 0; k < numTrafos; ++k )
244 const Transformation &trafo = trafoBlock.transformation( k );
246 bool identity =
true;
247 for(
int i = 0; i < dim; ++i )
248 for(
int j = 0; j < dim; ++j )
249 identity &= (abs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
251 DUNE_THROW( DGFException,
"YaspGrid can only handle shifts as periodic face transformations." );
255 for(
int i = 0; i < dim; ++i )
257 if( abs( trafo.shift[ i ] ) < 1e-10 )
262 if( (numDirs != 1) || (abs( abs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
264 std::cerr <<
"Tranformation '" << trafo
265 <<
"' does not map boundaries on boundaries." << std::endl;
272 dgf::YaspGridParameterBlock grdParam( gridin );
274 grid_ =
new YaspGrid< dim , EquidistantCoordinates<ctype, dim> >( lang, anz, per, grdParam.overlap(), comm );
276 boundaryDomainBlock_ =
new dgf::BoundaryDomBlock( gridin, dimension );
279 template <
typename ctype,
int dim>
280 struct DGFGridInfo< YaspGrid<dim , EquidistantCoordinates<ctype, dim> > > {
282 static double refineWeight() {
return std::pow(0.5,dim);}
288 template <
typename ctype,
int dim>
297 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
300 explicit DGFGridFactory ( std::istream &input,
303 generate( input, comm );
306 explicit DGFGridFactory (
const std::string &filename,
309 std::ifstream input( filename.c_str() );
310 generate( input, comm );
315 delete boundaryDomainBlock_;
323 template <
class Intersection>
324 bool wasInserted(
const Intersection &intersection)
const
329 template <
class Intersection>
332 if( boundaryDomainBlock_->isactive() )
334 std::vector< Point > corners;
335 getCorners( intersection.
geometry(), corners );
336 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
346 template<
int codim >
347 int numParameters ()
const
353 bool haveBoundaryParameters ()
const
355 return boundaryDomainBlock_->hasParameter();
358 template<
class GG,
class II >
362 if( haveBoundaryParameters() )
364 std::vector< Point > corners;
365 getCorners( intersection.
geometry(), corners );
366 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
368 return data->parameter();
376 template<
class Entity >
377 std::vector<double> ¶meter ( [[maybe_unused]]
const Entity &entity )
383 void generate( std::istream &gridin, MPICommunicatorType comm );
385 template<
class Geometry >
386 static void getCorners (
const Geometry &geometry, std::vector< Point > &corners )
388 corners.resize( geometry.
corners() );
389 for(
int i = 0; i < geometry.
corners(); ++i )
392 for(
int j = 0; j < dimension; ++j )
393 corners[ i ][ j ] = corner[ j ];
398 dgf::BoundaryDomBlock *boundaryDomainBlock_;
399 std::vector<double> emptyParam;
403 template<
typename ctype,
int dim >
404 inline void DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > >
405 ::generate ( std::istream &gridin, MPICommunicatorType comm )
408 dgf::IntervalBlock intervalBlock( gridin );
410 if( !intervalBlock.isactive() )
413 if( intervalBlock.numIntervals() != 1 )
416 if( intervalBlock.dimw() != dim )
419 "Cannot read an interval of dimension "
420 << intervalBlock.dimw()
421 <<
" into a YaspGrid< " << dim <<
" >." );
424 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
426 FieldVector<ctype, dim> lower;
427 FieldVector<ctype, dim> upper;
428 std::array<int,dim> anz;
429 for(
int i = 0; i < dim; ++i )
431 lower[ i ] = interval.p[ 0 ][ i ];
432 upper[ i ] = interval.p[ 1 ][ i ];
433 anz[ i ] = interval.n[ i ];
436 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
437 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
439 const int numTrafos = trafoBlock.numTransformations();
440 for(
int k = 0; k < numTrafos; ++k )
442 const Transformation &trafo = trafoBlock.transformation( k );
444 bool identity =
true;
445 for(
int i = 0; i < dim; ++i )
446 for(
int j = 0; j < dim; ++j )
447 identity &= (abs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
449 DUNE_THROW( DGFException,
"YaspGrid can only handle shifts as periodic face transformations." );
453 for(
int currentDir = 0; currentDir < dim; ++currentDir )
455 if( abs( trafo.shift[ currentDir ] ) > 1e-10 )
462 || (abs( abs( trafo.shift[ dir ] ) - abs( upper[ dir ] - lower[ dir ] ) ) >= 1e-10) )
464 std::cerr <<
"Tranformation '" << trafo
465 <<
"' does not map boundaries on boundaries." << std::endl;
474 dgf::YaspGridParameterBlock grdParam( gridin );
476 grid_ =
new YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >
477 ( lower, upper, anz,
periodic, grdParam.overlap(), comm );
479 boundaryDomainBlock_ =
new dgf::BoundaryDomBlock( gridin, dimension );
482 template <
typename ctype,
int dim>
483 struct DGFGridInfo< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > > {
485 static double refineWeight() {
return std::pow(0.5,dim);}
exception class for IO errors in the DGF parser
Definition: dgfexception.hh:14
Wrapper class for entities.
Definition: entity.hh:64
Container for equidistant coordinates in a YaspGrid.
Definition: coordinates.hh:27
Container for equidistant coordinates in a YaspGrid with non-trivial origin.
Definition: coordinates.hh:129
vector space out of a tensor product of fields.
Definition: fvector.hh:95
Wrapper class for geometries.
Definition: geometry.hh:67
GlobalCoordinate corner(int i) const
Obtain a corner of the geometry.
Definition: geometry.hh:156
int corners() const
Return the number of corners of the reference element.
Definition: geometry.hh:142
@ dimension
The dimension of the grid.
Definition: grid.hh:386
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: intersection.hh:162
Geometry geometry() const
geometrical information about the intersection in global coordinates.
Definition: intersection.hh:321
int indexInInside() const
Local index of codim 1 entity in the inside() entity where intersection is contained in.
Definition: intersection.hh:344
MPI_Comm MPICommunicator
The type of the mpi communicator.
Definition: mpihelper.hh:188
static MPICommunicator getCommunicator()
get the default communicator
Definition: mpihelper.hh:196
[ provides Dune::Grid ]
Definition: yaspgrid.hh:160
Common Grid parameters.
Definition: gridparameter.hh:33
Grid parameters for YaspGrid.
Definition: dgfyasp.hh:38
YaspGridParameterBlock(std::istream &in)
constructor taking istream
Definition: dgfyasp.hh:44
int overlap() const
get dimension of world found in block
Definition: dgfyasp.hh:72
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
auto periodic(RawPreBasisIndicator &&rawPreBasisIndicator, PIS &&periodicIndexSet)
Create a pre-basis factory that can create a periodic pre-basis.
Definition: periodicbasis.hh:195
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:159
Dune namespace.
Definition: alignedallocator.hh:11
static const type & defaultValue()
default constructor
Definition: parser.hh:26
std::string type
type of additional boundary parameters
Definition: parser.hh:23
static double refineWeight()
static int refineStepsForHalf()
number of globalRefine steps needed to refuce h by 0.5