Dune Core Modules (2.4.1)

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 <int dim>
85 struct DGFGridFactory< YaspGrid<dim /*, EquidistantCoordinates<double, dim> */> >
86 {
87 typedef YaspGrid<dim> Grid;
88 const static int dimension = Grid::dimension;
89 typedef MPIHelper::MPICommunicator MPICommunicatorType;
90
91 private:
92 typedef FieldVector< double, 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 generate( input, comm );
107 }
108
109 ~DGFGridFactory ()
110 {
111 delete boundaryDomainBlock_;
112 }
113
114 Grid *grid() const
115 {
116 return grid_;
117 }
118
119 template <class Intersection>
120 bool wasInserted(const Intersection &intersection) const
121 {
122 return false;
123 }
124
125 template <class Intersection>
126 int boundaryId(const Intersection &intersection) const
127 {
128 if( boundaryDomainBlock_->isactive() )
129 {
130 std::vector< Point > corners;
131 getCorners( intersection.geometry(), corners );
132 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
133 if( data )
134 return data->id();
135 else
136 return intersection.indexInInside();
137 }
138 else
139 return intersection.indexInInside();
140 }
141
142 template< int codim >
143 int numParameters () const
144 {
145 return 0;
146 }
147
148 // return true if boundary parameters found
149 bool haveBoundaryParameters () const
150 {
151 return boundaryDomainBlock_->hasParameter();
152 }
153
154 template< class GG, class II >
155 const typename DGFBoundaryParameter::type &
156 boundaryParameter ( const Intersection< GG, II > & intersection ) const
157 {
158 if( haveBoundaryParameters() )
159 {
160 std::vector< Point > corners;
161 getCorners( intersection.geometry(), corners );
162 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
163 if( data )
164 return data->parameter();
165 else
166 return DGFBoundaryParameter::defaultValue();
167 }
168 else
169 return DGFBoundaryParameter::defaultValue();
170 }
171
172 template< class Entity >
173 std::vector<double> &parameter ( const Entity &entity )
174 {
175 return emptyParam;
176 }
177
178 private:
179 void generate( std::istream &gridin, MPICommunicatorType comm );
180
181 template< class Geometry >
182 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
183 {
184 corners.resize( geometry.corners() );
185 for( int i = 0; i < geometry.corners(); ++i )
186 {
187 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
188 for( int j = 0; j < dimension; ++j )
189 corners[ i ][ j ] = corner[ j ];
190 }
191 }
192
193 Grid *grid_;
194 dgf::BoundaryDomBlock *boundaryDomainBlock_;
195 std::vector<double> emptyParam;
196 };
197
198 // generate YaspGrid from the provided DGF
199 template< int dim >
200 inline void DGFGridFactory< YaspGrid< dim /*, EquidistantCoordinates<double, dim> */> >
201 ::generate ( std::istream &gridin, MPICommunicatorType comm )
202 {
203 dgf::IntervalBlock intervalBlock( gridin );
204
205 if( !intervalBlock.isactive() )
206 DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
207
208 if( intervalBlock.numIntervals() != 1 )
209 DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
210
211 if( intervalBlock.dimw() != dim )
212 {
213 DUNE_THROW( DGFException,
214 "Cannot read an interval of dimension " << intervalBlock.dimw()
215 << " into a YaspGrid< " << dim << " >." );
216 }
217
218 const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
219
220 FieldVector<double,dim> lang;
221 std::array<int,dim> anz;
222 for( int i = 0; i < dim; ++i )
223 {
224 // check that start point is 0.0
225 if( std::abs( interval.p[ 0 ][ i ] ) > 1e-10 )
226 {
227 DUNE_THROW( DGFException,
228 "YaspGrid cannot handle grids with non-zero left lower corner." );
229 }
230
231 lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
232 anz[ i ] = interval.n[ i ];
233 }
234
235 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
236 dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
237 std::bitset< dim > per;
238 const int numTrafos = trafoBlock.numTransformations();
239 for( int k = 0; k < numTrafos; ++k )
240 {
241 const Transformation &trafo = trafoBlock.transformation( k );
242
243 bool identity = true;
244 for( int i = 0; i < dim; ++i )
245 for( int j = 0; j < dim; ++j )
246 identity &= (std::abs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
247 if( !identity )
248 DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
249
250 int numDirs = 0;
251 int dir = -1;
252 for( int i = 0; i < dim; ++i )
253 {
254 if( std::abs( trafo.shift[ i ] ) < 1e-10 )
255 continue;
256 dir = i;
257 ++numDirs;
258 }
259 if( (numDirs != 1) || (std::abs( std::abs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
260 {
261 std::cerr << "Tranformation '" << trafo
262 << "' does not map boundaries on boundaries." << std::endl;
263 }
264 else
265 per[ dir ] = true;
266 }
267
268 // get grid parameters
269 dgf::YaspGridParameterBlock grdParam( gridin );
270
271 grid_ = new YaspGrid< dim >( lang, anz, per, grdParam.overlap(), comm );
272
273 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
274 }
275
276 template <int dim>
277 struct DGFGridInfo< YaspGrid<dim /*, EquidistantCoordinates<double, dim> */> > {
278 static int refineStepsForHalf() {return 1;}
279 static double refineWeight() {return std::pow(0.5,dim);}
280 };
281
285 template <int dim>
286 struct DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> > >
287 {
288 typedef YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> > Grid;
289 const static int dimension = Grid::dimension;
290 typedef MPIHelper::MPICommunicator MPICommunicatorType;
291
292 private:
293 typedef FieldVector< double, dimension > Point;
294 typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
295
296 public:
297 explicit DGFGridFactory ( std::istream &input,
298 MPICommunicatorType comm = MPIHelper::getCommunicator() )
299 {
300 generate( input, comm );
301 }
302
303 explicit DGFGridFactory ( const std::string &filename,
304 MPICommunicatorType comm = MPIHelper::getCommunicator() )
305 {
306 std::ifstream input( filename.c_str() );
307 generate( input, comm );
308 }
309
310 ~DGFGridFactory ()
311 {
312 delete boundaryDomainBlock_;
313 }
314
315 Grid *grid() const
316 {
317 return grid_;
318 }
319
320 template <class Intersection>
321 bool wasInserted(const Intersection &intersection) const
322 {
323 return false;
324 }
325
326 template <class Intersection>
327 int boundaryId(const Intersection &intersection) const
328 {
329 if( boundaryDomainBlock_->isactive() )
330 {
331 std::vector< Point > corners;
332 getCorners( intersection.geometry(), corners );
333 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
334 if( data )
335 return data->id();
336 else
337 return intersection.indexInInside();
338 }
339 else
340 return intersection.indexInInside();
341 }
342
343 template< int codim >
344 int numParameters () const
345 {
346 return 0;
347 }
348
349 // return true if boundary parameters found
350 bool haveBoundaryParameters () const
351 {
352 return boundaryDomainBlock_->hasParameter();
353 }
354
355 template< class GG, class II >
356 const typename DGFBoundaryParameter::type &
357 boundaryParameter ( const Intersection< GG, II > & intersection ) const
358 {
359 if( haveBoundaryParameters() )
360 {
361 std::vector< Point > corners;
362 getCorners( intersection.geometry(), corners );
363 const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
364 if( data )
365 return data->parameter();
366 else
367 return DGFBoundaryParameter::defaultValue();
368 }
369 else
370 return DGFBoundaryParameter::defaultValue();
371 }
372
373 template< class Entity >
374 std::vector<double> &parameter ( const Entity &entity )
375 {
376 return emptyParam;
377 }
378
379 private:
380 void generate( std::istream &gridin, MPICommunicatorType comm );
381
382 template< class Geometry >
383 static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
384 {
385 corners.resize( geometry.corners() );
386 for( int i = 0; i < geometry.corners(); ++i )
387 {
388 const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
389 for( int j = 0; j < dimension; ++j )
390 corners[ i ][ j ] = corner[ j ];
391 }
392 }
393
394 Grid *grid_;
395 dgf::BoundaryDomBlock *boundaryDomainBlock_;
396 std::vector<double> emptyParam;
397 };
398
399 // generate YaspGrid from the provided DGF
400 template< int dim >
401 inline void DGFGridFactory< YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> > >
402 ::generate ( std::istream &gridin, MPICommunicatorType comm )
403 {
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<double,dim> lower;
423 FieldVector<double,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
432 typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
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 &= (std::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( std::abs( trafo.shift[ currentDir ] ) > 1e-10 )
452 {
453 dir = currentDir;
454 ++numDirs;
455 }
456 }
457 if ( (numDirs != 1)
458 || (std::abs( std::abs( trafo.shift[ dir ] ) - std::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<double, dim> >
473 ( lower, upper, anz, periodic, grdParam.overlap(), comm );
474
475 boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
476 }
477
478 template <int dim>
479 struct DGFGridInfo< YaspGrid<dim, EquidistantOffsetCoordinates<double, dim> > > {
480 static int refineStepsForHalf() {return 1;}
481 static double refineWeight() {return std::pow(0.5,dim);}
482 };
483
484}
485#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:62
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:94
Dune::Geometry
Wrapper class for geometries.
Definition: geometry.hh:66
Dune::Geometry::corner
GlobalCoordinate corner(int i) const
Obtain a corner of the geometry.
Definition: geometry.hh:156
Dune::Geometry::corners
int corners() const
Return the number of corners of the reference element.
Definition: geometry.hh:142
Dune::Grid::dimension
@ dimension
The dimension of the grid.
Definition: grid.hh:402
Dune::Intersection
Intersection of a mesh entities of codimension 0 ("elements") with a "neighboring" element or with th...
Definition: intersection.hh:161
Dune::Intersection::geometry
Geometry geometry() const
geometrical information about the intersection in global coordinates.
Definition: intersection.hh:373
Dune::Intersection::indexInInside
int indexInInside() const
Local index of codim 1 entity in the inside() entity where intersection is contained in.
Definition: intersection.hh:393
Dune::MPIHelper::MPICommunicator
MPI_Comm MPICommunicator
The type of the mpi communicator.
Definition: mpihelper.hh:173
Dune::MPIHelper::getCommunicator
static MPICommunicator getCommunicator()
get the default communicator
Definition: mpihelper.hh:181
Dune::YaspGrid
[ provides Dune::Grid ]
Definition: yaspgrid.hh:166
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:243
Dune::dwarn
DWarnType dwarn(std::cerr)
Stream for warnings indicating problems.
Definition: stdstreams.hh:159
Dune
Dune namespace.
Definition: alignment.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 12, 23:30, 2024)