dune-geometry  2.2.0
cornermapping.hh
Go to the documentation of this file.
00001 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH
00002 #define DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH
00003 
00004 #include <dune/geometry/genericgeometry/topologytypes.hh>
00005 #include <dune/geometry/genericgeometry/referenceelements.hh>
00006 #include <dune/geometry/genericgeometry/matrixhelper.hh>
00007 
00008 namespace Dune
00009 {
00010 
00011   namespace GenericGeometry
00012   {
00013 
00014     // External Forward Declarations
00015     // -----------------------------
00016 
00017     template< class CT, unsigned int dim, unsigned int dimW >
00018     struct MappingTraits;
00019 
00020 
00021 
00022     // GenericCornerMapping
00023     // --------------------
00024 
00025     template< class Topology, class Traits, bool affine, unsigned int offset = 0 >
00026     class GenericCornerMapping;
00027 
00028     template< class Traits, bool affine, unsigned int offset >
00029     class GenericCornerMapping < Point, Traits, affine, offset > 
00030     {
00031       typedef Point Topology;
00032      
00033     public:
00034       static const unsigned int dim = Topology::dimension;
00035       static const unsigned int dimW = Traits::dimWorld;
00036 
00037       typedef typename Traits::FieldType FieldType;           
00038       typedef typename Traits::LocalCoordinate LocalCoordinate;
00039       typedef typename Traits::GlobalCoordinate GlobalCoordinate;
00040       typedef typename Traits::JacobianTransposedType JacobianTransposedType;
00041 
00042       static const bool alwaysAffine = true;
00043 
00044       template< class CoordStorage >
00045       static const GlobalCoordinate &origin ( const CoordStorage &coords )
00046       {
00047         dune_static_assert( CoordStorage::size, "Invalid offset." );
00048         return coords[ offset ];
00049       }
00050 
00051       template< class CoordStorage >
00052       static void phi_set ( const CoordStorage &coords,
00053                             const LocalCoordinate &x,
00054                             const FieldType &factor,
00055                             GlobalCoordinate &p )
00056       {
00057         const GlobalCoordinate &y = origin( coords );
00058         for( unsigned int i = 0; i < dimW; ++i )
00059           p[ i ] = factor * y[ i ];
00060       }
00061       
00062       template< class CoordStorage >
00063       static void phi_add ( const CoordStorage &coords,
00064                             const LocalCoordinate &x,
00065                             const FieldType &factor,
00066                             GlobalCoordinate &p )
00067       {
00068         const GlobalCoordinate &y = origin( coords );
00069         for( unsigned int i = 0; i < dimW; ++i )
00070           p[ i ] += factor * y[ i ];
00071       }
00072 
00073       template< class CoordStorage >
00074       static bool Dphi_set ( const CoordStorage &coords,
00075                              const LocalCoordinate &x,
00076                              const FieldType &factor,
00077                              JacobianTransposedType &J )
00078       {
00079         return true;
00080       }
00081       
00082       template< class CoordStorage >
00083       static bool Dphi_add ( const CoordStorage &coords,
00084                              const LocalCoordinate &x,
00085                              const FieldType &factor,
00086                              JacobianTransposedType &J )
00087       {
00088         return true;
00089       }
00090     };
00091 
00092 
00093     template< class BaseTopology, class Traits, bool affine, unsigned int offset >
00094     class GenericCornerMapping< Prism< BaseTopology >, Traits, affine, offset > 
00095     {
00096       typedef Prism< BaseTopology > Topology;
00097       
00098       typedef GenericCornerMapping< BaseTopology, Traits, affine, offset >
00099         BottomMapping;
00100       typedef GenericCornerMapping
00101         < BaseTopology, Traits, affine, offset + BaseTopology::numCorners >
00102         TopMapping;
00103 
00104     public:
00105       static const unsigned int dim = Topology::dimension;
00106       static const unsigned int dimW = Traits::dimWorld;
00107 
00108       typedef typename Traits::FieldType FieldType;           
00109       typedef typename Traits::LocalCoordinate LocalCoordinate;
00110       typedef typename Traits::GlobalCoordinate GlobalCoordinate;
00111       typedef typename Traits::JacobianTransposedType JacobianTransposedType;
00112 
00113       static const bool alwaysAffine = ((dim < 2) || affine);
00114 
00115       template< class CoordStorage >
00116       static const GlobalCoordinate &origin ( const CoordStorage &coords )
00117       {
00118         return BottomMapping::origin( coords );
00119       }
00120 
00121       template< class CoordStorage >
00122       static void phi_set ( const CoordStorage &coords,
00123                             const LocalCoordinate &x,
00124                             const FieldType &factor,
00125                             GlobalCoordinate &p )
00126       {
00127         const FieldType xn = x[ dim-1 ];
00128         const FieldType cxn = FieldType( 1 ) - xn;
00129         BottomMapping::phi_set( coords, x, factor * cxn, p );
00130         TopMapping::phi_add( coords, x, factor * xn, p );
00131       }
00132       
00133       template< class CoordStorage >
00134       static void phi_add ( const CoordStorage &coords,
00135                             const LocalCoordinate &x,
00136                             const FieldType &factor,
00137                             GlobalCoordinate &p )
00138       {
00139         const FieldType xn = x[ dim-1 ];
00140         const FieldType cxn = FieldType( 1 ) - xn;
00141         BottomMapping::phi_add( coords, x, factor * cxn, p );
00142         TopMapping::phi_add( coords, x, factor * xn, p );
00143       }
00144     
00145       template< class CoordStorage >
00146       static bool Dphi_set ( const CoordStorage &coords,
00147                              const LocalCoordinate &x,
00148                              const FieldType &factor,
00149                              JacobianTransposedType &J )
00150       {
00151         const FieldType xn = x[ dim-1 ];
00152         bool isAffine = true;
00153         if( alwaysAffine )
00154         {
00155           const FieldType cxn = FieldType( 1 ) - xn;
00156           BottomMapping::Dphi_set( coords, x, factor * cxn, J );
00157           TopMapping::Dphi_add( coords, x, factor * xn, J );
00158         }
00159         else
00160         {
00161           JacobianTransposedType Jtop;
00162           isAffine &= BottomMapping::Dphi_set( coords, x, factor, J );
00163           isAffine &= TopMapping::Dphi_set( coords, x, factor, Jtop );
00164         
00165           FieldType norm = FieldType( 0 );
00166           for( unsigned int i = 0; i < dim-1; ++i )
00167           {
00168             Jtop[ i ] -= J[ i ];
00169             norm += Jtop[ i ].two_norm2();
00170             J[ i ].axpy( xn, Jtop[ i ] );
00171           }
00172           isAffine &= (norm < 1e-12);
00173         }
00174         BottomMapping::phi_set( coords, x, -factor, J[ dim-1 ] );
00175         TopMapping::phi_add( coords, x, factor, J[ dim-1 ] );
00176         return isAffine;
00177       }
00178       
00179       template< class CoordStorage >
00180       static bool Dphi_add ( const CoordStorage &coords,
00181                              const LocalCoordinate &x,
00182                              const FieldType &factor,
00183                              JacobianTransposedType &J )
00184       {
00185         const FieldType xn = x[ dim-1 ];
00186         bool isAffine = true;
00187         if( alwaysAffine )
00188         {
00189           const FieldType cxn = FieldType( 1 ) - xn;
00190           BottomMapping::Dphi_add( coords, x, factor * cxn, J );
00191           TopMapping::Dphi_add( coords, x, factor * xn, J );
00192         }
00193         else
00194         {
00195           JacobianTransposedType Jbottom, Jtop;
00196           isAffine &= BottomMapping::Dphi_set( coords, x, FieldType( 1 ), Jbottom );
00197           isAffine &= TopMapping::Dphi_set( coords, x, FieldType( 1 ), Jtop );
00198 
00199           FieldType norm = FieldType( 0 );
00200           for( unsigned int i = 0; i < dim-1; ++i )
00201           {
00202             Jtop[ i ] -= Jbottom[ i ];
00203             norm += Jtop[ i ].two_norm2();
00204             J[ i ].axpy( factor, Jbottom[ i ] );
00205             J[ i ].axpy( factor*xn, Jtop[ i ] );
00206           }
00207           isAffine &= (norm < 1e-12);
00208         }
00209         BottomMapping::phi_add( coords, x, -factor, J[ dim-1 ] );
00210         TopMapping::phi_add( coords, x, factor, J[ dim-1 ] );
00211         return isAffine;
00212       }
00213     };
00214     
00215 
00216     template< class BaseTopology, class Traits, bool affine, unsigned int offset >
00217     class GenericCornerMapping < Pyramid< BaseTopology >, Traits, affine, offset >
00218     {
00219       typedef Pyramid< BaseTopology > Topology;
00220       
00221       typedef GenericCornerMapping< BaseTopology, Traits, affine, offset >
00222         BottomMapping;
00223       typedef GenericCornerMapping
00224         < Point, Traits, affine, offset + BaseTopology::numCorners >
00225         TopMapping;
00226 
00227     public:
00228       static const unsigned int dim = Topology::dimension;
00229       static const unsigned int dimW = Traits::dimWorld;
00230 
00231       typedef typename Traits::FieldType FieldType;           
00232       typedef typename Traits::LocalCoordinate LocalCoordinate;
00233       typedef typename Traits::GlobalCoordinate GlobalCoordinate;
00234       typedef typename Traits::JacobianTransposedType JacobianTransposedType;
00235 
00236       static const bool alwaysAffine = (BottomMapping::alwaysAffine || affine);
00237     
00238       template< class CoordStorage >
00239       static const GlobalCoordinate &origin ( const CoordStorage &coords )
00240       {
00241         return BottomMapping::origin( coords );
00242       }
00243 
00244       template< class CoordStorage >
00245       static void phi_set ( const CoordStorage &coords,
00246                             const LocalCoordinate &x,
00247                             const FieldType &factor,
00248                             GlobalCoordinate &p )
00249       {
00250         const FieldType xn = x[ dim-1 ];
00251         if( alwaysAffine )
00252         {
00253           const GlobalCoordinate &top = TopMapping::origin( coords );
00254           const GlobalCoordinate &bottom = BottomMapping::origin( coords );
00255           
00256           BottomMapping::phi_set( coords, x, factor, p );
00257           for( unsigned int i = 0; i < dimW; ++i )
00258             p[ i ] += (factor * xn) * (top[ i ] - bottom[ i ]);
00259         }
00260         else
00261         {
00262           TopMapping::phi_set( coords, x, factor * xn, p );
00263           const FieldType cxn = FieldType( 1 ) - xn;
00264           if( cxn > 1e-12 )
00265           {
00266             const FieldType icxn = FieldType( 1 ) / cxn;
00267             LocalCoordinate xb;
00268             for( unsigned int i = 0; i < dim-1; ++i )
00269               xb[ i ] = icxn * x[ i ];
00270 
00271             BottomMapping::phi_add( coords, xb, factor * cxn, p );
00272           }
00273         }
00274       }
00275       
00276       template< class CoordStorage >
00277       static void phi_add ( const CoordStorage &coords,
00278                             const LocalCoordinate &x,
00279                             const FieldType &factor,
00280                             GlobalCoordinate &p )
00281       {
00282         const FieldType xn = x[ dim-1 ];
00283         if( alwaysAffine )
00284         {
00285           const GlobalCoordinate &top = TopMapping::origin( coords );
00286           const GlobalCoordinate &bottom = BottomMapping::origin( coords );
00287           
00288           BottomMapping::phi_add( coords, x, factor, p );
00289           for( unsigned int i = 0; i < dimW; ++i )
00290             p[ i ] += (factor * xn) * (top[ i ] - bottom[ i ]);
00291         }
00292         else
00293         {
00294           TopMapping::phi_add( coords, x, factor * xn, p );
00295           const FieldType cxn = FieldType( 1 ) - xn;
00296           if( cxn > 1e-12 )
00297           {
00298             const FieldType icxn = FieldType( 1 ) / cxn;
00299             LocalCoordinate xb;
00300             for( unsigned int i = 0; i < dim-1; ++i )
00301               xb[ i ] = icxn * x[ i ];
00302 
00303             BottomMapping::phi_add( coords, xb, factor * cxn, p );
00304           }
00305         }
00306       }
00307 
00308       template< class CoordStorage >
00309       static bool Dphi_set ( const CoordStorage &coords,
00310                              const LocalCoordinate &x,
00311                              const FieldType &factor,
00312                              JacobianTransposedType &J )
00313       {
00314         GlobalCoordinate &q = J[ dim-1 ];
00315         bool isAffine;
00316         if( alwaysAffine )
00317         {
00318           const GlobalCoordinate &top = TopMapping::origin( coords );
00319           const GlobalCoordinate &bottom = BottomMapping::origin( coords );
00320 
00321           isAffine = BottomMapping::Dphi_set( coords, x, factor, J );
00322           for( unsigned int i = 0; i < dimW; ++i )
00323             q[ i ] = factor * (top[ i ] - bottom[ i ]);
00324         }
00325         else
00326         {
00327           const FieldType xn = x[ dim-1 ];
00328           const FieldType cxn = FieldType( 1 ) - xn;
00329           const FieldType icxn = FieldType( 1 ) / cxn;
00330           LocalCoordinate xb;
00331           for( unsigned int i = 0; i < dim-1; ++i )
00332             xb[ i ] = icxn * x[ i ];
00333           isAffine = BottomMapping::Dphi_set( coords, xb, factor, J );
00334           
00335           TopMapping::phi_set( coords, x, factor, q );
00336           BottomMapping::phi_add( coords, xb, -factor, q );
00337           xb *= factor;
00338           for( unsigned int j = 0; j < dim-1; ++j )
00339           {
00340             for( unsigned int i = 0; i < dimW; ++i )
00341               q[ i ] += J[ j ][ i ] * xb[ j ];
00342           }
00343         }
00344         return isAffine;
00345       }
00346       
00347       template< class CoordStorage >
00348       static bool Dphi_add ( const CoordStorage &coords,
00349                              const LocalCoordinate &x,
00350                              const FieldType &factor,
00351                              JacobianTransposedType &J )
00352       {
00353         GlobalCoordinate &q = J[ dim-1 ];
00354         bool isAffine;
00355         if( alwaysAffine )
00356         {
00357           const GlobalCoordinate &top = TopMapping::origin( coords );
00358           const GlobalCoordinate &bottom = BottomMapping::origin( coords );
00359 
00360           isAffine = BottomMapping::Dphi_add( coords, x, factor, J );
00361           for( unsigned int i = 0; i < dimW; ++i )
00362             q[ i ] += factor * (top[ i ] - bottom[ i ]);
00363         }
00364         else
00365         {
00366           const FieldType xn = x[ dim-1 ];
00367           const FieldType cxn = FieldType( 1 ) - xn;
00368           const FieldType icxn = FieldType( 1 ) / cxn;
00369           LocalCoordinate xb;
00370           for( unsigned int i = 0; i < dim-1; ++i )
00371             xb[ i ] = icxn * x[ i ];
00372           isAffine = BottomMapping::Dphi_add( coords, xb, factor, J );
00373           
00374           TopMapping::phi_add( coords, x, factor, q );
00375           BottomMapping::phi_add( coords, xb, -factor, q );
00376           xb *= factor;
00377           for( unsigned int j = 0; j < dim-1; ++j )
00378           {
00379             for( unsigned int i = 0; i < dimW; ++i )
00380               q[ i ] += J[ j ][ i ] * xb[ j ];
00381           }
00382         }
00383         return isAffine;
00384       }
00385     };
00386 
00387 
00388 
00389     // SubMappingCoords
00390     // ----------------
00391     
00392     template< class Mapping, unsigned int codim >
00393     class SubMappingCoords
00394     {
00395       typedef typename Mapping::GlobalCoordinate GlobalCoordinate;
00396       typedef typename Mapping::ReferenceElement ReferenceElement;
00397 
00398       static const unsigned int dimension = ReferenceElement::dimension;
00399 
00400       const Mapping &mapping_;
00401       const unsigned int i_;
00402 
00403     public:
00404       SubMappingCoords ( const Mapping &mapping, unsigned int i )
00405       : mapping_( mapping ), i_( i )
00406       {}
00407 
00408       const GlobalCoordinate &operator[] ( unsigned int j ) const
00409       {
00410         const unsigned int k
00411           = ReferenceElement::template subNumbering< codim, dimension - codim >( i_, j );
00412         return mapping_.corner( k );
00413       }
00414     };
00415 
00416 
00417 
00418     // CoordStorage
00419     // ------------
00420 
00425     template< class CoordTraits, class Topology, unsigned int dimW >
00426     class CoordStorage
00427     {
00428       typedef CoordStorage< CoordTraits, Topology, dimW > This;
00429 
00430     public:
00431       static const unsigned int size = Topology::numCorners;
00432 
00433       static const unsigned int dimWorld = dimW;
00434 
00435       typedef typename CoordTraits::template Vector< dimWorld >::type GlobalCoordinate;
00436 
00437       template< class SubTopology >
00438       struct SubStorage
00439       {
00440         typedef CoordStorage< CoordTraits, SubTopology, dimWorld > type;
00441       };
00442 
00443       template< class CoordVector >
00444       explicit CoordStorage ( const CoordVector &coords )
00445       {
00446         for( unsigned int i = 0; i < size; ++i )
00447           coords_[ i ] = coords[ i ];
00448       }
00449 
00450       const GlobalCoordinate &operator[] ( unsigned int i ) const
00451       {
00452         return coords_[ i ];
00453       }
00454 
00455     private:
00456       GlobalCoordinate coords_[ size ];
00457     };
00458 
00459 
00460 
00461     // CoordPointerStorage
00462     // -------------------
00463 
00468     template< class CoordTraits, class Topology, unsigned int dimW >
00469     class CoordPointerStorage
00470     {
00471       typedef CoordPointerStorage< CoordTraits, Topology, dimW > This;
00472 
00473     public:
00474       static const unsigned int size = Topology::numCorners;
00475 
00476       static const unsigned int dimWorld = dimW;
00477 
00478       typedef typename CoordTraits::template Vector< dimWorld >::type GlobalCoordinate;
00479 
00480       template< class SubTopology >
00481       struct SubStorage
00482       {
00483         typedef CoordPointerStorage< CoordTraits, SubTopology, dimWorld > type;
00484       };
00485 
00486       template< class CoordVector >
00487       explicit CoordPointerStorage ( const CoordVector &coords )
00488       {
00489         for( unsigned int i = 0; i < size; ++i )
00490           coords_[ i ] = &(coords[ i ]);
00491       }
00492 
00493       const GlobalCoordinate &operator[] ( unsigned int i ) const
00494       {
00495         return *(coords_[ i ]);
00496       }
00497 
00498     private:
00499       const GlobalCoordinate *coords_[ size ];
00500     };
00501 
00502 
00503     
00504     // CornerMapping
00505     // -------------
00506 
00512     template< class CoordTraits, class Topo, unsigned int dimW,
00513               class CStorage = CoordStorage< CoordTraits, Topo, dimW >,
00514               bool affine = false >
00515     class CornerMapping
00516     {
00517       typedef CornerMapping< CoordTraits, Topo, dimW, CStorage, affine > This;
00518       
00519     public:
00520       typedef Topo Topology;
00521       typedef CStorage CornerStorage;
00522       typedef MappingTraits< CoordTraits, Topology::dimension, dimW > Traits;
00523 
00524       static const unsigned int dimension = Traits::dimension;
00525       static const unsigned int dimWorld = Traits::dimWorld;
00526 
00527       typedef typename Traits::FieldType FieldType;           
00528       typedef typename Traits::LocalCoordinate LocalCoordinate;
00529       typedef typename Traits::GlobalCoordinate GlobalCoordinate;
00530       typedef typename Traits::JacobianType JacobianType;
00531       typedef typename Traits::JacobianTransposedType JacobianTransposedType;
00532 
00533       typedef GenericGeometry::ReferenceElement< Topology, FieldType > ReferenceElement;
00534 
00535       template< unsigned int codim, unsigned int i >
00536       struct SubTopology
00537       {
00538         typedef typename GenericGeometry::SubTopology< Topo, codim, i >::type Topology;
00539         typedef typename CStorage::template SubStorage< Topology >::type CornerStorage;
00540         typedef CornerMapping< CoordTraits, Topology, dimWorld, CornerStorage, affine > Trace;
00541       };
00542 
00543     private:
00544       typedef GenericGeometry::GenericCornerMapping< Topology, Traits, affine > GenericMapping;
00545 
00546     public:
00547       static const bool alwaysAffine = GenericMapping::alwaysAffine;
00548 
00549       template< class CoordVector >
00550       explicit CornerMapping ( const CoordVector &coords )
00551       : coords_( coords )
00552       {}
00553 
00554       const GlobalCoordinate &corner ( int i ) const
00555       {
00556         return coords_[ i ];
00557       }
00558 
00559       void global ( const LocalCoordinate &x, GlobalCoordinate &y ) const
00560       {
00561         GenericMapping::phi_set( coords_, x, FieldType( 1 ), y );
00562       }
00563 
00564       bool jacobianTransposed ( const LocalCoordinate &x, 
00565                                 JacobianTransposedType &JT ) const
00566       {
00567         return GenericMapping::Dphi_set( coords_, x, FieldType( 1 ), JT );
00568       }
00569 
00570       template< unsigned int codim, unsigned int i >
00571       typename SubTopology< codim, i >::Trace trace () const
00572       {
00573         typedef typename SubTopology< codim, i >::Trace Trace;
00574         typedef SubMappingCoords< This, codim > CoordVector;
00575         return Trace( CoordVector( *this, i ) );
00576       }
00577 
00578     protected:
00579       CornerStorage coords_;
00580     };
00581 
00582   } // namespace GenericGeometry
00583 
00584 } // namespace Dune
00585 
00586 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CORNERMAPPING_HH