dune-grid  2.2.0
alugrid/2d/geometry.hh
Go to the documentation of this file.
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