5#ifndef DUNE_DGFPARSERYASP_HH
6#define DUNE_DGFPARSERYASP_HH
8#include <dune/grid/common/intersection.hh>
10#include "dgfparser.hh"
18 template<
class Gr
idImp,
class IntersectionImp >
51 if( findtoken(
"overlap" ) )
54 if( getnextentry(x) ) _overlap = x;
57 dwarn <<
"GridParameterBlock: found keyword `overlap' but no value, defaulting to `" << _overlap <<
"' !\n";
67 dwarn <<
"YaspGridParameterBlock: Parameter 'overlap' not specified, "
68 <<
"defaulting to '" << _overlap <<
"'." << std::endl;
86 template <
typename ctype,
int dim>
95 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
98 explicit DGFGridFactory ( std::istream &input,
101 generate( input, comm );
104 explicit DGFGridFactory (
const std::string &filename,
107 std::ifstream input( filename.c_str() );
110 generate( input, comm );
115 delete boundaryDomainBlock_;
123 template <
class Intersection>
124 bool wasInserted(
const Intersection &intersection)
const
129 template <
class Intersection>
132 if( boundaryDomainBlock_->isactive() )
134 std::vector< Point > corners;
135 getCorners( intersection.
geometry(), corners );
136 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
146 template<
int codim >
147 int numParameters ()
const
153 bool haveBoundaryParameters ()
const
155 return boundaryDomainBlock_->hasParameter();
158 template<
class GG,
class II >
162 if( haveBoundaryParameters() )
164 std::vector< Point > corners;
165 getCorners( intersection.
geometry(), corners );
166 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
168 return data->parameter();
176 template<
class Entity >
177 std::vector<double> ¶meter (
const Entity & )
183 void generate( std::istream &gridin, MPICommunicatorType comm );
185 template<
class Geometry >
186 static void getCorners (
const Geometry &geometry, std::vector< Point > &corners )
188 corners.resize( geometry.
corners() );
189 for(
int i = 0; i < geometry.
corners(); ++i )
192 for(
int j = 0; j < dimension; ++j )
193 corners[ i ][ j ] = corner[ j ];
198 dgf::BoundaryDomBlock *boundaryDomainBlock_;
199 std::vector<double> emptyParam;
203 template<
typename ctype,
int dim >
204 inline void DGFGridFactory< YaspGrid< dim , EquidistantCoordinates<ctype, dim> > >
205 ::generate ( std::istream &gridin, MPICommunicatorType comm )
208 dgf::IntervalBlock intervalBlock( gridin );
210 if( !intervalBlock.isactive() )
213 if( intervalBlock.numIntervals() != 1 )
216 if( intervalBlock.dimw() != dim )
219 "Cannot read an interval of dimension " << intervalBlock.dimw()
220 <<
" into a YaspGrid< " << dim <<
" >." );
223 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
225 FieldVector<ctype, dim> lang;
226 std::array<int,dim> anz;
227 for(
int i = 0; i < dim; ++i )
230 if( abs( interval.p[ 0 ][ i ] ) > 1e-10 )
233 "YaspGrid cannot handle grids with non-zero left lower corner." );
236 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
237 anz[ i ] = interval.n[ i ];
240 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
241 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
242 std::bitset< dim > per;
243 const int numTrafos = trafoBlock.numTransformations();
244 for(
int k = 0; k < numTrafos; ++k )
246 const Transformation &trafo = trafoBlock.transformation( k );
248 bool identity =
true;
249 for(
int i = 0; i < dim; ++i )
250 for(
int j = 0; j < dim; ++j )
251 identity &= (abs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
253 DUNE_THROW( DGFException,
"YaspGrid can only handle shifts as periodic face transformations." );
257 for(
int i = 0; i < dim; ++i )
259 if( abs( trafo.shift[ i ] ) < 1e-10 )
264 if( (numDirs != 1) || (abs( abs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
266 std::cerr <<
"Transformation '" << trafo
267 <<
"' does not map boundaries on boundaries." << std::endl;
274 dgf::YaspGridParameterBlock grdParam( gridin );
276 grid_ =
new YaspGrid< dim , EquidistantCoordinates<ctype, dim> >( lang, anz, per, grdParam.overlap(), comm );
278 boundaryDomainBlock_ =
new dgf::BoundaryDomBlock( gridin, dimension );
284 template <
typename ctype,
int dim>
293 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
296 explicit DGFGridFactory ( std::istream &input,
299 generate( input, comm );
302 explicit DGFGridFactory (
const std::string &filename,
305 std::ifstream input( filename.c_str() );
306 generate( input, comm );
311 delete boundaryDomainBlock_;
319 template <
class Intersection>
320 bool wasInserted(
const Intersection &intersection)
const
325 template <
class Intersection>
328 if( boundaryDomainBlock_->isactive() )
330 std::vector< Point > corners;
331 getCorners( intersection.
geometry(), corners );
332 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
342 template<
int codim >
343 int numParameters ()
const
349 bool haveBoundaryParameters ()
const
351 return boundaryDomainBlock_->hasParameter();
354 template<
class GG,
class II >
358 if( haveBoundaryParameters() )
360 std::vector< Point > corners;
361 getCorners( intersection.
geometry(), corners );
362 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
364 return data->parameter();
372 template<
class Entity >
373 std::vector<double> ¶meter ( [[maybe_unused]]
const Entity &entity )
379 void generate( std::istream &gridin, MPICommunicatorType comm );
381 template<
class Geometry >
382 static void getCorners (
const Geometry &geometry, std::vector< Point > &corners )
384 corners.resize( geometry.
corners() );
385 for(
int i = 0; i < geometry.
corners(); ++i )
388 for(
int j = 0; j < dimension; ++j )
389 corners[ i ][ j ] = corner[ j ];
394 dgf::BoundaryDomBlock *boundaryDomainBlock_;
395 std::vector<double> emptyParam;
399 template<
typename ctype,
int dim >
400 inline void DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > >
401 ::generate ( std::istream &gridin, MPICommunicatorType comm )
404 dgf::IntervalBlock intervalBlock( gridin );
406 if( !intervalBlock.isactive() )
409 if( intervalBlock.numIntervals() != 1 )
412 if( intervalBlock.dimw() != dim )
415 "Cannot read an interval of dimension "
416 << intervalBlock.dimw()
417 <<
" into a YaspGrid< " << dim <<
" >." );
420 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
422 FieldVector<ctype, dim> lower;
423 FieldVector<ctype, dim> upper;
424 std::array<int,dim> anz;
425 for(
int i = 0; i < dim; ++i )
427 lower[ i ] = interval.p[ 0 ][ i ];
428 upper[ i ] = interval.p[ 1 ][ i ];
429 anz[ i ] = interval.n[ i ];
432 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
433 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
435 const int numTrafos = trafoBlock.numTransformations();
436 for(
int k = 0; k < numTrafos; ++k )
438 const Transformation &trafo = trafoBlock.transformation( k );
440 bool identity =
true;
441 for(
int i = 0; i < dim; ++i )
442 for(
int j = 0; j < dim; ++j )
443 identity &= (abs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
445 DUNE_THROW( DGFException,
"YaspGrid can only handle shifts as periodic face transformations." );
449 for(
int currentDir = 0; currentDir < dim; ++currentDir )
451 if( abs( trafo.shift[ currentDir ] ) > 1e-10 )
458 || (abs( abs( trafo.shift[ dir ] ) - abs( upper[ dir ] - lower[ dir ] ) ) >= 1e-10) )
460 std::cerr <<
"Transformation '" << trafo
461 <<
"' does not map boundaries on boundaries." << std::endl;
470 dgf::YaspGridParameterBlock grdParam( gridin );
472 grid_ =
new YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >
473 ( lower, upper, anz,
periodic, grdParam.overlap(), comm );
475 boundaryDomainBlock_ =
new dgf::BoundaryDomBlock( gridin, dimension );
483 template<
class ctype,
int dim >
484 class DGFGridFactory<
Dune::
YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > >
488 template<
typename In >
489 DGFGridFactory (
const In & ) {}
492 throw std::invalid_argument(
"Tensor product coordinates for YaspGrid are currently not supported via the DGFGridFactory" );
496 template <
typename Coordinates,
int dim>
499 static double refineWeight() {
return std::pow(0.5,dim);}
exception class for IO errors in the DGF parser
Definition: dgfexception.hh:16
Wrapper class for entities.
Definition: entity.hh:66
Container for equidistant coordinates in a YaspGrid.
Definition: coordinates.hh:29
Container for equidistant coordinates in a YaspGrid with non-trivial origin.
Definition: coordinates.hh:131
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Wrapper class for geometries.
Definition: geometry.hh:71
GlobalCoordinate corner(int i) const
Obtain a corner of the geometry.
Definition: geometry.hh:219
int corners() const
Return the number of corners of the reference element.
Definition: geometry.hh:205
static constexpr int dimension
The dimension of the grid.
Definition: grid.hh:387
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: intersection.hh:164
Geometry geometry() const
geometrical information about the intersection in global coordinates.
Definition: intersection.hh:323
int indexInInside() const
Local index of codim 1 entity in the inside() entity where intersection is contained in.
Definition: intersection.hh:346
MPI_Comm MPICommunicator
The type of the mpi communicator.
Definition: mpihelper.hh:192
static MPICommunicator getCommunicator()
get the default communicator
Definition: mpihelper.hh:200
[ provides Dune::Grid ]
Definition: yaspgrid.hh:166
Common Grid parameters.
Definition: gridparameter.hh:35
Grid parameters for YaspGrid.
Definition: dgfyasp.hh:40
YaspGridParameterBlock(std::istream &in)
constructor taking istream
Definition: dgfyasp.hh:46
int overlap() const
get dimension of world found in block
Definition: dgfyasp.hh:74
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
auto periodic(RawPreBasisIndicator &&rawPreBasisIndicator, PIS &&periodicIndexSet)
Create a pre-basis factory that can create a periodic pre-basis.
Definition: periodicbasis.hh:203
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:162
Dune namespace.
Definition: alignedallocator.hh:13
static const type & defaultValue()
default constructor
Definition: parser.hh:28
std::string type
type of additional boundary parameters
Definition: parser.hh:25
Some simple static information for a given GridType.
Definition: dgfparser.hh:56
static double refineWeight()
static int refineStepsForHalf()
number of globalRefine steps needed to refuce h by 0.5