dune-grid  2.2.0
albertagrid/geometry.hh
Go to the documentation of this file.
00001 #ifndef DUNE_ALBERTA_GEOMETRY_HH
00002 #define DUNE_ALBERTA_GEOMETRY_HH
00003 
00004 #include <dune/geometry/genericgeometry/geometry.hh>
00005 
00006 #include <dune/grid/common/geometry.hh>
00007 #include <dune/grid/albertagrid/misc.hh>
00008 #include <dune/grid/albertagrid/elementinfo.hh>
00009 
00010 #if HAVE_ALBERTA
00011 
00012 namespace Dune
00013 {
00014 
00015   // Forward Declarations
00016   // --------------------
00017 
00018   template< int dim, int dimworld >
00019   class AlbertaGrid;
00020 
00021 
00022 
00023   // AlbertaGridCoordinateReader
00024   // ---------------------------
00025 
00026   template< int codim, class GridImp >
00027   struct AlbertaGridCoordinateReader
00028   {
00029     typedef typename remove_const< GridImp >::type Grid;
00030 
00031     static const int dimension = Grid::dimension;
00032     static const int codimension = codim;
00033     static const int mydimension = dimension - codimension;
00034     static const int coorddimension = Grid::dimensionworld;
00035 
00036     typedef Alberta::Real ctype;
00037 
00038     typedef Alberta::ElementInfo< dimension > ElementInfo;
00039     typedef FieldVector< ctype, coorddimension > Coordinate;
00040 
00041     AlbertaGridCoordinateReader ( const GridImp &grid,
00042                                   const ElementInfo &elementInfo,
00043                                   int subEntity )
00044     : grid_( grid ),
00045       elementInfo_( elementInfo ),
00046       subEntity_( subEntity )
00047     {}
00048 
00049     const ElementInfo &elementInfo () const
00050     {
00051       return elementInfo_;
00052     }
00053 
00054     void coordinate ( int i, Coordinate &x ) const
00055     {
00056       assert( !elementInfo_ == false );
00057       assert( (i >= 0) && (i <= mydimension) );
00058 
00059       const int k = mapVertices( subEntity_, i );
00060       const Alberta::GlobalVector &coord = grid_.getCoord( elementInfo_, k );
00061       for( int j = 0; j < coorddimension; ++j )
00062         x[ j ] = coord[ j ];
00063     }
00064 
00065     bool hasDeterminant () const
00066     {
00067       return false;
00068     }
00069 
00070     ctype determinant () const
00071     {
00072       assert( hasDeterminant() );
00073       return ctype( 0 );
00074     }
00075 
00076   private:
00077     static int mapVertices ( int subEntity, int i )
00078     {
00079       return Alberta::MapVertices< dimension, codimension >::apply( subEntity, i );
00080     }
00081 
00082     const Grid &grid_;
00083     const ElementInfo &elementInfo_;
00084     const int subEntity_;
00085   };
00086 
00087 
00088 
00089   // AlbertaGridCoordStorage
00090   // -----------------------
00091 
00092   template< class CoordTraits, class Topology, unsigned int dimW >
00093   class AlbertaGridCornerStorage
00094   {
00095     typedef AlbertaGridCornerStorage< CoordTraits, Topology, dimW > This;
00096 
00097   public:
00098     static const unsigned int size = Topology::numCorners;
00099 
00100     static const unsigned int dimWorld = dimW;
00101 
00102     typedef typename CoordTraits::template Vector< dimWorld >::type
00103       GlobalCoordinate;
00104 
00105     template< class SubTopology >
00106     struct SubStorage
00107     {
00108       typedef AlbertaGridCornerStorage< CoordTraits, SubTopology, dimWorld > type;
00109     };
00110 
00111   private:
00112     GlobalCoordinate coords_[ size ];
00113 
00114   public:
00115     template< class CoordReader >
00116     explicit AlbertaGridCornerStorage ( const CoordReader &coordReader )
00117     {
00118       for( unsigned int i = 0; i < size; ++i )
00119         coordReader.coordinate( i, coords_[ i ] );
00120     }
00121 
00122     template< class Mapping, unsigned int codim >
00123     explicit AlbertaGridCornerStorage ( const GenericGeometry::SubMappingCoords< Mapping, codim > &coords )
00124     {
00125       for( unsigned int i = 0; i < size; ++i )
00126         coords_[ i ] = coords[ i ];
00127     }
00128 
00129     const GlobalCoordinate &operator[] ( unsigned int i ) const
00130     {
00131       return coords_[ i ];
00132     }
00133   };
00134 
00135 
00136 
00137   // AlbertaGridGeometryTraits
00138   // -------------------------
00139 
00140   template< class GridImp, int cdim >
00141   struct AlbertaGridGeometryTraits
00142   {
00143     typedef typename remove_const< GridImp >::type Grid;
00144 
00145     typedef GenericGeometry::DuneCoordTraits< Alberta::Real > CoordTraits;
00146 
00147     static const int dimGrid = Grid::dimension;
00148     static const int dimWorld = cdim;
00149 
00150     static const bool hybrid = false;
00151     static const unsigned int topologyId = GenericGeometry::SimplexTopology< dimGrid >::type::id;
00152 
00153     template< class Topology >
00154     struct Mapping
00155     {
00156       typedef AlbertaGridCornerStorage< CoordTraits, Topology, dimWorld > CornerStorage;
00157       typedef GenericGeometry::CornerMapping< CoordTraits, Topology, dimWorld, CornerStorage > type;
00158     };
00159 
00160     struct Caching
00161     {
00162       static const GenericGeometry::EvaluationType evaluateJacobianTransposed = GenericGeometry::ComputeOnDemand;
00163       static const GenericGeometry::EvaluationType evaluateJacobianInverseTransposed = GenericGeometry::ComputeOnDemand;
00164       static const GenericGeometry::EvaluationType evaluateIntegrationElement = GenericGeometry::ComputeOnDemand;
00165       static const GenericGeometry::EvaluationType evaluateNormal = GenericGeometry::ComputeOnDemand;
00166     };
00167     
00168     struct UserData {};
00169   };
00170 
00171 
00172 
00173   // AlbertaGridGeometry
00174   // -------------------
00175 
00176 #if DUNE_ALBERTA_USE_GENERICGEOMETRY
00177   template< int mydim, int cdim, class GridImp >
00178   class AlbertaGridGeometry
00179   : public GenericGeometry::BasicGeometry
00180     < mydim, AlbertaGridGeometryTraits< GridImp, cdim > >
00181   { 
00182     typedef AlbertaGridGeometry< mydim, cdim, GridImp > This;
00183     typedef GenericGeometry::BasicGeometry
00184       < mydim, AlbertaGridGeometryTraits< GridImp, cdim > >
00185       Base;
00186 
00187     public:
00189     AlbertaGridGeometry ()
00190     : Base ()
00191     {}
00192 
00193     AlbertaGridGeometry ( const This &other )
00194     : Base ( other )
00195     {}
00196 
00197     template< class CoordReader >
00198     AlbertaGridGeometry ( const CoordReader &coordReader )
00199     : Base( GeometryType( GenericGeometry::SimplexTopology< mydim >::type::id, mydim ), coordReader )
00200     {}
00201 
00202     template< class CoordReader >
00203     void build ( const CoordReader &coordReader )
00204     {
00205       (*this) = AlbertaGridGeometry( coordReader );
00206     }
00207   };
00208 #endif // #if DUNE_ALBERTA_USE_GENERICGEOMETRY
00209 
00210 #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
00211 
00223   template< int mydim, int cdim, class GridImp >
00224   class AlbertaGridGeometry
00225   { 
00226     typedef AlbertaGridGeometry< mydim, cdim, GridImp > This;
00227 
00228     // remember type of the grid
00229     typedef GridImp Grid;
00230 
00231     // dimension of barycentric coordinates
00232     static const int dimbary = mydim + 1;
00233 
00234   public:
00236     typedef Alberta::Real ctype;
00237 
00238     static const int dimension = Grid :: dimension;
00239     static const int mydimension = mydim;
00240     static const int codimension = dimension - mydimension;
00241     static const int coorddimension = cdim;
00242 
00243     typedef FieldVector< ctype, mydimension > LocalCoordinate;
00244     typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
00245 
00246     typedef FieldMatrix< ctype, mydimension, coorddimension >
00247       JacobianTransposed;
00248     typedef FieldMatrix< ctype, coorddimension, mydimension >
00249       JacobianInverseTransposed;
00250     typedef JacobianInverseTransposed  Jacobian;
00251 
00252   private:
00253     static const int numCorners = mydimension + 1;
00254 
00255     typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
00256 
00257   public:
00258     AlbertaGridGeometry ()
00259     {
00260       invalidate();
00261     }
00262 
00263     template< class CoordReader >
00264     AlbertaGridGeometry ( const CoordReader &coordReader )
00265     {
00266       build( coordReader );
00267     }
00268 
00270     GeometryType type () const
00271     {
00272       typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
00273       return GeometryType( Topology() );
00274     }
00275 
00277     bool affine () const { return true; }
00278 
00280     int corners () const
00281     {
00282       return numCorners;
00283     }
00284 
00286     GlobalCoordinate corner ( const int i ) const
00287     {
00288       assert( (i >= 0) && (i < corners()) );
00289       return coord_[ i ];
00290     }
00291 
00293     GlobalCoordinate center () const
00294     {
00295       return centroid_;
00296     }
00297 
00299     const GlobalCoordinate &operator[] ( const int i ) const
00300     {
00301       assert( (i >= 0) && (i < corners()) );
00302       return coord_[ i ];
00303     }
00304 
00306     GlobalCoordinate global ( const LocalCoordinate &local ) const;
00307   
00309     LocalCoordinate local ( const GlobalCoordinate &global ) const;
00310 
00316     ctype integrationElement () const
00317     {
00318       assert( calcedDet_ );
00319       return elDet_;
00320     }
00321 
00323     ctype integrationElement ( const LocalCoordinate &local ) const
00324     {
00325       return integrationElement();
00326     }
00327 
00329     ctype volume () const
00330     {
00331       return integrationElement() / ctype( Factorial< mydimension >::factorial );
00332     }
00333 
00339     const JacobianTransposed &jacobianTransposed () const;
00340 
00342     const JacobianTransposed &
00343     jacobianTransposed ( const LocalCoordinate &local ) const
00344     {
00345       return jacobianTransposed();
00346     }
00347  
00353     const JacobianInverseTransposed &jacobianInverseTransposed () const;
00354 
00356     const JacobianInverseTransposed &
00357     jacobianInverseTransposed ( const LocalCoordinate &local ) const
00358     {
00359       return jacobianInverseTransposed();
00360     }
00361 
00362     //***********************************************************************
00363     //  Methods that not belong to the Interface, but have to be public
00364     //***********************************************************************
00365 
00367     void invalidate ()
00368     {
00369       builtJT_ = false;
00370       builtJTInv_ = false;
00371       calcedDet_ = false;
00372     }
00373 
00375     template< class CoordReader >
00376     void build ( const CoordReader &coordReader );
00377 
00378     void print ( std::ostream &out ) const;
00379 
00380   private:
00381     // calculates the volume of the element 
00382     ctype elDeterminant () const
00383     {
00384       return std::abs( Alberta::determinant( jacobianTransposed() ) );
00385     }
00386 
00388     CoordMatrix coord_;
00389 
00391     GlobalCoordinate centroid_;
00392 
00393     // storage for the transposed of the jacobian
00394     mutable JacobianTransposed jT_;
00395 
00396     // storage for the transposed inverse of the jacboian
00397     mutable JacobianInverseTransposed jTInv_;
00398 
00399     // has jT_ been computed, yet?
00400     mutable bool builtJT_;
00401     // has jTInv_ been computed, yet?
00402     mutable bool builtJTInv_;
00403 
00404     mutable bool calcedDet_; 
00405     mutable ctype elDet_; 
00406   };
00407 #endif // #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
00408 
00409 
00410 
00411   // AlbertaGridGlobalGeometry
00412   // -------------------------
00413 
00414   template< int mydim, int cdim, class GridImp >
00415   class AlbertaGridGlobalGeometry
00416   : public AlbertaGridGeometry< mydim, cdim, GridImp >
00417   {
00418     typedef AlbertaGridGlobalGeometry< mydim, cdim, GridImp > This;
00419     typedef AlbertaGridGeometry< mydim, cdim, GridImp > Base;
00420 
00421   public:
00422     AlbertaGridGlobalGeometry ()
00423     : Base()
00424     {}
00425 
00426     template< class CoordReader >
00427     AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
00428     : Base( coordReader )
00429     {}
00430   };
00431 
00432 
00433 #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
00434 #if !DUNE_ALBERTA_CACHE_COORDINATES
00435   template< int dim, int cdim >
00436   class AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > >
00437   { 
00438     typedef AlbertaGridGlobalGeometry< dim, cdim, const AlbertaGrid< dim, cdim > > This;
00439 
00440     // remember type of the grid
00441     typedef AlbertaGrid< dim, cdim > Grid;
00442 
00443     // dimension of barycentric coordinates
00444     static const int dimbary = dim + 1;
00445 
00446     typedef Alberta::ElementInfo< dim > ElementInfo;
00447 
00448   public:
00450     typedef Alberta::Real ctype;
00451 
00452     static const int dimension = Grid::dimension;
00453     static const int mydimension = dimension;
00454     static const int codimension = dimension - mydimension;
00455     static const int coorddimension = cdim;
00456 
00457     typedef FieldVector< ctype, mydimension > LocalCoordinate;
00458     typedef FieldVector< ctype, coorddimension > GlobalCoordinate;
00459 
00460     typedef FieldMatrix< ctype, mydimension, coorddimension >
00461       JacobianTransposed;
00462     typedef FieldMatrix< ctype, coorddimension, mydimension >
00463       JacobianInverseTransposed;
00464     typedef JacobianInverseTransposed  Jacobian;
00465 
00466   private:
00467     static const int numCorners = mydimension + 1;
00468 
00469     typedef FieldMatrix< ctype, numCorners, coorddimension > CoordMatrix;
00470 
00471   public:
00472     AlbertaGridGlobalGeometry ()
00473     {
00474       invalidate();
00475     }
00476 
00477     template< class CoordReader >
00478     AlbertaGridGlobalGeometry ( const CoordReader &coordReader )
00479     {
00480       build( coordReader );
00481     }
00482 
00484     GeometryType type () const
00485     {
00486       typedef typename GenericGeometry::SimplexTopology< mydimension >::type Topology;
00487       return GeometryType( Topology() );
00488     }
00489 
00491     int corners () const
00492     {
00493       return numCorners;
00494     }
00495 
00497     GlobalCoordinate corner ( const int i ) const
00498     {
00499       assert( (i >= 0) && (i < corners()) );
00500       const Alberta::GlobalCoordinate &x = elementInfo_.coordinate( i );
00501       GlobalCoordinate y;
00502       for( int j = 0; j < coorddimension; ++j )
00503         y[ j ] = x[ j ];
00504       return y;
00505     }
00506 
00508     const GlobalCoordinate &operator[] ( const int i ) const
00509     {
00510       return reinterpret_cast< const GlobalCoordinate & >( elementInfo_.coordinate( i ) );
00511     }
00512 
00514     GlobalCoordinate center () const
00515     {
00516       GlobalCoordinate centroid_ = corner( 0 );
00517       for( int i = 1; i < numCorners; ++i )
00518         centroid_ += corner( i );
00519       centroid_ *= ctype( 1 ) / ctype( numCorners );
00520       return centroid_;
00521     }
00522 
00524     GlobalCoordinate global ( const LocalCoordinate &local ) const;
00525   
00527     LocalCoordinate local ( const GlobalCoordinate &global ) const;
00528 
00534     ctype integrationElement () const
00535     {
00536       return elementInfo_.geometryCache().integrationElement();
00537     }
00538 
00540     ctype integrationElement ( const LocalCoordinate &local ) const
00541     {
00542       return integrationElement();
00543     }
00544 
00546     ctype volume () const
00547     {
00548       return integrationElement() / ctype( Factorial< mydimension >::factorial );
00549     }
00550 
00556     const JacobianTransposed &jacobianTransposed () const
00557     {
00558       return elementInfo_.geometryCache().jacobianTransposed();
00559     }
00560 
00562     const JacobianTransposed &
00563     jacobianTransposed ( const LocalCoordinate &local ) const
00564     {
00565       return jacobianTransposed();
00566     }
00567  
00573     const JacobianInverseTransposed &jacobianInverseTransposed () const
00574     {
00575       return elementInfo_.geometryCache().jacobianInverseTransposed();
00576     }
00577 
00579     const JacobianInverseTransposed &
00580     jacobianInverseTransposed ( const LocalCoordinate &local ) const
00581     {
00582       return jacobianInverseTransposed();
00583     }
00584 
00585     //***********************************************************************
00586     //  Methods that not belong to the Interface, but have to be public
00587     //***********************************************************************
00588 
00590     void invalidate ()
00591     {
00592       elementInfo_ = ElementInfo(); 
00593     }
00594 
00596     template< class CoordReader >
00597     void build ( const CoordReader &coordReader )
00598     {
00599       elementInfo_ = coordReader.elementInfo();
00600     }
00601 
00602   private:
00603     ElementInfo elementInfo_;
00604   };
00605 #endif // #if !DUNE_ALBERTA_CACHE_COORDINATES
00606 #endif // #if !DUNE_ALBERTA_USE_GENERICGEOMETRY
00607 
00608 
00609 
00610   // AlbertaGridLocalGeometryProvider
00611   // --------------------------------
00612 
00613   template< class Grid >
00614   class AlbertaGridLocalGeometryProvider
00615   {
00616     typedef AlbertaGridLocalGeometryProvider< Grid > This;
00617 
00618   public:
00619     typedef typename Grid::ctype ctype;
00620     
00621     static const int dimension = Grid::dimension;
00622 
00623     template< int codim >
00624     struct Codim
00625     {
00626       typedef AlbertaGridGeometry< dimension-codim, dimension, Grid > LocalGeometry;
00627     };
00628 
00629     typedef typename Codim< 0 >::LocalGeometry LocalElementGeometry;
00630     typedef typename Codim< 1 >::LocalGeometry LocalFaceGeometry;
00631 
00632     static const int numChildren = 2;
00633     static const int numFaces = dimension + 1;
00634 
00635     static const int minFaceTwist = Alberta::Twist< dimension, dimension-1 >::minTwist;
00636     static const int maxFaceTwist = Alberta::Twist< dimension, dimension-1 >::maxTwist;
00637     static const int numFaceTwists = maxFaceTwist - minFaceTwist + 1;
00638 
00639   private:
00640     struct GeoInFatherCoordReader;
00641     struct FaceCoordReader;
00642 
00643     AlbertaGridLocalGeometryProvider ()
00644     {
00645       buildGeometryInFather();
00646       buildFaceGeometry();
00647     }
00648 
00649     ~AlbertaGridLocalGeometryProvider ()
00650     {
00651       for( int child = 0; child < numChildren; ++child )
00652       {
00653         delete geometryInFather_[ child ][ 0 ];
00654         delete geometryInFather_[ child ][ 1 ];
00655       }
00656 
00657       for( int i = 0; i < numFaces; ++i )
00658       {
00659         for( int j = 0; j < numFaceTwists; ++j )
00660           delete faceGeometry_[ i ][ j ];
00661       }
00662     }
00663 
00664     void buildGeometryInFather();
00665     void buildFaceGeometry();
00666 
00667   public:
00668     const LocalElementGeometry &
00669     geometryInFather ( int child, const int orientation = 1 ) const
00670     {
00671       assert( (child >= 0) && (child < numChildren) );
00672       assert( (orientation == 1) || (orientation == -1) );
00673       return *geometryInFather_[ child ][ (orientation + 1) / 2 ];
00674     }
00675 
00676     const LocalFaceGeometry &
00677     faceGeometry ( int face, int twist = 0 ) const
00678     {
00679       assert( (face >= 0) && (face < numFaces) );
00680       assert( (twist >= minFaceTwist) && (twist <= maxFaceTwist) );
00681       return *faceGeometry_[ face ][ twist - minFaceTwist ];
00682     }
00683 
00684     static const This &instance ()
00685     {
00686       static This theInstance;
00687       return theInstance;
00688     }
00689 
00690   private:
00691     template< int codim >
00692     static int mapVertices ( int subEntity, int i )
00693     {
00694       return Alberta::MapVertices< dimension, codim >::apply( subEntity, i );
00695     }
00696 
00697     const LocalElementGeometry *geometryInFather_[ numChildren ][ 2 ];
00698     const LocalFaceGeometry *faceGeometry_[ numFaces ][ numFaceTwists ];
00699   };
00700 
00701 
00702 
00703   // FacadeOptions
00704   // -------------
00705 
00706   namespace FacadeOptions
00707   {
00708 
00709     template< int mydim, int cdim, class Grid >
00710     struct StoreGeometryReference< mydim, cdim, Grid, AlbertaGridGlobalGeometry >
00711     {
00712       static const bool v = false;
00713     };
00714 
00715   } // namespace FacadeOptions
00716 
00717 } // namespace Dune
00718 
00719 #endif // #if HAVE_ALBERTA
00720 
00721 #endif // #ifndef DUNE_ALBERTA_GEOMETRY_HH