dune-grid  2.2.1
albertagrid/geometry.hh
Go to the documentation of this file.
1 #ifndef DUNE_ALBERTA_GEOMETRY_HH
2 #define DUNE_ALBERTA_GEOMETRY_HH
3 
4 #include <dune/geometry/genericgeometry/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  // AlbertaGridCoordStorage
90  // -----------------------
91 
92  template< class CoordTraits, class Topology, unsigned int dimW >
94  {
96 
97  public:
98  static const unsigned int size = Topology::numCorners;
99 
100  static const unsigned int dimWorld = dimW;
101 
102  typedef typename CoordTraits::template Vector< dimWorld >::type
104 
105  template< class SubTopology >
106  struct SubStorage
107  {
109  };
110 
111  private:
112  GlobalCoordinate coords_[ size ];
113 
114  public:
115  template< class CoordReader >
116  explicit AlbertaGridCornerStorage ( const CoordReader &coordReader )
117  {
118  for( unsigned int i = 0; i < size; ++i )
119  coordReader.coordinate( i, coords_[ i ] );
120  }
121 
122  template< class Mapping, unsigned int codim >
123  explicit AlbertaGridCornerStorage ( const GenericGeometry::SubMappingCoords< Mapping, codim > &coords )
124  {
125  for( unsigned int i = 0; i < size; ++i )
126  coords_[ i ] = coords[ i ];
127  }
128 
129  const GlobalCoordinate &operator[] ( unsigned int i ) const
130  {
131  return coords_[ i ];
132  }
133  };
134 
135 
136 
137  // AlbertaGridGeometryTraits
138  // -------------------------
139 
140  template< class GridImp, int cdim >
142  {
143  typedef typename remove_const< GridImp >::type Grid;
144 
145  typedef GenericGeometry::DuneCoordTraits< Alberta::Real > CoordTraits;
146 
147  static const int dimGrid = Grid::dimension;
148  static const int dimWorld = cdim;
149 
150  static const bool hybrid = false;
151  static const unsigned int topologyId = GenericGeometry::SimplexTopology< dimGrid >::type::id;
152 
153  template< class Topology >
154  struct Mapping
155  {
157  typedef GenericGeometry::CornerMapping< CoordTraits, Topology, dimWorld, CornerStorage > type;
158  };
159 
160  struct Caching
161  {
162  static const GenericGeometry::EvaluationType evaluateJacobianTransposed = GenericGeometry::ComputeOnDemand;
163  static const GenericGeometry::EvaluationType evaluateJacobianInverseTransposed = GenericGeometry::ComputeOnDemand;
164  static const GenericGeometry::EvaluationType evaluateIntegrationElement = GenericGeometry::ComputeOnDemand;
165  static const GenericGeometry::EvaluationType evaluateNormal = GenericGeometry::ComputeOnDemand;
166  };
167 
168  struct UserData {};
169  };
170 
171 
172 
173  // AlbertaGridGeometry
174  // -------------------
175 
176 #if DUNE_ALBERTA_USE_GENERICGEOMETRY
177  template< int mydim, int cdim, class GridImp >
178  class AlbertaGridGeometry
179  : public GenericGeometry::BasicGeometry
180  < mydim, AlbertaGridGeometryTraits< GridImp, cdim > >
181  {
183  typedef GenericGeometry::BasicGeometry
185  Base;
186 
187  public:
190  : Base ()
191  {}
192 
193  AlbertaGridGeometry ( const This &other )
194  : Base ( other )
195  {}
196 
197  template< class CoordReader >
198  AlbertaGridGeometry ( const CoordReader &coordReader )
199  : Base( GeometryType( GenericGeometry::SimplexTopology< mydim >::type::id, mydim ), coordReader )
200  {}
201 
202  template< class CoordReader >
203  void build ( const CoordReader &coordReader )
204  {
205  (*this) = AlbertaGridGeometry( coordReader );
206  }
207  };
208 #endif // #if DUNE_ALBERTA_USE_GENERICGEOMETRY
209 
210 #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
211 
223  template< int mydim, int cdim, class GridImp >
225  {
227 
228  // remember type of the grid
229  typedef GridImp Grid;
230 
231  // dimension of barycentric coordinates
232  static const int dimbary = mydim + 1;
233 
234  public:
237 
238  static const int dimension = Grid :: dimension;
239  static const int mydimension = mydim;
240  static const int codimension = dimension - mydimension;
241  static const int coorddimension = cdim;
242 
243  typedef FieldVector< ctype, mydimension > LocalCoordinate;
244  typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
245 
246  typedef FieldMatrix< ctype, mydimension, coorddimension >
248  typedef FieldMatrix< ctype, coorddimension, mydimension >
251 
252  private:
253  static const int numCorners = mydimension + 1;
254 
255  typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
256 
257  public:
259  {
260  invalidate();
261  }
262 
263  template< class CoordReader >
264  AlbertaGridGeometry ( const CoordReader &coordReader )
265  {
266  build( coordReader );
267  }
268 
271  {
272  typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
273  return GeometryType( Topology() );
274  }
275 
277  bool affine () const { return true; }
278 
280  int corners () const
281  {
282  return numCorners;
283  }
284 
286  GlobalCoordinate corner ( const int i ) const
287  {
288  assert( (i >= 0) && (i < corners()) );
289  return coord_[ i ];
290  }
291 
294  {
295  return centroid_;
296  }
297 
299  const GlobalCoordinate &operator[] ( const int i ) const
300  {
301  assert( (i >= 0) && (i < corners()) );
302  return coord_[ i ];
303  }
304 
306  GlobalCoordinate global ( const LocalCoordinate &local ) const;
307 
309  LocalCoordinate local ( const GlobalCoordinate &global ) const;
310 
317  {
318  assert( calcedDet_ );
319  return elDet_;
320  }
321 
324  {
325  return integrationElement();
326  }
327 
329  ctype volume () const
330  {
332  }
333 
339  const JacobianTransposed &jacobianTransposed () const;
340 
342  const JacobianTransposed &
344  {
345  return jacobianTransposed();
346  }
347 
354 
358  {
359  return jacobianInverseTransposed();
360  }
361 
362  //***********************************************************************
363  // Methods that not belong to the Interface, but have to be public
364  //***********************************************************************
365 
367  void invalidate ()
368  {
369  builtJT_ = false;
370  builtJTInv_ = false;
371  calcedDet_ = false;
372  }
373 
375  template< class CoordReader >
376  void build ( const CoordReader &coordReader );
377 
378  void print ( std::ostream &out ) const;
379 
380  private:
381  // calculates the volume of the element
382  ctype elDeterminant () const
383  {
385  }
386 
388  CoordMatrix coord_;
389 
391  GlobalCoordinate centroid_;
392 
393  // storage for the transposed of the jacobian
394  mutable JacobianTransposed jT_;
395 
396  // storage for the transposed inverse of the jacboian
397  mutable JacobianInverseTransposed jTInv_;
398 
399  // has jT_ been computed, yet?
400  mutable bool builtJT_;
401  // has jTInv_ been computed, yet?
402  mutable bool builtJTInv_;
403 
404  mutable bool calcedDet_;
405  mutable ctype elDet_;
406  };
407 #endif // #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
408 
409 
410 
411  // AlbertaGridGlobalGeometry
412  // -------------------------
413 
414  template< int mydim, int cdim, class GridImp >
416  : public AlbertaGridGeometry< mydim, cdim, GridImp >
417  {
420 
421  public:
423  : Base()
424  {}
425 
426  template< class CoordReader >
427  AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
428  : Base( coordReader )
429  {}
430  };
431 
432 
433 #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
434 #if !DUNE_ALBERTA_CACHE_COORDINATES
435  template< int dim, int cdim >
436  class AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > >
437  {
439 
440  // remember type of the grid
442 
443  // dimension of barycentric coordinates
444  static const int dimbary = dim + 1;
445 
447 
448  public:
451 
452  static const int dimension = Grid::dimension;
453  static const int mydimension = dimension;
454  static const int codimension = dimension - mydimension;
455  static const int coorddimension = cdim;
456 
457  typedef FieldVector< ctype, mydimension > LocalCoordinate;
458  typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
459 
460  typedef FieldMatrix< ctype, mydimension, coorddimension >
462  typedef FieldMatrix< ctype, coorddimension, mydimension >
465 
466  private:
467  static const int numCorners = mydimension + 1;
468 
469  typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
470 
471  public:
473  {
474  invalidate();
475  }
476 
477  template< class CoordReader >
478  AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
479  {
480  build( coordReader );
481  }
482 
485  {
486  typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
487  return GeometryType( Topology() );
488  }
489 
491  int corners () const
492  {
493  return numCorners;
494  }
495 
497  GlobalCoordinate corner ( const int i ) const
498  {
499  assert( (i >= 0) && (i < corners()) );
500  const Alberta::GlobalCoordinate &x = elementInfo_.coordinate( i );
502  for( int j = 0; j < coorddimension; ++j )
503  y[ j ] = x[ j ];
504  return y;
505  }
506 
508  const GlobalCoordinate &operator[] ( const int i ) const
509  {
510  return reinterpret_cast< const GlobalCoordinate & >( elementInfo_.coordinate( i ) );
511  }
512 
515  {
516  GlobalCoordinate centroid_ = corner( 0 );
517  for( int i = 1; i < numCorners; ++i )
518  centroid_ += corner( i );
519  centroid_ *= ctype( 1 ) / ctype( numCorners );
520  return centroid_;
521  }
522 
524  GlobalCoordinate global ( const LocalCoordinate &local ) const;
525 
527  LocalCoordinate local ( const GlobalCoordinate &global ) const;
528 
535  {
536  return elementInfo_.geometryCache().integrationElement();
537  }
538 
541  {
542  return integrationElement();
543  }
544 
546  ctype volume () const
547  {
549  }
550 
557  {
558  return elementInfo_.geometryCache().jacobianTransposed();
559  }
560 
562  const JacobianTransposed &
564  {
565  return jacobianTransposed();
566  }
567 
574  {
575  return elementInfo_.geometryCache().jacobianInverseTransposed();
576  }
577 
581  {
582  return jacobianInverseTransposed();
583  }
584 
585  //***********************************************************************
586  // Methods that not belong to the Interface, but have to be public
587  //***********************************************************************
588 
590  void invalidate ()
591  {
592  elementInfo_ = ElementInfo();
593  }
594 
596  template< class CoordReader >
597  void build ( const CoordReader &coordReader )
598  {
599  elementInfo_ = coordReader.elementInfo();
600  }
601 
602  private:
603  ElementInfo elementInfo_;
604  };
605 #endif // #if !DUNE_ALBERTA_CACHE_COORDINATES
606 #endif // #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
607 
608 
609 
610  // AlbertaGridLocalGeometryProvider
611  // --------------------------------
612 
613  template< class Grid >
615  {
617 
618  public:
619  typedef typename Grid::ctype ctype;
620 
621  static const int dimension = Grid::dimension;
622 
623  template< int codim >
624  struct Codim
625  {
627  };
628 
631 
632  static const int numChildren = 2;
633  static const int numFaces = dimension + 1;
634 
635  static const int minFaceTwist = Alberta::Twist< dimension, dimension-1 >::minTwist;
636  static const int maxFaceTwist = Alberta::Twist< dimension, dimension-1 >::maxTwist;
637  static const int numFaceTwists = maxFaceTwist - minFaceTwist + 1;
638 
639  private:
640  struct GeoInFatherCoordReader;
641  struct FaceCoordReader;
642 
644  {
645  buildGeometryInFather();
646  buildFaceGeometry();
647  }
648 
650  {
651  for( int child = 0; child < numChildren; ++child )
652  {
653  delete geometryInFather_[ child ][ 0 ];
654  delete geometryInFather_[ child ][ 1 ];
655  }
656 
657  for( int i = 0; i < numFaces; ++i )
658  {
659  for( int j = 0; j < numFaceTwists; ++j )
660  delete faceGeometry_[ i ][ j ];
661  }
662  }
663 
664  void buildGeometryInFather();
665  void buildFaceGeometry();
666 
667  public:
668  const LocalElementGeometry &
669  geometryInFather ( int child, const int orientation = 1 ) const
670  {
671  assert( (child >= 0) && (child < numChildren) );
672  assert( (orientation == 1) || (orientation == -1) );
673  return *geometryInFather_[ child ][ (orientation + 1) / 2 ];
674  }
675 
676  const LocalFaceGeometry &
677  faceGeometry ( int face, int twist = 0 ) const
678  {
679  assert( (face >= 0) && (face < numFaces) );
680  assert( (twist >= minFaceTwist) && (twist <= maxFaceTwist) );
681  return *faceGeometry_[ face ][ twist - minFaceTwist ];
682  }
683 
684  static const This &instance ()
685  {
686  static This theInstance;
687  return theInstance;
688  }
689 
690  private:
691  template< int codim >
692  static int mapVertices ( int subEntity, int i )
693  {
695  }
696 
697  const LocalElementGeometry *geometryInFather_[ numChildren ][ 2 ];
698  const LocalFaceGeometry *faceGeometry_[ numFaces ][ numFaceTwists ];
699  };
700 
701 
702 
703  // FacadeOptions
704  // -------------
705 
706  namespace FacadeOptions
707  {
708 
709  template< int mydim, int cdim, class Grid >
711  {
712  static const bool v = false;
713  };
714 
715  } // namespace FacadeOptions
716 
717 } // namespace Dune
718 
719 #endif // #if HAVE_ALBERTA
720 
721 #endif // #ifndef DUNE_ALBERTA_GEOMETRY_HH