dune-grid
2.2.0
|
00001 #ifndef DUNE_ALU2DGRIDGEOMETRY_HH 00002 #define DUNE_ALU2DGRIDGEOMETRY_HH 00003 00004 // Dune includes 00005 #include <dune/common/misc.hh> 00006 #include <dune/grid/common/grid.hh> 00007 #include <dune/geometry/genericgeometry/topologytypes.hh> 00008 00009 #include <dune/grid/alugrid/2d/alu2dinclude.hh> 00010 #include <dune/grid/alugrid/3d/mappings.hh> 00011 00012 namespace Dune 00013 { 00014 00015 // Forward declarations 00016 template<int cd, int dim, class GridImp> 00017 class ALU2dGridEntity; 00018 template<int cd, class GridImp > 00019 class ALU2dGridEntityPointer; 00020 template<int mydim, int cdim, class GridImp> 00021 class ALU2dGridGeometry; 00022 template< int dim, int dimworld, ALU2DSPACE ElementType eltype > 00023 class ALU2dGrid; 00024 00025 00026 template< int mydim, int cdim, ALU2DSPACE ElementType eltype > 00027 class MyALU2dGridGeometryImpl; 00028 00029 // geometry implementation for vertices 00030 template< int cdim, ALU2DSPACE ElementType eltype > 00031 class MyALU2dGridGeometryImpl< 0, cdim, eltype > 00032 { 00033 typedef LinearMapping< cdim, 0 > MappingType; 00034 00035 typedef typename MappingType::ctype ctype; 00036 00037 typedef typename MappingType::map_t map_t; 00038 typedef typename MappingType::world_t world_t; 00039 00040 typedef typename MappingType::matrix_t matrix_t; 00041 typedef typename MappingType::inv_t inv_t; 00042 00043 MappingType mapping_; 00044 bool valid_ ; 00045 00046 const MappingType& mapping() const 00047 { 00048 assert( valid() ); 00049 return mapping_; 00050 } 00051 00052 public: 00053 MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {} 00054 00055 // returns true if goemetry info is valid 00056 bool valid () const { return valid_; } 00057 00058 // reset geometry status 00059 void invalidate() { valid_ = false ; } 00060 00061 bool affine() const 00062 { 00063 assert( valid() ); 00064 return mapping().affine(); 00065 } 00066 00067 int corners () const 00068 { 00069 return 1; 00070 } 00071 00072 GeometryType type () const 00073 { 00074 return GeometryType( 00075 (eltype == ALU2DSPACE triangle ? 00076 GenericGeometry :: SimplexTopology< 0 > :: type :: id : 00077 GenericGeometry :: CubeTopology < 0 > :: type :: id), 00078 0 ); 00079 } 00080 00081 void map2world ( const map_t &m, world_t &w ) const 00082 { 00083 return mapping().map2world( m, w ); 00084 } 00085 00086 void world2map ( const world_t &w, map_t &m ) const 00087 { 00088 return mapping().world2map( w, m ); 00089 } 00090 00091 const matrix_t &jacobianTransposed ( const map_t &m ) const 00092 { 00093 return mapping().jacobianTransposed( m ); 00094 } 00095 00096 const inv_t &jacobianInverseTransposed ( const map_t &m ) const 00097 { 00098 return mapping().jacobianInverseTransposed( m ); 00099 } 00100 00101 ctype det ( const map_t &m ) const 00102 { 00103 return mapping().det( m ); 00104 } 00105 00106 // update geometry coordinates 00107 template< class Vector > 00108 void update ( const Vector &p0 ) 00109 { 00110 mapping_.buildMapping( p0 ); 00111 valid_ = true ; 00112 } 00113 }; 00114 00115 // geometry implementation for lines 00116 template< int cdim, ALU2DSPACE ElementType eltype > 00117 class MyALU2dGridGeometryImpl< 1, cdim, eltype > 00118 { 00119 static const int ncorners = 2; 00120 00121 typedef LinearMapping< cdim, 1 > MappingType; 00122 00123 typedef typename MappingType::ctype ctype; 00124 00125 typedef typename MappingType::map_t map_t; 00126 typedef typename MappingType::world_t world_t; 00127 00128 typedef typename MappingType::matrix_t matrix_t; 00129 typedef typename MappingType::inv_t inv_t; 00130 00131 MappingType mapping_; 00132 bool valid_; 00133 00134 const MappingType& mapping() const 00135 { 00136 assert( valid() ); 00137 return mapping_; 00138 } 00139 00140 public: 00141 MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {} 00142 00143 // returns true if goemetry info is valid 00144 bool valid () const { return valid_; } 00145 00146 // reset geometry status 00147 void invalidate() { valid_ = false ; } 00148 00149 bool affine() const 00150 { 00151 return mapping().affine(); 00152 } 00153 00154 int corners () const 00155 { 00156 return ncorners; 00157 } 00158 00159 GeometryType type () const 00160 { 00161 return GeometryType( 00162 (eltype == ALU2DSPACE triangle ? 00163 GenericGeometry :: SimplexTopology< 1 > :: type :: id : 00164 GenericGeometry :: CubeTopology < 1 > :: type :: id), 00165 1 ); 00166 } 00167 00168 void map2world ( const map_t &m, world_t &w ) const 00169 { 00170 return mapping().map2world( m, w ); 00171 } 00172 00173 void world2map ( const world_t &w, map_t &m ) const 00174 { 00175 return mapping().world2map( w, m ); 00176 } 00177 00178 const matrix_t &jacobianTransposed ( const map_t &m ) const 00179 { 00180 return mapping().jacobianTransposed( m ); 00181 } 00182 00183 const inv_t &jacobianInverseTransposed ( const map_t &m ) const 00184 { 00185 return mapping().jacobianInverseTransposed( m ); 00186 } 00187 00188 ctype det ( const map_t &m ) const 00189 { 00190 return mapping().det( m ); 00191 } 00192 00193 // update geometry in father coordinates 00194 template< class Geo, class LocalGeo > 00195 void updateLocal ( const Geo &geo, const LocalGeo &localGeo ) 00196 { 00197 assert( localGeo.corners() == ncorners ); 00198 // compute the local coordinates in father refelem 00199 FieldMatrix< alu2d_ctype, ncorners, cdim > coord; 00200 for( int i = 0; i < ncorners; ++i ) 00201 { 00202 // calculate coordinate 00203 coord[ i ] = geo.local( localGeo.corner( i ) ); 00204 // to avoid rounding errors 00205 for( int j = 0; j < cdim; ++j ) 00206 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]); 00207 } 00208 mapping_.buildMapping( coord[ 0 ], coord[ 1 ] ); 00209 valid_ = true ; 00210 } 00211 00212 // update geometry coordinates 00213 template< class Vector > 00214 void update ( const Vector &p0, const Vector &p1 ) 00215 { 00216 mapping_.buildMapping( p0, p1 ); 00217 valid_ = true ; 00218 } 00219 }; 00220 00221 // geometry implementation for triangles 00222 template< int cdim > 00223 class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE triangle > 00224 { 00225 static const int ncorners = 3; 00226 00227 typedef LinearMapping< cdim, 2 > MappingType; 00228 00229 typedef typename MappingType::ctype ctype; 00230 00231 typedef typename MappingType::map_t map_t; 00232 typedef typename MappingType::world_t world_t; 00233 00234 typedef typename MappingType::matrix_t matrix_t; 00235 typedef typename MappingType::inv_t inv_t; 00236 00237 MappingType mapping_; 00238 bool valid_; 00239 00240 const MappingType& mapping() const 00241 { 00242 assert( valid() ); 00243 return mapping_; 00244 } 00245 00246 public: 00247 MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {} 00248 00249 // returns true if goemetry info is valid 00250 bool valid () const { return valid_; } 00251 00252 // reset geometry status 00253 void invalidate() { valid_ = false ; } 00254 00255 bool affine () const 00256 { 00257 return mapping().affine(); 00258 } 00259 00260 int corners () const 00261 { 00262 return ncorners; 00263 } 00264 00265 GeometryType type () const 00266 { 00267 return GeometryType( GenericGeometry :: SimplexTopology< 2 > :: type :: id , 2 ); 00268 } 00269 00270 void map2world ( const map_t &m, world_t &w ) const 00271 { 00272 return mapping().map2world( m, w ); 00273 } 00274 00275 void world2map ( const world_t &w, map_t &m ) const 00276 { 00277 return mapping().world2map( w, m ); 00278 } 00279 00280 const matrix_t &jacobianTransposed ( const map_t &m ) const 00281 { 00282 return mapping().jacobianTransposed( m ); 00283 } 00284 00285 const inv_t &jacobianInverseTransposed ( const map_t &m ) const 00286 { 00287 return mapping().jacobianInverseTransposed( m ); 00288 } 00289 00290 ctype det ( const map_t &m ) const 00291 { 00292 return mapping().det( m ); 00293 } 00294 00295 // update geometry in father coordinates 00296 template< class Geo, class LocalGeo > 00297 void updateLocal ( const Geo &geo, const LocalGeo &localGeo ) 00298 { 00299 assert( localGeo.corners() == ncorners ); 00300 // compute the local coordinates in father refelem 00301 FieldMatrix< alu2d_ctype, ncorners, cdim > coord; 00302 for( int i = 0; i < ncorners; ++i ) 00303 { 00304 // calculate coordinate 00305 coord[ i ] = geo.local( localGeo.corner( i ) ); 00306 // to avoid rounding errors 00307 for( int j = 0; j < cdim; ++j ) 00308 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]); 00309 } 00310 mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] ); 00311 valid_ = true ; 00312 } 00313 00314 template< class HElement > 00315 void update ( const HElement &item ) 00316 { 00317 mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(), 00318 item.getVertex( 2 )->coord() ); 00319 valid_ = true ; 00320 } 00321 }; 00322 00323 // geometry implementation for quadrilaterals 00324 template< int cdim > 00325 class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE quadrilateral > 00326 { 00327 static const int ncorners = 4; 00328 00329 typedef BilinearMapping< cdim > MappingType; 00330 00331 typedef typename MappingType::ctype ctype; 00332 00333 typedef typename MappingType::map_t map_t; 00334 typedef typename MappingType::world_t world_t; 00335 00336 typedef typename MappingType::matrix_t matrix_t; 00337 typedef typename MappingType::inv_t inv_t; 00338 00339 MappingType mapping_; 00340 bool valid_ ; 00341 00342 const MappingType& mapping() const 00343 { 00344 assert( valid() ); 00345 return mapping_; 00346 } 00347 00348 public: 00349 MyALU2dGridGeometryImpl() : mapping_(), valid_( false ) {} 00350 00351 // returns true if goemetry info is valid 00352 bool valid () const { return valid_; } 00353 00354 // reset geometry status 00355 void invalidate() { valid_ = false ; } 00356 00357 bool affine () const 00358 { 00359 return mapping().affine(); 00360 } 00361 00362 int corners () const 00363 { 00364 return ncorners; 00365 } 00366 00367 GeometryType type () const 00368 { 00369 return GeometryType( GeometryType::cube, 2 ); 00370 } 00371 00372 void map2world ( const map_t &m, world_t &w ) const 00373 { 00374 return mapping().map2world( m, w ); 00375 } 00376 00377 void world2map ( const world_t &w, map_t &m ) const 00378 { 00379 return mapping().world2map( w, m ); 00380 } 00381 00382 const matrix_t &jacobianTransposed ( const map_t &m ) const 00383 { 00384 return mapping().jacobianTransposed( m ); 00385 } 00386 00387 const inv_t &jacobianInverseTransposed ( const map_t &m ) const 00388 { 00389 return mapping().jacobianInverseTransposed( m ); 00390 } 00391 00392 ctype det ( const map_t &m ) const 00393 { 00394 return mapping().det( m ); 00395 } 00396 00397 // update geometry in father coordinates 00398 template< class Geo, class LocalGeo > 00399 void updateLocal ( const Geo &geo, const LocalGeo &localGeo ) 00400 { 00401 assert( localGeo.corners() == ncorners ); 00402 // compute the local coordinates in father refelem 00403 FieldMatrix< alu2d_ctype, ncorners, cdim > coord; 00404 for( int i = 0; i < ncorners; ++i ) 00405 { 00406 // calculate coordinate 00407 coord[ i ] = geo.local( localGeo.corner( i ) ); 00408 // to avoid rounding errors 00409 for( int j = 0; j < cdim; ++j ) 00410 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]); 00411 } 00412 mapping_.buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] ); 00413 valid_ = true ; 00414 } 00415 00416 template< class HElement > 00417 void update ( const HElement &item ) 00418 { 00419 mapping_.buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(), 00420 item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() ); 00421 valid_ = true ; 00422 } 00423 }; 00424 00425 // geometry implementation for triangles 00426 template< int cdim > 00427 class MyALU2dGridGeometryImpl< 2, cdim, ALU2DSPACE mixed > 00428 { 00429 typedef Dune::LinearMapping< cdim, 2 > LinearMapping; 00430 typedef Dune::BilinearMapping< cdim > BilinearMapping; 00431 00432 typedef typename LinearMapping::ctype ctype; 00433 00434 typedef typename LinearMapping::map_t map_t; 00435 typedef typename LinearMapping::world_t world_t; 00436 00437 typedef typename LinearMapping::matrix_t matrix_t; 00438 typedef typename LinearMapping::inv_t inv_t; 00439 00440 static const int lms = sizeof( LinearMapping ); 00441 static const int bms = sizeof( BilinearMapping ); 00442 00443 char mapping_[ lms > bms ? lms : bms ]; 00444 int corners_; 00445 bool valid_ ; 00446 00447 public: 00448 MyALU2dGridGeometryImpl () : corners_( 0 ), valid_( false ) {} 00449 00450 MyALU2dGridGeometryImpl ( const MyALU2dGridGeometryImpl &other ) 00451 : corners_( other.corners() ), valid_( other.valid_ ) 00452 { 00453 if( corners_ == 3 ) 00454 new( &mapping_ ) LinearMapping( other.linearMapping() ); 00455 if( corners_ == 4 ) 00456 new( &mapping_ ) BilinearMapping( other.bilinearMapping() ); 00457 } 00458 00459 // returns true if goemetry info is valid 00460 bool valid () const { return valid_; } 00461 00462 // reset geometry status 00463 void invalidate() { valid_ = false ; } 00464 00465 bool affine () const 00466 { 00467 return (corners() == 3 ? linearMapping().affine() : bilinearMapping().affine()); 00468 } 00469 00470 int corners () const 00471 { 00472 return corners_; 00473 } 00474 00475 GeometryType type () const 00476 { 00477 return GeometryType( (corners_ == 3 ? 00478 GenericGeometry :: SimplexTopology< 2 > :: type :: id : 00479 GenericGeometry :: CubeTopology < 2 > :: type :: id), 2); 00480 } 00481 00482 void map2world ( const map_t &m, world_t &w ) const 00483 { 00484 if( corners() == 3 ) 00485 linearMapping().map2world( m, w ); 00486 else 00487 bilinearMapping().map2world( m, w ); 00488 } 00489 00490 void world2map ( const world_t &w, map_t &m ) const 00491 { 00492 if( corners() == 3 ) 00493 linearMapping().world2map( w, m ); 00494 else 00495 bilinearMapping().world2map( w, m ); 00496 } 00497 00498 const matrix_t &jacobianTransposed ( const map_t &m ) const 00499 { 00500 return (corners() == 3 ? linearMapping().jacobianTransposed( m ) : bilinearMapping().jacobianTransposed( m )); 00501 } 00502 00503 const inv_t &jacobianInverseTransposed ( const map_t &m ) const 00504 { 00505 return (corners() == 3 ? linearMapping().jacobianInverseTransposed( m ) : bilinearMapping().jacobianInverseTransposed( m )); 00506 } 00507 00508 ctype det ( const map_t &m ) const 00509 { 00510 return (corners() == 3 ? linearMapping().det( m ) : bilinearMapping().det( m )); 00511 } 00512 00513 // update geometry in father coordinates 00514 template< class Geo, class LocalGeo > 00515 void updateLocal ( const Geo &geo, const LocalGeo &localGeo ) 00516 { 00517 const int corners = localGeo.corners(); 00518 00519 // compute the local coordinates in father refelem 00520 FieldMatrix< alu2d_ctype, 4, cdim > coord; 00521 for( int i = 0; i < corners; ++i ) 00522 { 00523 // calculate coordinate 00524 coord[ i ] = geo.local( localGeo.corner( i ) ); 00525 // to avoid rounding errors 00526 for( int j = 0; j < cdim; ++j ) 00527 coord[ i ][ j ] = (coord[ i ][ j ] < 1e-14 ? 0 : coord[ i ][ j ]); 00528 } 00529 00530 updateMapping( corners ); 00531 if( corners == 3 ) 00532 linearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ] ); 00533 else 00534 bilinearMapping().buildMapping( coord[ 0 ], coord[ 1 ], coord[ 2 ], coord[ 3 ] ); 00535 00536 valid_ = true ; 00537 } 00538 00539 template< class HElement > 00540 void update ( const HElement &item ) 00541 { 00542 const int corners = item.numvertices(); 00543 updateMapping( corners ); 00544 if( corners == 3 ) 00545 linearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(), 00546 item.getVertex( 2 )->coord() ); 00547 else 00548 bilinearMapping().buildMapping( item.getVertex( 0 )->coord(), item.getVertex( 1 )->coord(), 00549 item.getVertex( 3 )->coord(), item.getVertex( 2 )->coord() ); 00550 00551 valid_ = true ; 00552 } 00553 00554 private: 00555 MyALU2dGridGeometryImpl &operator= ( const MyALU2dGridGeometryImpl &other ); 00556 00557 const LinearMapping &linearMapping () const 00558 { 00559 assert( valid() ); 00560 return static_cast< const LinearMapping * >( &mapping_ ); 00561 } 00562 00563 LinearMapping &linearMapping () 00564 { 00565 assert( valid() ); 00566 return static_cast< LinearMapping * >( &mapping_ ); 00567 } 00568 00569 const BilinearMapping &bilinearMapping () const 00570 { 00571 assert( valid() ); 00572 return static_cast< const BilinearMapping * >( &mapping_ ); 00573 } 00574 00575 BilinearMapping &bilinearMapping () 00576 { 00577 assert( valid() ); 00578 return static_cast< BilinearMapping * >( &mapping_ ); 00579 } 00580 00581 void updateMapping ( const int corners ) 00582 { 00583 assert( (corners == 3) || (corners == 4) ); 00584 if( corners != corners_ ) 00585 { 00586 destroyMapping(); 00587 corners = corners_; 00588 if( corners == 3 ) 00589 new( &mapping_ ) LinearMapping; 00590 else 00591 new( &mapping_ ) BilinearMapping; 00592 } 00593 } 00594 00595 void destroyMapping () 00596 { 00597 if( corners() == 3 ) 00598 linearMapping().~LinearMapping(); 00599 else if( corners() == 4 ) 00600 bilinearMapping().~BilinearMapping(); 00601 } 00602 }; 00603 00604 00605 //********************************************************************** 00606 // 00607 // --ALU2dGridGeometry 00608 // --Geometry 00609 //********************************************************************** 00622 00623 00624 template< int mydim, int cdim, class GridImp > 00625 class ALU2dGridGeometry 00626 : public GeometryDefaultImplementation< mydim, cdim, GridImp, ALU2dGridGeometry > 00627 { 00628 static const ALU2DSPACE ElementType eltype = GridImp::elementType; 00629 00631 typedef typename GridImp::template Codim<0>::Geometry Geometry; 00633 typedef ALU2dGridGeometry<mydim,cdim,GridImp> GeometryImp; 00635 enum { dimbary=mydim+1}; 00636 00637 typedef typename ALU2dImplTraits< GridImp::dimensionworld, eltype >::HElementType HElementType ; 00638 typedef typename ALU2dImplInterface< 0, GridImp::dimensionworld, eltype >::Type VertexType; 00639 00640 // type of specialized geometry implementation 00641 typedef MyALU2dGridGeometryImpl< mydim, cdim, eltype > GeometryImplType; 00642 00643 public: 00644 typedef FieldVector< alu2d_ctype, cdim > GlobalCoordinate; 00645 typedef FieldVector< alu2d_ctype, mydim > LocalCoordinate; 00646 00647 public: 00650 ALU2dGridGeometry(); 00651 00654 const GeometryType type () const { return geoImpl_.type(); } 00655 00657 int corners () const { return geoImpl_.corners(); } 00658 00660 const GlobalCoordinate &operator[] ( int i ) const; 00661 00663 GlobalCoordinate corner ( int i ) const; 00664 00667 GlobalCoordinate global ( const LocalCoordinate& local ) const; 00668 00671 LocalCoordinate local (const GlobalCoordinate& global) const; 00672 00674 alu2d_ctype integrationElement (const LocalCoordinate& local) const; 00675 00677 alu2d_ctype volume () const; 00678 00680 bool affine() const { return geoImpl_.affine(); } 00681 00683 const FieldMatrix<alu2d_ctype,cdim,mydim>& jacobianInverseTransposed (const LocalCoordinate& local) const; 00684 00686 const FieldMatrix<alu2d_ctype,mydim,cdim>& jacobianTransposed (const LocalCoordinate& local) const; 00687 00688 //*********************************************************************** 00690 //*********************************************************************** 00692 // method for elements 00693 bool buildGeom(const HElementType & item); 00694 // method for edges 00695 bool buildGeom(const HElementType & item, const int aluFace); 00696 // method for vertices 00697 bool buildGeom(const VertexType & item, const int ); 00698 00701 template <class GeometryType, class LocalGeomType > 00702 bool buildLocalGeom(const GeometryType & geo , const LocalGeomType & lg); 00703 00705 bool buildLocalGeometry(const int faceNumber, const int twist,const int coorns); 00706 00708 GlobalCoordinate& getCoordVec (int i); 00709 00711 void print (std::ostream& ss) const; 00712 00714 inline bool buildGeomInFather(const Geometry &fatherGeom, 00715 const Geometry & myGeom ); 00716 00717 // returns true if geometry information is valid 00718 inline bool valid() const { return geoImpl_.valid(); } 00719 00720 // invalidate geometry information 00721 inline void invalidate() const { geoImpl_.invalidate(); } 00722 00723 protected: 00724 // return reference coordinates of the alu triangle 00725 static std::pair< FieldMatrix< alu2d_ctype, 4, 2 >, FieldVector< alu2d_ctype, 4 > > 00726 calculateReferenceCoords ( const int corners ); 00727 00728 // implementation of coord and mapping 00729 mutable GeometryImplType geoImpl_; 00730 00731 // determinant 00732 mutable alu2d_ctype det_; 00733 }; 00734 00735 } // end namespace Dune 00736 00737 #include "geometry_imp.cc" 00738 00739 #endif