Dune Core Modules (2.6.0)

dgfyasp.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_DGFPARSERYASP_HH
4#define DUNE_DGFPARSERYASP_HH
5
6#include <dune/grid/common/intersection.hh>
7#include <dune/grid/yaspgrid.hh>
8#include "dgfparser.hh"
9
10namespace Dune
11{
12
13 // External Forward Declarations
14 // -----------------------------
15
16 template< class GridImp, class IntersectionImp >
17 class Intersection;
18
19
20 namespace dgf
21 {
22
36 class YaspGridParameterBlock
37 : public GridParameterBlock
38 {
39 protected:
40 int _overlap; // overlap for YaspGrid
41
42 public:
44 YaspGridParameterBlock( std::istream &in )
45 : GridParameterBlock( in ),
46 _overlap( 0 ) // default value
47 {
48 // check overlap
49 if( findtoken( "overlap" ) )
50 {
51 int x;
52 if( getnextentry(x) ) _overlap = x;
53 else
54 {
55 dwarn << "GridParameterBlock: found keyword `overlap' but no value, defaulting to `" << _overlap <<"' !\n";
56 }
57
58 if (_overlap < 0)
59 {
60 DUNE_THROW(DGFException,"Negative overlap specified!");
61 }
62 }
63 else
64 {
65 dwarn << "YaspGridParameterBlock: Parameter 'overlap' not specified, "
66 << "defaulting to '" << _overlap << "'." << std::endl;
67 }
68
69 }
70
72 int overlap () const
73 {
74 return _overlap;
75 }
76
77 };
78
79 }
80
84 template <typename ctype, int dim>
85 struct DGFGridFactory< YaspGrid<dim, EquidistantCoordinates<ctype, dim> > >
86 {
87 typedef YaspGrid<dim, EquidistantCoordinates<ctype, dim> > Grid;
88 const static int dimension = Grid::dimension;
89 typedef MPIHelper::MPICommunicator MPICommunicatorType;
90
91 private:
92 typedef FieldVector< ctype, dimension > Point;
93 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
94
95 public:
96 explicit DGFGridFactory ( std::istream &input,
97 MPICommunicatorType comm = MPIHelper::getCommunicator() )
98 {
99 generate( input, comm );
100 }
101
102 explicit DGFGridFactory ( const std::string &filename,
103 MPICommunicatorType comm = MPIHelper::getCommunicator() )
104 {
105 std::ifstream input( filename.c_str() );
106 if( !input )
107 DUNE_THROW( DGFException, "Error: Macrofile '" << filename << "' not found" );
108 generate( input, comm );
109 }
110
111 ~DGFGridFactory ()
112 {
113 delete boundaryDomainBlock_;
114 }
115
116 Grid *grid() const
117 {
118 return grid_;
119 }
120
121 template <class Intersection>
122 bool wasInserted(const Intersection &intersection) const
123 {
124 return false;
125 }
126
127 template <class Intersection>
128 int boundaryId(const Intersection &intersection) const
129 {
130 if( boundaryDomainBlock_->isactive() )
131 {
132 std::vector< Point > corners;
133 getCorners( intersection.geometry(), corners );
134 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
135 if( data )
136 return data->id();
137 else
138 return intersection.indexInInside();
139 }
140 else
141 return intersection.indexInInside();
142 }
143
144 template< int codim >
145 int numParameters () const
146 {
147 return 0;
148 }
149
150 // return true if boundary parameters found
151 bool haveBoundaryParameters () const
152 {
153 return boundaryDomainBlock_->hasParameter();
154 }
155
156 template< class GG, class II >
157 const typename DGFBoundaryParameter::type &
158 boundaryParameter ( const Intersection< GG, II > & intersection ) const
159 {
160 if( haveBoundaryParameters() )
161 {
162 std::vector< Point > corners;
163 getCorners( intersection.geometry(), corners );
164 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
165 if( data )
166 return data->parameter();
167 else
168 return DGFBoundaryParameter::defaultValue();
169 }
170 else
171 return DGFBoundaryParameter::defaultValue();
172 }
173
174 template< class Entity >
175 std::vector<double> &parameter ( const Entity &entity )
176 {
177 return emptyParam;
178 }
179
180 private:
181 void generate( std::istream &gridin, MPICommunicatorType comm );
182
183 template< class Geometry >
184 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
185 {
186 corners.resize( geometry.corners() );
187 for( int i = 0; i < geometry.corners(); ++i )
188 {
189 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
190 for( int j = 0; j < dimension; ++j )
191 corners[ i ][ j ] = corner[ j ];
192 }
193 }
194
195 Grid *grid_;
196 dgf::BoundaryDomBlock *boundaryDomainBlock_;
197 std::vector<double> emptyParam;
198 };
199
200 // generate YaspGrid from the provided DGF
201 template< typename ctype, int dim >
202 inline void DGFGridFactory< YaspGrid< dim , EquidistantCoordinates<ctype, dim> > >
203 ::generate ( std::istream &gridin, MPICommunicatorType comm )
204 {
205 using std::abs;
206 dgf::IntervalBlock intervalBlock( gridin );
207
208 if( !intervalBlock.isactive() )
209 DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
210
211 if( intervalBlock.numIntervals() != 1 )
212 DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
213
214 if( intervalBlock.dimw() != dim )
215 {
216 DUNE_THROW( DGFException,
217 "Cannot read an interval of dimension " << intervalBlock.dimw()
218 << " into a YaspGrid< " << dim << " >." );
219 }
220
221 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
222
223 FieldVector<ctype, dim> lang;
224 std::array<int,dim> anz;
225 for( int i = 0; i < dim; ++i )
226 {
227 // check that start point is 0.0
228 if( abs( interval.p[ 0 ][ i ] ) > 1e-10 )
229 {
230 DUNE_THROW( DGFException,
231 "YaspGrid cannot handle grids with non-zero left lower corner." );
232 }
233
234 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
235 anz[ i ] = interval.n[ i ];
236 }
237
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 )
243 {
244 const Transformation &trafo = trafoBlock.transformation( k );
245
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);
250 if( !identity )
251 DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
252
253 int numDirs = 0;
254 int dir = -1;
255 for( int i = 0; i < dim; ++i )
256 {
257 if( abs( trafo.shift[ i ] ) < 1e-10 )
258 continue;
259 dir = i;
260 ++numDirs;
261 }
262 if( (numDirs != 1) || (abs( abs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
263 {
264 std::cerr << "Tranformation '" << trafo
265 << "' does not map boundaries on boundaries." << std::endl;
266 }
267 else
268 per[ dir ] = true;
269 }
270
271 // get grid parameters
272 dgf::YaspGridParameterBlock grdParam( gridin );
273
274 grid_ = new YaspGrid< dim , EquidistantCoordinates<ctype, dim> >( lang, anz, per, grdParam.overlap(), comm );
275
276 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
277 }
278
279 template <typename ctype, int dim>
280 struct DGFGridInfo< YaspGrid<dim , EquidistantCoordinates<ctype, dim> > > {
281 static int refineStepsForHalf() {return 1;}
282 static double refineWeight() {return std::pow(0.5,dim);}
283 };
284
288 template <typename ctype, int dim>
289 struct DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > >
290 {
291 typedef YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > Grid;
292 const static int dimension = Grid::dimension;
293 typedef MPIHelper::MPICommunicator MPICommunicatorType;
294
295 private:
296 typedef FieldVector< ctype, dimension > Point;
297 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
298
299 public:
300 explicit DGFGridFactory ( std::istream &input,
301 MPICommunicatorType comm = MPIHelper::getCommunicator() )
302 {
303 generate( input, comm );
304 }
305
306 explicit DGFGridFactory ( const std::string &filename,
307 MPICommunicatorType comm = MPIHelper::getCommunicator() )
308 {
309 std::ifstream input( filename.c_str() );
310 generate( input, comm );
311 }
312
313 ~DGFGridFactory ()
314 {
315 delete boundaryDomainBlock_;
316 }
317
318 Grid *grid() const
319 {
320 return grid_;
321 }
322
323 template <class Intersection>
324 bool wasInserted(const Intersection &intersection) const
325 {
326 return false;
327 }
328
329 template <class Intersection>
330 int boundaryId(const Intersection &intersection) const
331 {
332 if( boundaryDomainBlock_->isactive() )
333 {
334 std::vector< Point > corners;
335 getCorners( intersection.geometry(), corners );
336 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
337 if( data )
338 return data->id();
339 else
340 return intersection.indexInInside();
341 }
342 else
343 return intersection.indexInInside();
344 }
345
346 template< int codim >
347 int numParameters () const
348 {
349 return 0;
350 }
351
352 // return true if boundary parameters found
353 bool haveBoundaryParameters () const
354 {
355 return boundaryDomainBlock_->hasParameter();
356 }
357
358 template< class GG, class II >
359 const typename DGFBoundaryParameter::type &
360 boundaryParameter ( const Intersection< GG, II > & intersection ) const
361 {
362 if( haveBoundaryParameters() )
363 {
364 std::vector< Point > corners;
365 getCorners( intersection.geometry(), corners );
366 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
367 if( data )
368 return data->parameter();
369 else
370 return DGFBoundaryParameter::defaultValue();
371 }
372 else
373 return DGFBoundaryParameter::defaultValue();
374 }
375
376 template< class Entity >
377 std::vector<double> &parameter ( const Entity &entity )
378 {
379 return emptyParam;
380 }
381
382 private:
383 void generate( std::istream &gridin, MPICommunicatorType comm );
384
385 template< class Geometry >
386 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
387 {
388 corners.resize( geometry.corners() );
389 for( int i = 0; i < geometry.corners(); ++i )
390 {
391 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
392 for( int j = 0; j < dimension; ++j )
393 corners[ i ][ j ] = corner[ j ];
394 }
395 }
396
397 Grid *grid_;
398 dgf::BoundaryDomBlock *boundaryDomainBlock_;
399 std::vector<double> emptyParam;
400 };
401
402 // generate YaspGrid from the provided DGF
403 template< typename ctype, int dim >
404 inline void DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > >
405 ::generate ( std::istream &gridin, MPICommunicatorType comm )
406 {
407 using std::abs;
408 dgf::IntervalBlock intervalBlock( gridin );
409
410 if( !intervalBlock.isactive() )
411 DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
412
413 if( intervalBlock.numIntervals() != 1 )
414 DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
415
416 if( intervalBlock.dimw() != dim )
417 {
418 DUNE_THROW( DGFException,
419 "Cannot read an interval of dimension "
420 << intervalBlock.dimw()
421 << " into a YaspGrid< " << dim << " >." );
422 }
423
424 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
425
426 FieldVector<ctype, dim> lower;
427 FieldVector<ctype, dim> upper;
428 std::array<int,dim> anz;
429 for( int i = 0; i < dim; ++i )
430 {
431 lower[ i ] = interval.p[ 0 ][ i ];
432 upper[ i ] = interval.p[ 1 ][ i ];
433 anz[ i ] = interval.n[ i ];
434 }
435
436 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
437 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
438 std::bitset< dim > periodic;
439 const int numTrafos = trafoBlock.numTransformations();
440 for( int k = 0; k < numTrafos; ++k )
441 {
442 const Transformation &trafo = trafoBlock.transformation( k );
443
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);
448 if( !identity )
449 DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
450
451 int numDirs = 0;
452 int dir = -1;
453 for( int currentDir = 0; currentDir < dim; ++currentDir )
454 {
455 if( abs( trafo.shift[ currentDir ] ) > 1e-10 )
456 {
457 dir = currentDir;
458 ++numDirs;
459 }
460 }
461 if ( (numDirs != 1)
462 || (abs( abs( trafo.shift[ dir ] ) - abs( upper[ dir ] - lower[ dir ] ) ) >= 1e-10) )
463 {
464 std::cerr << "Tranformation '" << trafo
465 << "' does not map boundaries on boundaries." << std::endl;
466 }
467 else
468 {
469 periodic[ dir ] = true;
470 }
471 }
472
473 // get grid parameters
474 dgf::YaspGridParameterBlock grdParam( gridin );
475
476 grid_ = new YaspGrid< dim, EquidistantOffsetCoordinates<ctype, dim> >
477 ( lower, upper, anz, periodic, grdParam.overlap(), comm );
478
479 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
480 }
481
482 template <typename ctype, int dim>
483 struct DGFGridInfo< YaspGrid<dim, EquidistantOffsetCoordinates<ctype, dim> > > {
484 static int refineStepsForHalf() {return 1;}
485 static double refineWeight() {return std::pow(0.5,dim);}
486 };
487
488}
489#endif // #ifndef DUNE_DGFPARSERYASP_HH
Dune::DGFException
exception class for IO errors in the DGF parser
Definition: dgfexception.hh:14
Dune::Entity
Wrapper class for entities.
Definition: entity.hh:64
Dune::EquidistantCoordinates
Container for equidistant coordinates in a YaspGrid.
Definition: coordinates.hh:27
Dune::EquidistantOffsetCoordinates
Container for equidistant coordinates in a YaspGrid with non-trivial origin.
Definition: coordinates.hh:125
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Dune::Geometry
Wrapper class for geometries.
Definition: geometry.hh:67
Dune::Geometry::corner
GlobalCoordinate corner(int i) const
Obtain a corner of the geometry.
Definition: geometry.hh:153
Dune::Geometry::corners
int corners() const
Return the number of corners of the reference element.
Definition: geometry.hh:139
Dune::Grid::dimension
@ dimension
The dimension of the grid.
Definition: grid.hh:387
Dune::Intersection
Intersection of a mesh entity of codimension 0 ("element") with a "neighboring" element or with the d...
Definition: intersection.hh:162
Dune::Intersection::geometry
Geometry geometry() const
geometrical information about the intersection in global coordinates.
Definition: intersection.hh:333
Dune::Intersection::indexInInside
int indexInInside() const
Local index of codim 1 entity in the inside() entity where intersection is contained in.
Definition: intersection.hh:356
Dune::MPIHelper::MPICommunicator
MPI_Comm MPICommunicator
The type of the mpi communicator.
Definition: mpihelper.hh:178
Dune::MPIHelper::getCommunicator
static MPICommunicator getCommunicator()
get the default communicator
Definition: mpihelper.hh:186
Dune::YaspGrid
[ provides Dune::Grid ]
Definition: yaspgrid.hh:165
Dune::dgf::GridParameterBlock
Common Grid parameters.
Definition: gridparameter.hh:33
Dune::dgf::YaspGridParameterBlock
Grid parameters for YaspGrid.
Definition: dgfyasp.hh:38
Dune::dgf::YaspGridParameterBlock::YaspGridParameterBlock
YaspGridParameterBlock(std::istream &in)
constructor taking istream
Definition: dgfyasp.hh:44
Dune::dgf::YaspGridParameterBlock::overlap
int overlap() const
get dimension of world found in block
Definition: dgfyasp.hh:72
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::dwarn
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:159
Dune
Dune namespace.
Definition: alignedallocator.hh:10
Dune::DGFBoundaryParameter::defaultValue
static const type & defaultValue()
default constructor
Definition: parser.hh:26
Dune::DGFBoundaryParameter::type
std::string type
type of additional boundary parameters
Definition: parser.hh:23
Dune::DGFGridInfo::refineWeight
static double refineWeight()
Dune::DGFGridInfo::refineStepsForHalf
static int refineStepsForHalf()
number of globalRefine steps needed to refuce h by 0.5
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)