dune-grid
2.2.0
|
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