dune-grid  2.3.1
albertagrid/geometry.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_ALBERTA_GEOMETRY_HH
4 #define DUNE_ALBERTA_GEOMETRY_HH
5 
9 
10 #if HAVE_ALBERTA
11 
12 namespace Dune
13 {
14 
15  // Forward Declarations
16  // --------------------
17 
18  template< int dim, int dimworld >
19  class AlbertaGrid;
20 
21 
22 
23  // AlbertaGridCoordinateReader
24  // ---------------------------
25 
26  template< int codim, class GridImp >
28  {
29  typedef typename remove_const< GridImp >::type Grid;
30 
31  static const int dimension = Grid::dimension;
32  static const int codimension = codim;
33  static const int mydimension = dimension - codimension;
35 
37 
39  typedef FieldVector< ctype, coorddimension > Coordinate;
40 
41  AlbertaGridCoordinateReader ( const GridImp &grid,
42  const ElementInfo &elementInfo,
43  int subEntity )
44  : grid_( grid ),
45  elementInfo_( elementInfo ),
46  subEntity_( subEntity )
47  {}
48 
49  const ElementInfo &elementInfo () const
50  {
51  return elementInfo_;
52  }
53 
54  void coordinate ( int i, Coordinate &x ) const
55  {
56  assert( !elementInfo_ == false );
57  assert( (i >= 0) && (i <= mydimension) );
58 
59  const int k = mapVertices( subEntity_, i );
60  const Alberta::GlobalVector &coord = grid_.getCoord( elementInfo_, k );
61  for( int j = 0; j < coorddimension; ++j )
62  x[ j ] = coord[ j ];
63  }
64 
65  bool hasDeterminant () const
66  {
67  return false;
68  }
69 
70  ctype determinant () const
71  {
72  assert( hasDeterminant() );
73  return ctype( 0 );
74  }
75 
76  private:
77  static int mapVertices ( int subEntity, int i )
78  {
79  return Alberta::MapVertices< dimension, codimension >::apply( subEntity, i );
80  }
81 
82  const Grid &grid_;
83  const ElementInfo &elementInfo_;
84  const int subEntity_;
85  };
86 
87 
88 
89  // AlbertaGridGeometry
90  // -------------------
91 
104  template< int mydim, int cdim, class GridImp >
106  {
108 
109  // remember type of the grid
110  typedef GridImp Grid;
111 
112  // dimension of barycentric coordinates
113  static const int dimbary = mydim + 1;
114 
115  public:
118 
119  static const int dimension = Grid :: dimension;
120  static const int mydimension = mydim;
121  static const int codimension = dimension - mydimension;
122  static const int coorddimension = cdim;
123 
124  typedef FieldVector< ctype, mydimension > LocalCoordinate;
125  typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
126 
127  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
128  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
129 
130  private:
131  static const int numCorners = mydimension + 1;
132 
133  typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
134 
135  public:
137  {
138  invalidate();
139  }
140 
141  template< class CoordReader >
142  AlbertaGridGeometry ( const CoordReader &coordReader )
143  {
144  build( coordReader );
145  }
146 
149  {
150  typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
151  return GeometryType( Topology() );
152  }
153 
155  bool affine () const { return true; }
156 
158  int corners () const
159  {
160  return numCorners;
161  }
162 
164  GlobalCoordinate corner ( const int i ) const
165  {
166  assert( (i >= 0) && (i < corners()) );
167  return coord_[ i ];
168  }
169 
172  {
173  return centroid_;
174  }
175 
179  const GlobalCoordinate &operator[] ( const int i ) const
180  DUNE_DEPRECATED_MSG("Use corner( int i) instead.")
181  {
182  assert( (i >= 0) && (i < corners()) );
183  return coord_[ i ];
184  }
185 
187  GlobalCoordinate global ( const LocalCoordinate &local ) const;
188 
190  LocalCoordinate local ( const GlobalCoordinate &global ) const;
191 
198  {
199  assert( calcedDet_ );
200  return elDet_;
201  }
202 
205  {
206  return integrationElement();
207  }
208 
210  ctype volume () const
211  {
212  return integrationElement() / ctype( Factorial< mydimension >::factorial );
213  }
214 
220  const JacobianTransposed &jacobianTransposed () const;
221 
223  const JacobianTransposed &
225  {
226  return jacobianTransposed();
227  }
228 
235 
239  {
240  return jacobianInverseTransposed();
241  }
242 
243  //***********************************************************************
244  // Methods that not belong to the Interface, but have to be public
245  //***********************************************************************
246 
248  void invalidate ()
249  {
250  builtJT_ = false;
251  builtJTInv_ = false;
252  calcedDet_ = false;
253  }
254 
256  template< class CoordReader >
257  void build ( const CoordReader &coordReader );
258 
259  void print ( std::ostream &out ) const;
260 
261  private:
262  // calculates the volume of the element
263  ctype elDeterminant () const
264  {
266  }
267 
269  CoordMatrix coord_;
270 
272  GlobalCoordinate centroid_;
273 
274  // storage for the transposed of the jacobian
275  mutable JacobianTransposed jT_;
276 
277  // storage for the transposed inverse of the jacboian
278  mutable JacobianInverseTransposed jTInv_;
279 
280  // has jT_ been computed, yet?
281  mutable bool builtJT_;
282  // has jTInv_ been computed, yet?
283  mutable bool builtJTInv_;
284 
285  mutable bool calcedDet_;
286  mutable ctype elDet_;
287  };
288 
289 
290 
291  // AlbertaGridGlobalGeometry
292  // -------------------------
293 
294  template< int mydim, int cdim, class GridImp >
296  : public AlbertaGridGeometry< mydim, cdim, GridImp >
297  {
300 
301  public:
303  : Base()
304  {}
305 
306  template< class CoordReader >
307  AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
308  : Base( coordReader )
309  {}
310  };
311 
312 
313 #if !DUNE_ALBERTA_CACHE_COORDINATES
314  template< int dim, int cdim >
315  class AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > >
316  {
318 
319  // remember type of the grid
321 
322  // dimension of barycentric coordinates
323  static const int dimbary = dim + 1;
324 
326 
327  public:
330 
331  static const int dimension = Grid::dimension;
332  static const int mydimension = dimension;
333  static const int codimension = dimension - mydimension;
334  static const int coorddimension = cdim;
335 
336  typedef FieldVector< ctype, mydimension > LocalCoordinate;
337  typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
338 
339  typedef FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed;
340  typedef FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed;
341 
342  private:
343  static const int numCorners = mydimension + 1;
344 
345  typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
346 
347  public:
349  {
350  invalidate();
351  }
352 
353  template< class CoordReader >
354  AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
355  {
356  build( coordReader );
357  }
358 
361  {
362  typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
363  return GeometryType( Topology() );
364  }
365 
367  int corners () const
368  {
369  return numCorners;
370  }
371 
373  GlobalCoordinate corner ( const int i ) const
374  {
375  assert( (i >= 0) && (i < corners()) );
376  const Alberta::GlobalCoordinate &x = elementInfo_.coordinate( i );
378  for( int j = 0; j < coorddimension; ++j )
379  y[ j ] = x[ j ];
380  return y;
381  }
382 
386  const GlobalCoordinate &operator[] ( const int i ) const
387  DUNE_DEPRECATED_MSG("Use corner( int i) instead.")
388  {
389  return reinterpret_cast< const GlobalCoordinate & >( elementInfo_.coordinate( i ) );
390  }
391 
394  {
395  GlobalCoordinate centroid_ = corner( 0 );
396  for( int i = 1; i < numCorners; ++i )
397  centroid_ += corner( i );
398  centroid_ *= ctype( 1 ) / ctype( numCorners );
399  return centroid_;
400  }
401 
403  GlobalCoordinate global ( const LocalCoordinate &local ) const;
404 
406  LocalCoordinate local ( const GlobalCoordinate &global ) const;
407 
414  {
415  return elementInfo_.geometryCache().integrationElement();
416  }
417 
420  {
421  return integrationElement();
422  }
423 
425  ctype volume () const
426  {
427  return integrationElement() / ctype( Factorial< mydimension >::factorial );
428  }
429 
436  {
437  return elementInfo_.geometryCache().jacobianTransposed();
438  }
439 
441  const JacobianTransposed &
443  {
444  return jacobianTransposed();
445  }
446 
453  {
454  return elementInfo_.geometryCache().jacobianInverseTransposed();
455  }
456 
460  {
461  return jacobianInverseTransposed();
462  }
463 
464  //***********************************************************************
465  // Methods that not belong to the Interface, but have to be public
466  //***********************************************************************
467 
469  void invalidate ()
470  {
471  elementInfo_ = ElementInfo();
472  }
473 
475  template< class CoordReader >
476  void build ( const CoordReader &coordReader )
477  {
478  elementInfo_ = coordReader.elementInfo();
479  }
480 
481  private:
482  ElementInfo elementInfo_;
483  };
484 #endif // #if !DUNE_ALBERTA_CACHE_COORDINATES
485 
486 
487 
488  // AlbertaGridLocalGeometryProvider
489  // --------------------------------
490 
491  template< class Grid >
493  {
495 
496  public:
497  typedef typename Grid::ctype ctype;
498 
499  static const int dimension = Grid::dimension;
500 
501  template< int codim >
502  struct Codim
503  {
505  };
506 
509 
510  static const int numChildren = 2;
511  static const int numFaces = dimension + 1;
512 
513  static const int minFaceTwist = Alberta::Twist< dimension, dimension-1 >::minTwist;
514  static const int maxFaceTwist = Alberta::Twist< dimension, dimension-1 >::maxTwist;
515  static const int numFaceTwists = maxFaceTwist - minFaceTwist + 1;
516 
517  private:
518  struct GeoInFatherCoordReader;
519  struct FaceCoordReader;
520 
522  {
523  buildGeometryInFather();
524  buildFaceGeometry();
525  }
526 
528  {
529  for( int child = 0; child < numChildren; ++child )
530  {
531  delete geometryInFather_[ child ][ 0 ];
532  delete geometryInFather_[ child ][ 1 ];
533  }
534 
535  for( int i = 0; i < numFaces; ++i )
536  {
537  for( int j = 0; j < numFaceTwists; ++j )
538  delete faceGeometry_[ i ][ j ];
539  }
540  }
541 
542  void buildGeometryInFather();
543  void buildFaceGeometry();
544 
545  public:
546  const LocalElementGeometry &
547  geometryInFather ( int child, const int orientation = 1 ) const
548  {
549  assert( (child >= 0) && (child < numChildren) );
550  assert( (orientation == 1) || (orientation == -1) );
551  return *geometryInFather_[ child ][ (orientation + 1) / 2 ];
552  }
553 
554  const LocalFaceGeometry &
555  faceGeometry ( int face, int twist = 0 ) const
556  {
557  assert( (face >= 0) && (face < numFaces) );
558  assert( (twist >= minFaceTwist) && (twist <= maxFaceTwist) );
559  return *faceGeometry_[ face ][ twist - minFaceTwist ];
560  }
561 
562  static const This &instance ()
563  {
564  static This theInstance;
565  return theInstance;
566  }
567 
568  private:
569  template< int codim >
570  static int mapVertices ( int subEntity, int i )
571  {
572  return Alberta::MapVertices< dimension, codim >::apply( subEntity, i );
573  }
574 
575  const LocalElementGeometry *geometryInFather_[ numChildren ][ 2 ];
576  const LocalFaceGeometry *faceGeometry_[ numFaces ][ numFaceTwists ];
577  };
578 
579 
580 
581  // FacadeOptions
582  // -------------
583 
584  namespace FacadeOptions
585  {
586 
587  template< int mydim, int cdim, class Grid >
589  {
590  static const bool v = false;
591  };
592 
593  } // namespace FacadeOptions
594 
595 } // namespace Dune
596 
597 #endif // #if HAVE_ALBERTA
598 
599 #endif // #ifndef DUNE_ALBERTA_GEOMETRY_HH