dune-grid  2.2.0
dgfyasp.hh
Go to the documentation of this file.
1 #ifndef DUNE_DGFPARSERYASP_HH
2 #define DUNE_DGFPARSERYASP_HH
3 
5 #include <dune/grid/yaspgrid.hh>
6 #include "dgfparser.hh"
7 
8 namespace Dune
9 {
10  // forward declaration
11  // -------------------
12 
13  template< class GridImp, template< class > class IntersectionImp >
14  class Intersection;
15 
16 
17  namespace dgf
18  {
19 
34  : public GridParameterBlock
35  {
36  protected:
37  int _overlap; // overlap for YaspGrid
38 
39  public:
41  YaspGridParameterBlock( std::istream &in )
42  : GridParameterBlock( in ),
43  _overlap( 0 ) // default value
44  {
45  // check overlap
46  if( findtoken( "overlap" ) )
47  {
48  int x;
49  if( getnextentry(x) ) _overlap = x;
50  else
51  {
52  dwarn << "GridParameterBlock: found keyword `overlap' but no value, defaulting to `" << _overlap <<"' !\n";
53  }
54 
55  if (_overlap < 0)
56  {
57  DUNE_THROW(DGFException,"Negative overlap specified!");
58  }
59  }
60  else
61  {
62  dwarn << "YaspGridParameterBlock: Parameter 'overlap' not specified, "
63  << "defaulting to '" << _overlap << "'." << std::endl;
64  }
65 
66  }
67 
69  int overlap () const
70  {
71  return _overlap;
72  }
73 
74  };
75 
76  }
77 
78  template <int dim>
79  struct DGFGridFactory< YaspGrid<dim> >
80  {
82  const static int dimension = Grid::dimension;
83  typedef MPIHelper::MPICommunicator MPICommunicatorType;
84 
85  private:
86  typedef FieldVector< double, dimension > Point;
88 
89  public:
90  explicit DGFGridFactory ( std::istream &input,
91  MPICommunicatorType comm = MPIHelper::getCommunicator() )
92  {
93  generate( input, comm );
94  }
95 
96  explicit DGFGridFactory ( const std::string &filename,
97  MPICommunicatorType comm = MPIHelper::getCommunicator() )
98  {
99  std::ifstream input( filename.c_str() );
100  generate( input, comm );
101  }
102 
104  {
105  delete boundaryDomainBlock_;
106  }
107 
108  Grid *grid() const
109  {
110  return grid_;
111  }
112 
113  template <class Intersection>
114  bool wasInserted(const Intersection &intersection) const
115  {
116  return false;
117  }
118  template <class Intersection>
119  int boundaryId(const Intersection &intersection) const
120  {
121  if( boundaryDomainBlock_->isactive() )
122  {
123  std::vector< Point > corners;
124  getCorners( intersection.geometry(), corners );
125  const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
126  if( data )
127  return data->id();
128  else
129  return intersection.indexInInside();
130  }
131  else
132  return intersection.indexInInside();
133  }
134 
135  template< int codim >
136  int numParameters () const
137  {
138  return 0;
139  }
140 
141  // return true if boundary parameters found
143  {
144  return boundaryDomainBlock_->hasParameter();
145  }
146 
147  template < class GG, template < class > class II >
148  const typename DGFBoundaryParameter::type &
149  boundaryParameter ( const Intersection< GG, II > & intersection ) const
150  {
151  if( haveBoundaryParameters() )
152  {
153  std::vector< Point > corners;
154  getCorners( intersection.geometry(), corners );
155  const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
156  if( data )
157  return data->parameter();
158  else
160  }
161  else
163  }
164 
165  template< class Entity >
166  std::vector<double> &parameter ( const Entity &entity )
167  {
168  return emptyParam;
169  }
170 
171  private:
172  void generate( std::istream &gridin, MPICommunicatorType comm );
173 
174  template< class Geometry >
175  static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
176  {
177  corners.resize( geometry.corners() );
178  for( int i = 0; i < geometry.corners(); ++i )
179  {
180  const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
181  for( int j = 0; j < dimension; ++j )
182  corners[ i ][ j ] = corner[ j ];
183  }
184  }
185 
186  Grid *grid_;
187  dgf::BoundaryDomBlock *boundaryDomainBlock_;
188  std::vector<double> emptyParam;
189  };
190 
191  template< int dim >
192  inline void DGFGridFactory< YaspGrid< dim > >
193  ::generate ( std::istream &gridin, MPICommunicatorType comm )
194  {
195  dgf::IntervalBlock intervalBlock( gridin );
196 
197  if( !intervalBlock.isactive() )
198  DUNE_THROW( DGFException, "YaspGrid can only be created from an interval block." );
199 
200  if( intervalBlock.numIntervals() != 1 )
201  DUNE_THROW( DGFException, "YaspGrid can only handle 1 interval block." );
202 
203  if( intervalBlock.dimw() != dim )
204  {
205  DUNE_THROW( DGFException,
206  "Cannot read an interval of dimension " << intervalBlock.dimw()
207  << "into a YaspGrid< " << dim << " >." );
208  }
209 
210  const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
211 
212  FieldVector<double,dim> lang;
213  FieldVector<int,dim> anz;
214  for( int i = 0; i < dim; ++i )
215  {
216  // check that start point is 0.0
217  if( fabs( interval.p[ 0 ][ i ] ) > 1e-10 )
218  {
219  DUNE_THROW( DGFException,
220  "YaspGrid cannot handle grids with non-zero left lower corner." );
221  }
222 
223  lang[ i ] = interval.p[ 1 ][ i ] - interval.p[ 0 ][ i ];
224  anz[ i ] = interval.n[ i ];
225  }
226 
227  typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
228  dgf::PeriodicFaceTransformationBlock trafoBlock( gridin, dim );
229  FieldVector< bool, dim > per( false );
230  const int numTrafos = trafoBlock.numTransformations();
231  for( int k = 0; k < numTrafos; ++k )
232  {
233  const Transformation &trafo = trafoBlock.transformation( k );
234 
235  bool identity = true;
236  for( int i = 0; i < dim; ++i )
237  for( int j = 0; j < dim; ++j )
238  identity &= (fabs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
239  if( !identity )
240  DUNE_THROW( DGFException, "YaspGrid can only handle shifts as periodic face transformations." );
241 
242  int numDirs = 0;
243  int dir = -1;
244  for( int i = 0; i < dim; ++i )
245  {
246  if( fabs( trafo.shift[ i ] ) < 1e-10 )
247  continue;
248  dir = i;
249  ++numDirs;
250  }
251  if( (numDirs != 1) || (fabs( fabs( trafo.shift[ dir ] ) - lang[ dir ] ) >= 1e-10) )
252  {
253  std::cerr << "Tranformation '" << trafo
254  << "' does not map boundaries on boundaries." << std::endl;
255  }
256  else
257  per[ dir ] = true;
258  }
259 
260  // get grid parameters
261  dgf::YaspGridParameterBlock grdParam( gridin );
262 
263 #if HAVE_MPI
264  grid_ = new YaspGrid< dim >( comm, lang, anz, per, grdParam.overlap() );
265 #else
266  grid_ = new YaspGrid< dim >( lang, anz, per, grdParam.overlap() );
267 #endif
268 
269  boundaryDomainBlock_ = new dgf::BoundaryDomBlock( gridin, dimension );
270  }
271 
272  template <int dim>
273  struct DGFGridInfo< YaspGrid<dim> > {
274  static int refineStepsForHalf() {return 1;}
275  static double refineWeight() {return pow(0.5,dim);}
276  };
277 
278 
279 }
280 #endif // #ifndef DUNE_DGFPARSERYASP_HH