dune-grid  2.3.0
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  {
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  {
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
Definition: albertagrid/geometry.hh:295
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Definition: albertagrid/geometry.hh:340
Alberta::Real ctype
Definition: albertagrid/geometry.hh:36
const JacobianTransposed & jacobianTransposed() const
transposed of the geometry mapping's Jacobian
Definition: geometry.cc:53
Codim< 0 >::LocalGeometry LocalElementGeometry
Definition: albertagrid/geometry.hh:507
AlbertaGridCoordinateReader(const GridImp &grid, const ElementInfo &elementInfo, int subEntity)
Definition: albertagrid/geometry.hh:41
GeometryType type() const
obtain the type of reference element
Definition: albertagrid/geometry.hh:148
Alberta::Real ctype
type of coordinates
Definition: albertagrid/geometry.hh:117
Definition: misc.hh:558
int corners() const
number of corner the geometry
Definition: albertagrid/geometry.hh:367
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Definition: albertagrid/geometry.hh:339
Definition: albertagrid/geometry.hh:27
GeometryType type() const
obtain the type of reference element
Definition: albertagrid/geometry.hh:360
const JacobianInverseTransposed & jacobianInverseTransposed() const
transposed inverse of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:452
static const int coorddimension
Definition: albertagrid/geometry.hh:122
static const int dimension
Definition: albertagrid/geometry.hh:31
GlobalCoordinate corner(const int i) const
obtain the i-th corner of this geometry
Definition: albertagrid/geometry.hh:373
AlbertaGridGeometry< dimension-codim, dimension, Grid > LocalGeometry
Definition: albertagrid/geometry.hh:504
ctype integrationElement(const LocalCoordinate &local) const
integration element of the geometry mapping
Definition: albertagrid/geometry.hh:204
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:178
static const int maxFaceTwist
Definition: albertagrid/geometry.hh:514
ctype volume() const
volume of geometry
Definition: albertagrid/geometry.hh:425
static const int mydimension
Definition: albertagrid/geometry.hh:33
static K determinant(const FieldMatrix< K, 0, m > &matrix)
Definition: algebra.hh:28
const GlobalCoordinate & operator[](const int i) const
deprecated way of obtaining the i-th corner
Definition: albertagrid/geometry.hh:179
AlbertaGridGlobalGeometry()
Definition: albertagrid/geometry.hh:302
geometry implementation for AlbertaGrid
Definition: albertagrid/geometry.hh:105
FieldVector< ctype, coorddimension > GlobalCoordinate
Definition: albertagrid/geometry.hh:125
AlbertaGridGlobalGeometry(const CoordReader &coordReader)
Definition: albertagrid/geometry.hh:354
const JacobianInverseTransposed & jacobianInverseTransposed() const
transposed inverse of the geometry mapping's Jacobian
Definition: geometry.cc:71
static const int coorddimension
Definition: albertagrid/geometry.hh:34
GlobalCoordinate corner(const int i) const
obtain the i-th corner of this geometry
Definition: albertagrid/geometry.hh:164
remove_const< GridImp >::type Grid
Definition: albertagrid/geometry.hh:29
static const int numChildren
Definition: albertagrid/geometry.hh:510
ALBERTA REAL_D GlobalVector
Definition: misc.hh:47
GlobalCoordinate global(const LocalCoordinate &local) const
map a point from the refence element to the geometry
Definition: geometry.cc:32
Traits class determining whether the Dune::Geometry facade class stores the implementation object by ...
Definition: common/geometry.hh:46
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
transposed inverse of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:238
The dimension of the world the grid lives in.
Definition: common/grid.hh:406
AlbertaGridGeometry(const CoordReader &coordReader)
Definition: albertagrid/geometry.hh:142
void invalidate()
invalidate the geometry
Definition: albertagrid/geometry.hh:248
GlobalCoordinate center() const
return center of geometry
Definition: albertagrid/geometry.hh:393
static const This & instance()
Definition: albertagrid/geometry.hh:562
Definition: albertagrid/geometry.hh:492
static const int dimension
Definition: albertagrid/geometry.hh:499
The dimension of the grid.
Definition: common/grid.hh:400
[ provides Dune::Grid ]
Definition: agrid.hh:137
provides a wrapper for ALBERTA's el_info structure
ctype volume() const
volume of geometry
Definition: albertagrid/geometry.hh:210
static const bool v
Whether to store by reference.
Definition: common/geometry.hh:49
ctype integrationElement() const
integration element of the geometry mapping
Definition: albertagrid/geometry.hh:413
Grid abstract base classThis class is the base class for all grid implementations. Although no virtual functions are used we call it abstract since its methods do not contain an implementation but forward to the methods of the derived class via the Barton-Nackman trick.
Definition: common/grid.hh:386
Codim< 1 >::LocalGeometry LocalFaceGeometry
Definition: albertagrid/geometry.hh:508
Wrapper and interface classes for element geometries.
void build(const CoordReader &coordReader)
build the geometry from a coordinate reader
Definition: geometry.cc:88
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
transposed of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:442
static const int numFaceTwists
Definition: albertagrid/geometry.hh:515
const LocalFaceGeometry & faceGeometry(int face, int twist=0) const
Definition: albertagrid/geometry.hh:555
AlbertaGridGeometry()
Definition: albertagrid/geometry.hh:136
void print(std::ostream &out) const
Definition: geometry.cc:18
GlobalCoordinate center() const
return center of geometry
Definition: albertagrid/geometry.hh:171
const JacobianInverseTransposed & jacobianInverseTransposed(const LocalCoordinate &local) const
transposed inverse of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:459
ct ctype
Define type used for coordinates in grid module.
Definition: common/grid.hh:546
ctype integrationElement() const
integration element of the geometry mapping
Definition: albertagrid/geometry.hh:197
Alberta::ElementInfo< dimension > ElementInfo
Definition: albertagrid/geometry.hh:38
LocalCoordinate local(const GlobalCoordinate &global) const
map a point from the geometry to the reference element
Definition: geometry.cc:43
AlbertaGridGlobalGeometry(const CoordReader &coordReader)
Definition: albertagrid/geometry.hh:307
static const int codimension
Definition: albertagrid/geometry.hh:121
FieldMatrix< ctype, mydimension, coorddimension > JacobianTransposed
Definition: albertagrid/geometry.hh:127
const JacobianTransposed & jacobianTransposed() const
transposed of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:435
Alberta::Real ctype
type of coordinates
Definition: albertagrid/geometry.hh:329
const JacobianTransposed & jacobianTransposed(const LocalCoordinate &local) const
transposed of the geometry mapping's Jacobian
Definition: albertagrid/geometry.hh:224
ctype determinant() const
Definition: albertagrid/geometry.hh:70
Grid::ctype ctype
Definition: albertagrid/geometry.hh:497
ctype integrationElement(const LocalCoordinate &local) const
integration element of the geometry mapping
Definition: albertagrid/geometry.hh:419
FieldVector< ctype, coorddimension > Coordinate
Definition: albertagrid/geometry.hh:39
const LocalElementGeometry & geometryInFather(int child, const int orientation=1) const
Definition: albertagrid/geometry.hh:547
bool affine() const
returns always true since we only have simplices
Definition: albertagrid/geometry.hh:155
static const int dimension
Definition: albertagrid/geometry.hh:119
FieldVector< ctype, mydimension > LocalCoordinate
Definition: albertagrid/geometry.hh:124
ALBERTA REAL Real
Definition: misc.hh:45
bool hasDeterminant() const
Definition: albertagrid/geometry.hh:65
static const int mydimension
Definition: albertagrid/geometry.hh:120
void build(const CoordReader &coordReader)
build the geometry from a coordinate reader
Definition: albertagrid/geometry.hh:476
const ElementInfo & elementInfo() const
Definition: albertagrid/geometry.hh:49
FieldVector< ctype, mydimension > LocalCoordinate
Definition: albertagrid/geometry.hh:336
FieldMatrix< ctype, coorddimension, mydimension > JacobianInverseTransposed
Definition: albertagrid/geometry.hh:128
Definition: misc.hh:470
static const int minFaceTwist
Definition: albertagrid/geometry.hh:513
static const int codimension
Definition: albertagrid/geometry.hh:32
FieldVector< ctype, coorddimension > GlobalCoordinate
Definition: albertagrid/geometry.hh:337
void invalidate()
invalidate the geometry
Definition: albertagrid/geometry.hh:469
int corners() const
number of corner the geometry
Definition: albertagrid/geometry.hh:158
void coordinate(int i, Coordinate &x) const
Definition: albertagrid/geometry.hh:54
void abs(const DofVectorPointer< int > &dofVector)
Definition: dofvector.hh:332
static const int numFaces
Definition: albertagrid/geometry.hh:511
Definition: albertagrid/geometry.hh:502