dune-spgrid  2.7
dgfparser.hh
Go to the documentation of this file.
1 #ifndef DUNE_SPGRID_DGFPARSER_HH
2 #define DUNE_SPGRID_DGFPARSER_HH
3 
4 #include <dune/grid/io/file/dgfparser/dgfparser.hh>
5 #include <dune/grid/spgrid.hh>
6 
7 namespace Dune
8 {
9 
10  namespace dgf
11  {
12 
13  // SPGridParameterBlock
14  // --------------------
15 
16  template< int dim >
18  : public GridParameterBlock
19  {
20  typedef int MultiIndex[ dim ];
21 
22  explicit SPGridParameterBlock( std::istream &in );
23 
24  const MultiIndex &overlap () const
25  {
26  return overlap_;
27  }
28 
29  private:
30  MultiIndex overlap_;
31  };
32 
33 
34 
35  // Implementation of SPGridParameterBlock
36  // --------------------------------------
37 
38  template< int dim >
40  : GridParameterBlock( in )
41  {
42  for( int i = 0; i < dim; ++i )
43  overlap_[ i ] = 0;
44 
45  if( findtoken( "overlap" ) )
46  {
47  int i = 0;
48  for( int x; getnextentry( x ); ++i )
49  {
50  if( x < 0 )
51  DUNE_THROW( DGFException, "Negative overlap specified." );
52  if( i < dim )
53  overlap_[ i ] = x;
54  }
55 
56  if( i < 1 )
57  dwarn << "GridParameterBlock: Found keyword 'overlap' without valid value, defaulting to no overlap." << std::endl;
58  else if( i == 1 )
59  {
60  for( int i = 1; i < dim; ++i )
61  overlap_[ i ] = overlap_[ 0 ];
62  }
63  else if( i != dim )
64  DUNE_THROW( DGFException, "Invalid argument for parameter 'overlap' specified." );
65  }
66  else
67  dwarn << "GridParameterBlock: Parameter 'overlap' not specified, defaulting to no overlap." << std::endl;
68  }
69 
70  } // namespace dgf
71 
72 
73 
74  // DGFGridFactory< SPGrid >
75  // ------------------------
76 
77  template< class ct, int dim, template< int > class Ref, class Comm >
78  class DGFGridFactory< SPGrid< ct, dim, Ref, Comm > >
79  {
80  public:
82 
83  typedef MPIHelper::MPICommunicator MPICommunicatorType;
85 
86  static const int dimension = Grid::dimension;
87  typedef typename Grid::template Codim< 0 >::Entity Element;
88  typedef typename Grid::template Codim< dimension >::Entity Vertex;
89 
90  private:
91  typedef FieldVector< double, dimension > Point;
92 
93  typedef dgf::BoundaryDomBlock BoundaryDomainBlock;
94 
95  public:
96  explicit DGFGridFactory ( std::istream &input,
97  MPICommunicatorType comm = MPIHelper::getCommunicator() );
98 
99  explicit DGFGridFactory ( const std::string &filename,
100  MPICommunicatorType comm = MPIHelper::getCommunicator() );
101 
103  {
104  delete boundaryDomainBlock_;
105  }
106 
107  Grid *grid () const
108  {
109  return grid_;
110  }
111 
112  template< class Intersection >
113  bool wasInserted ( const Intersection &intersection ) const
114  {
115  return false;
116  }
117 
118  template< class Intersection >
119  int boundaryId ( const Intersection &intersection ) const
120  {
121  std::vector< Point > corners;
122  getCorners( intersection.geometry(), corners );
123  const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
124  if( data )
125  return data->id();
126  else
127  return intersection.indexInInside();
128  }
129 
131  {
132  return boundaryDomainBlock_->hasParameter();
133  }
134 
135  template< int codim >
136  int numParameters () const
137  {
138  return 0;
139  }
140 
141  template< class Intersection >
142  const typename DGFBoundaryParameter::type &
143  boundaryParameter ( const Intersection &intersection ) const
144  {
145  std::vector< Point > corners;
146  getCorners( intersection.geometry(), corners );
147  const dgf::DomainData *data = boundaryDomainBlock_->contains( corners );
148  if( data )
149  return data->parameter();
150  else
151  return DGFBoundaryParameter::defaultValue();
152  }
153 
154  template< class Entity >
155  std::vector< double > &parameter ( const Entity &entity )
156  {
157  DUNE_THROW( InvalidStateException,
158  "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
159  }
160 
161  private:
162  void generate ( std::istream &input, const CollectiveCommunication &comm );
163 
164  template< class Geometry >
165  static void getCorners ( const Geometry &geometry, std::vector< Point > &corners )
166  {
167  corners.resize( geometry.corners() );
168  for( int i = 0; i < geometry.corners(); ++i )
169  {
170  const typename Geometry::GlobalCoordinate corner = geometry.corner( i );
171  for( int j = 0; j < dimension; ++j )
172  corners[ i ][ j ] = corner[ j ];
173  }
174  }
175 
176  Grid *grid_;
177  BoundaryDomainBlock *boundaryDomainBlock_;
178  };
179 
180 
181 
182  // Implementation of DGFGridFactory< SPGrid >
183  // ------------------------------------------
184 
185  template< class ct, int dim, template< int > class Ref, class Comm >
186  inline DGFGridFactory< SPGrid< ct, dim, Ref, Comm > >
187  ::DGFGridFactory ( std::istream &input, MPICommunicatorType comm )
188  {
189  generate( input, SPCommunicationTraits< Comm >::comm( comm ) );
190  }
191 
192 
193  template< class ct, int dim, template< int > class Ref, class Comm >
194  inline DGFGridFactory< SPGrid< ct, dim, Ref, Comm > >
195  ::DGFGridFactory ( const std::string &filename, MPICommunicatorType comm )
196  {
197  std::ifstream input( filename.c_str() );
198  if( !input )
199  DUNE_THROW( DGFException, "Unable to open file: " << filename << "." );
200  generate( input, SPCommunicationTraits< Comm >::comm( comm ) );
201  input.close();
202  }
203 
204 
205  template< class ct, int dim, template< int > class Ref, class Comm >
206  inline void
207  DGFGridFactory< SPGrid< ct, dim, Ref, Comm > >
208  ::generate ( std::istream &input, const CollectiveCommunication &comm )
209  {
210  dgf::IntervalBlock intervalBlock( input );
211 
212  if( !intervalBlock.isactive() )
213  DUNE_THROW( DGFException, "DGF stream must contain an interval block to be used with SPGrid< ct, " << dim << " >." );
214  if( intervalBlock.numIntervals() != 1 )
215  DUNE_THROW( DGFException, "SPGrid< ct, " << dim << " > can only handle 1 interval block." );
216 
217  if( intervalBlock.dimw() != dim )
218  DUNE_THROW( DGFException, "SPGrid< ct, " << dim << " cannot handle an interval of dimension " << intervalBlock.dimw() << "." );
219  const dgf::IntervalBlock::Interval &interval = intervalBlock.get( 0 );
220 
221  FieldVector< ct, dim > a, b;
222  int cells[ dim ];
223  for( int i = 0; i < dim; ++i )
224  {
225  a[ i ] = interval.p[ 0 ][ i ];
226  b[ i ] = interval.p[ 1 ][ i ];
227  cells[ i ] = interval.n[ i ];
228  }
229 
230  typedef dgf::PeriodicFaceTransformationBlock::AffineTransformation Transformation;
231  dgf::PeriodicFaceTransformationBlock trafoBlock( input, dim );
232 
233  unsigned int periodic = 0;
234 
235  const int numTrafos = trafoBlock.numTransformations();
236  for( int k = 0; k < numTrafos; ++k )
237  {
238  const Transformation &trafo = trafoBlock.transformation( k );
239 
240  bool identity = true;
241  for( int i = 0; i < dim; ++i )
242  for( int j = 0; j < dim; ++j )
243  identity &= (fabs( (i == j ? 1.0 : 0.0) - trafo.matrix( i, j ) ) < 1e-10);
244  if( !identity )
245  DUNE_THROW( DGFException, "SPGrid< ct, " << dim << " > can only handle shifts as periodic face transformations." );
246 
247  int numDirs = 0;
248  int dir = -1;
249  for( int i = 0; i < dim; ++i )
250  {
251  if( fabs( trafo.shift[ i ] ) < 1e-10 )
252  continue;
253  dir = i;
254  ++numDirs;
255  }
256  const ct shift = fabs( trafo.shift[ dir ] );
257  const ct width = fabs( a[ dir ] - b[ dir ] );
258  if( (numDirs != 1) || (fabs( shift - width ) >= 1e-10) )
259  {
260  std::cerr << "Tranformation '" << trafo
261  << "' does not map boundaries on boundaries." << std::endl;
262  }
263  else
264  periodic |= (1 << dir);
265  }
266 
267  dgf::SPGridParameterBlock< dim > parameter( input );
268 
269  typedef typename Grid::Domain Domain;
270  std::vector< typename Domain::Cube > cubes;
271  cubes.push_back( typename Domain::Cube( a, b ) );
272  Domain domain( cubes, typename Domain::Topology( periodic ) );
273  grid_ = new Grid( domain, cells, parameter.overlap(), comm );
274 
275  boundaryDomainBlock_ = new BoundaryDomainBlock( input, dimension );
276  }
277 
278 
279 
280  // DGFGridInfo< SPGrid >
281  // ---------------------
282 
283  template< class ct, int dim, template< int > class Ref, class Comm >
284  struct DGFGridInfo< SPGrid< ct, dim, Ref, Comm > >
285  {
288 
289  static int refineStepsForHalf ( const RefinementPolicy &policy = RefinementPolicy() )
290  {
291  const unsigned int weight = policy.weight();
292  return (dim + weight - 1) / weight;
293  }
294 
295  static double refineWeight ( const RefinementPolicy &policy = RefinementPolicy() )
296  {
297  return 1.0 / double( 1 << policy.weight() );
298  }
299  };
300 
301 } // namespace Dune
302 
303 #endif // #ifndef DUNE_SPGRID_DGFPARSER_HH
Definition: iostream.hh:7
Definition: communication.hh:24
structured, parallel DUNE grid
Definition: grid.hh:136
Traits::RefinementPolicy RefinementPolicy
Definition: grid.hh:156
Traits::CollectiveCommunication CollectiveCommunication
Definition: grid.hh:168
Definition: dgfparser.hh:19
SPGridParameterBlock(std::istream &in)
Definition: dgfparser.hh:39
const MultiIndex & overlap() const
Definition: dgfparser.hh:24
int MultiIndex[dim]
Definition: dgfparser.hh:20
Grid * grid() const
Definition: dgfparser.hh:107
Grid::template Codim< dimension >::Entity Vertex
Definition: dgfparser.hh:88
const DGFBoundaryParameter::type & boundaryParameter(const Intersection &intersection) const
Definition: dgfparser.hh:143
Grid::template Codim< 0 >::Entity Element
Definition: dgfparser.hh:87
int boundaryId(const Intersection &intersection) const
Definition: dgfparser.hh:119
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgfparser.hh:83
bool wasInserted(const Intersection &intersection) const
Definition: dgfparser.hh:113
bool haveBoundaryParameters() const
Definition: dgfparser.hh:130
SPGrid< ct, dim, Ref, Comm > Grid
Definition: dgfparser.hh:81
int numParameters() const
Definition: dgfparser.hh:136
Grid::CollectiveCommunication CollectiveCommunication
Definition: dgfparser.hh:84
std::vector< double > & parameter(const Entity &entity)
Definition: dgfparser.hh:155
Grid::RefinementPolicy RefinementPolicy
Definition: dgfparser.hh:287
static int refineStepsForHalf(const RefinementPolicy &policy=RefinementPolicy())
Definition: dgfparser.hh:289
SPGrid< ct, dim, Ref, Comm > Grid
Definition: dgfparser.hh:286
static double refineWeight(const RefinementPolicy &policy=RefinementPolicy())
Definition: dgfparser.hh:295