dune-grid  2.4.1
dgfalu.hh
Go to the documentation of this file.
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_DGFPARSERALU_HH
4 #define DUNE_DGFPARSERALU_HH
5 
6 // only include if ALUGrid is used
7 #if HAVE_ALUGRID
8 
9 #include <dune/grid/alugrid.hh>
13 
15 
16 namespace Dune
17 {
18 
19  // External Forward Declarations
20  // -----------------------------
21 
22  template< class GridImp, class IntersectionImp >
23  class Intersection;
24 
25 
26 
27  // DGFGridInfo (specialization for ALUGrid)
28  // ----------------------------------------
29 
31  template<int dimg, int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
32  struct DGFGridInfo< Dune::ALUGrid< dimg, dimw, eltype, refinementtype, Comm > >
33  {
34  static int refineStepsForHalf () { return ( refinementtype == conforming ) ? dimg : 1; }
35  static double refineWeight () { return ( refinementtype == conforming ) ? 0.5 : 1.0/(std::pow( 2.0, double(dimg))); }
36  };
39  // DGFGridFactory for AluGrid
40  // --------------------------
41 
42  // template< int dim, int dimworld > // for a first version
43  template < class G >
44  struct DGFBaseFactory
45  {
46  typedef G Grid;
47  const static int dimension = Grid::dimension;
48  typedef MPIHelper::MPICommunicator MPICommunicatorType;
49  typedef typename Grid::template Codim<0>::Entity Element;
50  typedef typename Grid::template Codim<dimension>::Entity Vertex;
51  typedef Dune::GridFactory<Grid> GridFactory;
52 
53  DGFBaseFactory ()
54  : factory_( ),
55  dgf_( 0, 1 )
56  {}
57 
58  explicit DGFBaseFactory ( MPICommunicatorType comm )
59  : factory_(),
60  dgf_( rank(comm), size(comm) )
61  {}
62 
63  Grid *grid () const
64  {
65  return grid_;
66  }
67 
68  template< class Intersection >
69  bool wasInserted ( const Intersection &intersection ) const
70  {
71  return factory_.wasInserted( intersection );
72  }
73 
74  template< class GG, class II >
75  int boundaryId ( const Intersection< GG, II > & intersection ) const
76  {
77  typedef Dune::Intersection< GG, II > Intersection;
78  typename Intersection::EntityPointer inside = intersection.inside();
79  const typename Intersection::Entity & entity = *inside;
80  const int face = intersection.indexInInside();
81 
82  const ReferenceElement< double, dimension > & refElem =
83  ReferenceElements< double, dimension >::general( entity.type() );
84  int corners = refElem.size( face, 1, dimension );
85  std :: vector< unsigned int > bound( corners );
86  for( int i=0; i < corners; ++i )
87  {
88  const int k = refElem.subEntity( face, 1, i, dimension );
89  bound[ i ] = factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
90  }
91 
92  DuneGridFormatParser::facemap_t::key_type key( bound, false );
93  const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
94  if( pos != dgf_.facemap.end() )
95  return dgf_.facemap.find( key )->second.first;
96  else
97  return (intersection.boundary() ? 1 : 0);
98  }
99 
100  template< class GG, class II >
101  const typename DGFBoundaryParameter::type &
102  boundaryParameter ( const Intersection< GG, II > & intersection ) const
103  {
104  typedef Dune::Intersection< GG, II > Intersection;
105  typename Intersection::EntityPointer inside = intersection.inside();
106  const typename Intersection::Entity & entity = *inside;
107  const int face = intersection.indexInInside();
108 
109  const ReferenceElement< double, dimension > & refElem =
110  ReferenceElements< double, dimension >::general( entity.type() );
111  int corners = refElem.size( face, 1, dimension );
112  std :: vector< unsigned int > bound( corners );
113  for( int i=0; i < corners; ++i )
114  {
115  const int k = refElem.subEntity( face, 1, i, dimension );
116  bound[ i ] = factory_.insertionIndex( *entity.template subEntity< dimension >( k ) );
117  }
118 
119  DuneGridFormatParser::facemap_t::key_type key( bound, false );
120  const DuneGridFormatParser::facemap_t::const_iterator pos = dgf_.facemap.find( key );
121  if( pos != dgf_.facemap.end() )
122  return dgf_.facemap.find( key )->second.second;
123  else
125  }
126 
127  template< int codim >
128  int numParameters () const
129  {
130  if( codim == 0 )
131  return dgf_.nofelparams;
132  else if( codim == dimension )
133  return dgf_.nofvtxparams;
134  else
135  return 0;
136  }
137 
138  // return true if boundary parameters found
139  bool haveBoundaryParameters () const
140  {
141  return dgf_.haveBndParameters;
142  }
143 
144  std::vector< double > &parameter ( const Element &element )
145  {
146  if( numParameters< 0 >() <= 0 )
147  {
148  DUNE_THROW( InvalidStateException,
149  "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
150  }
151  return dgf_.elParams[ factory_.insertionIndex( element ) ];
152  }
153 
154  std::vector< double > &parameter ( const Vertex &vertex )
155  {
156  if( numParameters< dimension >() <= 0 )
157  {
158  DUNE_THROW( InvalidStateException,
159  "Calling DGFGridFactory::parameter is only allowed if there are parameters." );
160  }
161  return dgf_.vtxParams[ factory_.insertionIndex( vertex ) ];
162  }
163 
164  protected:
165  bool generateALUGrid( const ALUGridElementType eltype,
166  std::istream &file,
167  MPICommunicatorType communicator,
168  const std::string &filename );
169 
170  bool generateALU2dGrid( const ALUGridElementType eltype,
171  std::istream &file,
172  MPICommunicatorType communicator,
173  const std::string &filename );
174 
175  static Grid* callDirectly( const char* gridname,
176  const int rank,
177  const char *filename,
178  MPICommunicatorType communicator )
179  {
180  #if ALU3DGRID_PARALLEL
181  // in parallel runs add rank to filename
182  std :: stringstream tmps;
183  tmps << filename << "." << rank;
184  const std :: string &tmp = tmps.str();
185 
186  // if file exits then use it
187  if( fileExists( tmp.c_str() ) )
188  return new Grid( tmp.c_str(), communicator );
189  #endif
190  // for rank 0 we also check the normal file name
191  if( rank == 0 )
192  {
193  if( fileExists( filename ) )
194  return new Grid( filename , communicator );
195 
196  // only throw this exception on rank 0 because
197  // for the other ranks we can still create empty grids
198  DUNE_THROW( GridError, "Unable to create " << gridname << " from '"
199  << filename << "'." );
200  }
201  else
202  dwarn << "WARNING: P[" << rank << "]: Creating empty grid!" << std::endl;
203 
204  // return empty grid on all other processes
205  return new Grid( communicator );
206  }
207  static bool fileExists ( const char *fileName )
208  {
209  std :: ifstream testfile( fileName );
210  if( !testfile )
211  return false;
212  testfile.close();
213  return true;
214  }
215  static int rank( MPICommunicatorType MPICOMM )
216  {
217  int rank = 0;
218 #if HAVE_MPI
219  MPI_Comm_rank( MPICOMM, &rank );
220 #endif
221  return rank;
222  }
223  static int size( MPICommunicatorType MPICOMM )
224  {
225  int size = 1;
226 #if HAVE_MPI
227  MPI_Comm_size( MPICOMM, &size );
228 #endif
229  return size;
230  }
231  Grid *grid_;
232  GridFactory factory_;
233  DuneGridFormatParser dgf_;
234  };
235 
236  template < ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
237  struct DGFGridFactory< ALUGrid<3,3, eltype, refinementtype, Comm > > :
238  public DGFBaseFactory< ALUGrid<3,3, eltype, refinementtype, Comm > >
239  {
240  typedef ALUGrid<3,3, eltype, refinementtype, Comm > DGFGridType;
241  typedef DGFBaseFactory< DGFGridType > BaseType;
242  typedef typename BaseType :: MPICommunicatorType MPICommunicatorType;
243  protected:
244  using BaseType :: grid_;
245  using BaseType :: callDirectly;
246  public:
247  explicit DGFGridFactory ( std::istream &input,
248  MPICommunicatorType comm = MPIHelper::getCommunicator() )
249  : BaseType( comm )
250  {
251  input.clear();
252  input.seekg( 0 );
253  if( !input )
254  DUNE_THROW( DGFException, "Error resetting input stream." );
255  generate( input, comm );
256  }
257 
258  explicit DGFGridFactory ( const std::string &filename,
259  MPICommunicatorType comm = MPIHelper::getCommunicator())
260  : BaseType( comm )
261  {
262  std::ifstream input( filename.c_str() );
263  bool fileFound = input.is_open() ;
264  if( fileFound )
265  fileFound = generate( input, comm, filename );
266 
267  if( ! fileFound )
268  grid_ = callDirectly( "ALUGrid< 3, 3, eltype, ref, comm >", this->rank( comm ), filename.c_str(), comm );
269  }
270 
271  protected:
272  bool generate( std::istream &file, MPICommunicatorType comm, const std::string &filename = "" );
273  };
274 
275  template < int dimw, ALUGridElementType eltype, ALUGridRefinementType refinementtype, class Comm >
276  struct DGFGridFactory< ALUGrid<2, dimw, eltype, refinementtype, Comm > > :
277  public DGFBaseFactory< ALUGrid< 2, dimw, eltype, refinementtype, Comm > >
278  {
279  typedef ALUGrid< 2, dimw, eltype, refinementtype, Comm > DGFGridType;
280  typedef DGFBaseFactory< DGFGridType > BaseType;
281  typedef typename BaseType :: MPICommunicatorType MPICommunicatorType;
282  typedef typename BaseType::Grid Grid;
283  const static int dimension = Grid::dimension;
284  typedef typename BaseType::GridFactory GridFactory;
285 
286  explicit DGFGridFactory ( std::istream &input,
287  MPICommunicatorType comm = MPIHelper::getCommunicator() )
288  : BaseType( comm )
289  {
290  input.clear();
291  input.seekg( 0 );
292  if( !input )
293  DUNE_THROW(DGFException, "Error resetting input stream." );
294  generate( input, comm );
295  }
296 
297  explicit DGFGridFactory ( const std::string &filename,
298  MPICommunicatorType comm = MPIHelper::getCommunicator())
299  : BaseType( comm )
300  {
301  std::ifstream input( filename.c_str() );
302  if( !input )
303  DUNE_THROW( DGFException, "Macrofile '" << filename << "' not found." );
304  if( !generate( input, comm, filename ) )
305  {
306  if( BaseType::fileExists( filename.c_str() ) )
307  grid_ = new Grid( filename );
308  else
309  DUNE_THROW( GridError, "Unable to create a 2d ALUGrid from '" << filename << "'." );
310  }
311  }
312 
313  protected:
314  bool generate( std::istream &file, MPICommunicatorType comm, const std::string &filename = "" );
315  using BaseType::grid_;
316  using BaseType::factory_;
317  using BaseType::dgf_;
318  };
319 
320 }
321 
322 #include "dgfalu.cc"
323 
324 #endif // #if HAVE_ALUGRID
325 
326 #endif // #ifndef DUNE_DGFPARSERALU_HH
Intersection of a mesh entities of codimension 0 ("elements") with a "neighboring" element or with th...
Definition: albertagrid/dgfparser.hh:26
ALUGridElementType
basic element types for ALUGrid
Definition: alugrid/common/declaration.hh:18
Provides base classes for ALUGrid.
Definition: alugrid/common/declaration.hh:20
GridImp::template Codim< 0 >::Entity Entity
Type of entity that this Intersection belongs to.
Definition: common/intersection.hh:185
GridImp::template Codim< 0 >::EntityPointer EntityPointer
Pointer to the type of entities that this Intersection belongs to.
Definition: common/intersection.hh:188
G Grid
Definition: dgfgridfactory.hh:37
MPIHelper::MPICommunicator MPICommunicatorType
Definition: dgfgridfactory.hh:39
Definition: common.hh:179
std::string type
type of additional boundary parameters
Definition: parser.hh:23
static double refineWeight()
static int refineStepsForHalf()
number of globalRefine steps needed to refuce h by 0.5
static const type & defaultValue()
default constructor
Definition: parser.hh:26
The dimension of the grid.
Definition: common/grid.hh:402
Include standard header files.
Definition: agrid.hh:59
ALBERTA EL Element
Definition: misc.hh:51
DGFGridFactory(std::istream &input, MPICommunicatorType comm=MPIHelper::getCommunicator())
Definition: dgfgridfactory.hh:47