dune-grid 2.9.0
dgfyasp.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright (C) DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_DGFPARSERYASP_HH
6#define DUNE_DGFPARSERYASP_HH
7
10#include "dgfparser.hh"
11
12namespace Dune
13{
14
15 // External Forward Declarations
16 // -----------------------------
17
18 template< class GridImp, class IntersectionImp >
19 class Intersection;
20
21
22 namespace dgf
23 {
24
39 : public GridParameterBlock
40 {
41 protected:
42 int _overlap; // overlap for YaspGrid
43
44 public:
46 YaspGridParameterBlock( std::istream &in )
47 : GridParameterBlock( in ),
48 _overlap( 0 ) // default value
49 {
50 // check overlap
51 if( findtoken( "overlap" ) )
52 {
53 int x;
54 if( getnextentry(x) ) _overlap = x;
55 else
56 {
57 dwarn << "GridParameterBlock: found keyword `overlap' but no value, defaulting to `" << _overlap <<"' !\n";
58 }
59
60 if (_overlap < 0)
61 {
62 DUNE_THROW(DGFException,"Negative overlap specified!");
63 }
64 }
65 else
66 {
67 dwarn << "YaspGridParameterBlock: Parameter 'overlap' not specified, "
68 << "defaulting to '" << _overlap << "'." << std::endl;
69 }
70
71 }
72
74 int overlap () const
75 {
76 return _overlap;
77 }
78
79 };
80
81 }
82
86 template <typename ctype, int dim>
87 struct DGFGridFactory< YaspGrid<dim, EquidistantCoordinates<ctype, dim> > >
88 {
90 const static int dimension = Grid::dimension;
91 typedef MPIHelper::MPICommunicator MPICommunicatorType;
92
93 private:
94 typedef FieldVector< ctype, dimension > Point;
96
97 public:
98 explicit DGFGridFactory ( std::istream &input,
99 MPICommunicatorType comm = MPIHelper::getCommunicator() )
100 {
101 generate( input, comm );
102 }
103
104 explicit DGFGridFactory ( const std::string &filename,
105 MPICommunicatorType comm = MPIHelper::getCommunicator() )
106 {
107 std::ifstream input( filename.c_str() );
108 if( !input )
109 DUNE_THROW( DGFException, "Error: Macrofile '" << filename << "' not found" );
110 generate( input, comm );
111 }
112
114 {
115 delete boundaryDomainBlock_;
116 }
117
118 Grid *grid() const
119 {
120 return grid_;
121 }
122
123 template <class Intersection>
124 bool wasInserted(const Intersection &intersection) const
125 {
126 return false;
127 }
128
129 template <class Intersection>
130 int boundaryId(const Intersection &intersection) const
131 {
132 if( boundaryDomainBlock_->isactive() )
133 {
134 std::vector< Point > corners;
135 getCorners( intersection.geometry(), corners );
136 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
137 if( data )
138 return data->id();
139 else
140 return intersection.indexInInside();
141 }
142 else
143 return intersection.indexInInside();
144 }
145
146 template< int codim >
147 int numParameters () const
148 {
149 return 0;
150 }
151
152 // return true if boundary parameters found
154 {
155 return boundaryDomainBlock_->hasParameter();
156 }
157
158 template< class GG, class II >
159 const typename DGFBoundaryParameter::type &
160 boundaryParameter ( const Intersection< GG, II > & intersection ) const
161 {
163 {
164 std::vector< Point > corners;
165 getCorners( intersection.geometry(), corners );
166 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
167 if( data )
168 return data->parameter();
169 else
171 }
172 else
174 }
175
176 template< class Entity >
177 std::vector<double> &parameter ( const Entity & )
178 {
179 return emptyParam;
180 }
181
182 private:
183 void generate( std::istream &gridin, MPICommunicatorType comm );
184
185 template< class Geometry >
186 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
187 {
188 corners.resize( geometry.corners() );
189 for( int i = 0; i < geometry.corners(); ++i )
190 {
191 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
192 for( int j = 0; j < dimension; ++j )
193 corners[ i ][ j ] = corner[ j ];
194 }
195 }
196
197 Grid *grid_;
198 dgf::BoundaryDomBlock *boundaryDomainBlock_;
199 std::vector<double> emptyParam;
200 };
201
202 // generate YaspGrid from the provided DGF
203 template< typename ctype, int dim >
204 inline void DGFGridFactory< YaspGrid< dim , EquidistantCoordinates<ctype, dim> > >
205 ::generate ( std::istream &gridin, MPICommunicatorType comm )
206 {
207 using std::abs;
208 dgf::IntervalBlock intervalBlock( gridin );
209
210 if( !intervalBlock.isactive() )
211 DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
212
213 if( intervalBlock.numIntervals() != 1 )
214 DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
215
216 if( intervalBlock.dimw() != dim )
217 {
218 DUNE_THROW( DGFException,
219 "Cannot read an interval of dimension " << intervalBlock.dimw()
220 << " into a YaspGrid< " << dim << " >." );
221 }
222
223 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
224
225 FieldVector<ctype, dim> lang;
226 std::array<int,dim> anz;
227 for( int i = 0; i < dim; ++i )
228 {
229 // check that start point is 0.0
230 if( abs( interval.p[ 0 ][ i ] ) > 1e-10 )
231 {
232 DUNE_THROW( DGFException,
233 "YaspGrid cannot handle grids with non-zero left lower corner." );
234 }
235
236 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
237 anz[ i ] = interval.n[ i ];
238 }
239
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 )
245 {
246 const Transformation &trafo = trafoBlock.transformation( k );
247
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);
252 if( !identity )
253 DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
254
255 int numDirs = 0;
256 int dir = -1;
257 for( int i = 0; i < dim; ++i )
258 {
259 if( abs( trafo.shift[ i ] ) < 1e-10 )
260 continue;
261 dir = i;
262 ++numDirs;
263 }
264 if( (numDirs != 1) || (abs( abs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
265 {
266 std::cerr << "Tranformation '" << trafo
267 << "' does not map boundaries on boundaries." << std::endl;
268 }
269 else
270 per[ dir ] = true;
271 }
272
273 // get grid parameters
274 dgf::YaspGridParameterBlock grdParam( gridin );
275
276 grid_ = new YaspGrid< dim , EquidistantCoordinates<ctype, dim> >( lang, anz, per, grdParam.overlap(), comm );
277
278 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
279 }
280
284 template <typename ctype, int dim>
286 {
288 const static int dimension = Grid::dimension;
289 typedef MPIHelper::MPICommunicator MPICommunicatorType;
290
291 private:
292 typedef FieldVector< ctype, dimension > Point;
294
295 public:
296 explicit DGFGridFactory ( std::istream &input,
297 MPICommunicatorType comm = MPIHelper::getCommunicator() )
298 {
299 generate( input, comm );
300 }
301
302 explicit DGFGridFactory ( const std::string &filename,
303 MPICommunicatorType comm = MPIHelper::getCommunicator() )
304 {
305 std::ifstream input( filename.c_str() );
306 generate( input, comm );
307 }
308
310 {
311 delete boundaryDomainBlock_;
312 }
313
314 Grid *grid() const
315 {
316 return grid_;
317 }
318
319 template <class Intersection>
320 bool wasInserted(const Intersection &intersection) const
321 {
322 return false;
323 }
324
325 template <class Intersection>
326 int boundaryId(const Intersection &intersection) const
327 {
328 if( boundaryDomainBlock_->isactive() )
329 {
330 std::vector< Point > corners;
331 getCorners( intersection.geometry(), corners );
332 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
333 if( data )
334 return data->id();
335 else
336 return intersection.indexInInside();
337 }
338 else
339 return intersection.indexInInside();
340 }
341
342 template< int codim >
343 int numParameters () const
344 {
345 return 0;
346 }
347
348 // return true if boundary parameters found
350 {
351 return boundaryDomainBlock_->hasParameter();
352 }
353
354 template< class GG, class II >
355 const typename DGFBoundaryParameter::type &
356 boundaryParameter ( const Intersection< GG, II > & intersection ) const
357 {
359 {
360 std::vector< Point > corners;
361 getCorners( intersection.geometry(), corners );
362 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
363 if( data )
364 return data->parameter();
365 else
367 }
368 else
370 }
371
372 template< class Entity >
373 std::vector<double> &parameter ( [[maybe_unused]] const Entity &entity )
374 {
375 return emptyParam;
376 }
377
378 private:
379 void generate( std::istream &gridin, MPICommunicatorType comm );
380
381 template< class Geometry >
382 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
383 {
384 corners.resize( geometry.corners() );
385 for( int i = 0; i < geometry.corners(); ++i )
386 {
387 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
388 for( int j = 0; j < dimension; ++j )
389 corners[ i ][ j ] = corner[ j ];
390 }
391 }
392
393 Grid *grid_;
394 dgf::BoundaryDomBlock *boundaryDomainBlock_;
395 std::vector<double> emptyParam;
396 };
397
398 // generate YaspGrid from the provided DGF
399 template< typename ctype, int dim >
400 inline void DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > >
401 ::generate ( std::istream &gridin, MPICommunicatorType comm )
402 {
403 using std::abs;
404 dgf::IntervalBlock intervalBlock( gridin );
405
406 if( !intervalBlock.isactive() )
407 DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
408
409 if( intervalBlock.numIntervals() != 1 )
410 DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
411
412 if( intervalBlock.dimw() != dim )
413 {
414 DUNE_THROW( DGFException,
415 "Cannot read an interval of dimension "
416 << intervalBlock.dimw()
417 << " into a YaspGrid< " << dim << " >." );
418 }
419
420 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
421
422 FieldVector<ctype, dim> lower;
423 FieldVector<ctype, dim> upper;
424 std::array<int,dim> anz;
425 for( int i = 0; i < dim; ++i )
426 {
427 lower[ i ] = interval.p[ 0 ][ i ];
428 upper[ i ] = interval.p[ 1 ][ i ];
429 anz[ i ] = interval.n[ i ];
430 }
431
433 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
434 std::bitset< dim > periodic;
435 const int numTrafos = trafoBlock.numTransformations();
436 for( int k = 0; k < numTrafos; ++k )
437 {
438 const Transformation &trafo = trafoBlock.transformation( k );
439
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);
444 if( !identity )
445 DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
446
447 int numDirs = 0;
448 int dir = -1;
449 for( int currentDir = 0; currentDir < dim; ++currentDir )
450 {
451 if( abs( trafo.shift[ currentDir ] ) > 1e-10 )
452 {
453 dir = currentDir;
454 ++numDirs;
455 }
456 }
457 if ( (numDirs != 1)
458 || (abs( abs( trafo.shift[ dir ] ) - abs( upper[ dir ] - lower[ dir ] ) ) >= 1e-10) )
459 {
460 std::cerr << "Tranformation '" << trafo
461 << "' does not map boundaries on boundaries." << std::endl;
462 }
463 else
464 {
465 periodic[ dir ] = true;
466 }
467 }
468
469 // get grid parameters
470 dgf::YaspGridParameterBlock grdParam( gridin );
471
472 grid_ = new YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >
473 ( lower, upper, anz, periodic, grdParam.overlap(), comm );
474
475 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
476 }
477
483 template< class ctype, int dim >
484 class DGFGridFactory< Dune::YaspGrid<dim, Dune::TensorProductCoordinates<ctype, dim> > >
485 {
487 public:
488 template< typename In >
489 DGFGridFactory ( const In & ) {}
491 {
492 throw std::invalid_argument( "Tensor product coordinates for YaspGrid are currently not supported via the DGFGridFactory" );
493 }
494 };
495
496 template <typename Coordinates, int dim>
497 struct DGFGridInfo< YaspGrid<dim , Coordinates> > {
498 static int refineStepsForHalf() {return 1;}
499 static double refineWeight() {return std::pow(0.5,dim);}
500 };
501
502}
503#endif // #ifndef DUNE_DGFPARSERYASP_HH
Include standard header files.
Definition: agrid.hh:60
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:328
ALBERTA AFF_TRAFO AffineTransformation
Definition: misc.hh:52
Definition: dgfgridfactory.hh:38
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgfgridfactory.hh:41
G Grid
Definition: dgfgridfactory.hh:39
static const int dimension
Definition: dgfgridfactory.hh:40
bool haveBoundaryParameters() const
Definition: dgfgridfactory.hh:151
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: common/intersection.hh:164
Geometry geometry() const
geometrical information about the intersection in global coordinates.
Definition: common/intersection.hh:323
int indexInInside() const
Local index of codim 1 entity in the inside() entity where intersection is contained in.
Definition: common/intersection.hh:346
Wrapper class for entities.
Definition: common/entity.hh:66
Wrapper class for geometries.
Definition: common/geometry.hh:71
FieldVector< ctype, cdim > GlobalCoordinate
type of the global coordinates
Definition: common/geometry.hh:106
GlobalCoordinate corner(int i) const
Obtain a corner of the geometry.
Definition: common/geometry.hh:219
int corners() const
Return the number of corners of the reference element.
Definition: common/geometry.hh:205
static constexpr int dimension
The dimension of the grid.
Definition: common/grid.hh:387
bool findtoken(std ::string token)
Definition: basic.cc:123
bool getnextentry(ENTRY &entry)
Definition: basic.hh:63
Definition: boundarydom.hh:23
int id() const
Definition: boundarydom.hh:42
const BoundaryParameter & parameter() const
Definition: boundarydom.hh:54
Definition: boundarydom.hh:195
Common Grid parameters.
Definition: gridparameter.hh:35
exception class for IO errors in the DGF parser
Definition: dgfexception.hh:16
Some simple static information for a given GridType.
Definition: io/file/dgfparser/dgfparser.hh:56
Grid parameters for YaspGrid.
Definition: dgfyasp.hh:40
int _overlap
Definition: dgfyasp.hh:42
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
DGFGridFactory(std::istream &input, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgfyasp.hh:98
DGFGridFactory(const std::string &filename, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgfyasp.hh:104
bool wasInserted(const Intersection &intersection) const
Definition: dgfyasp.hh:124
std::vector< double > & parameter(const Entity &)
Definition: dgfyasp.hh:177
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgfyasp.hh:91
const DGFBoundaryParameter::type & boundaryParameter(const Intersection< GG, II > &intersection) const
Definition: dgfyasp.hh:160
YaspGrid< dim, EquidistantCoordinates< ctype, dim > > Grid
Definition: dgfyasp.hh:89
int boundaryId(const Intersection &intersection) const
Definition: dgfyasp.hh:130
const DGFBoundaryParameter::type & boundaryParameter(const Intersection< GG, II > &intersection) const
Definition: dgfyasp.hh:356
DGFGridFactory(const std::string &filename, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgfyasp.hh:302
int boundaryId(const Intersection &intersection) const
Definition: dgfyasp.hh:326
YaspGrid< dim, EquidistantOffsetCoordinates< ctype, dim > > Grid
Definition: dgfyasp.hh:287
bool wasInserted(const Intersection &intersection) const
Definition: dgfyasp.hh:320
std::vector< double > & parameter(const Entity &entity)
Definition: dgfyasp.hh:373
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgfyasp.hh:289
DGFGridFactory(std::istream &input, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgfyasp.hh:296
static double refineWeight()
Definition: dgfyasp.hh:499
static int refineStepsForHalf()
Definition: dgfyasp.hh:498
static const type & defaultValue()
default constructor
Definition: parser.hh:28
std::string type
type of additional boundary parameters
Definition: parser.hh:25
[ provides Dune::Grid ]
Definition: yaspgrid.hh:163
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