dune-geometry  2.2.0
cachedmapping.hh
Go to the documentation of this file.
00001 #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
00002 #define DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH
00003 
00004 #include <dune/geometry/type.hh>
00005 #include <dune/geometry/genericgeometry/topologytypes.hh>
00006 #include <dune/geometry/genericgeometry/referenceelements.hh>
00007 #include <dune/geometry/genericgeometry/matrixhelper.hh>
00008 #include <dune/geometry/genericgeometry/mapping.hh>
00009 //#include <dune/geometry/genericgeometry/traceprovider.hh>
00010 
00011 namespace Dune
00012 {
00013 
00014   namespace GenericGeometry
00015   {
00016 
00017     // Internal Forward Declarations
00018     // -----------------------------
00019 
00020     template< unsigned int, class >
00021     class CachedJacobianTransposed;
00022 
00023     template< unsigned int, class >
00024     class CachedJacobianInverseTransposed;
00025 
00026 
00027 
00028     // CachedStorage
00029     // -------------
00030 
00031     template< unsigned int dim, class GeometryTraits >
00032     class CachedStorage
00033     {
00034       friend class CachedJacobianTransposed< dim, GeometryTraits >;
00035 
00036     public:
00037       static const unsigned int dimension = dim;
00038       static const unsigned int dimWorld = GeometryTraits::dimWorld;
00039 
00040       typedef MappingTraits< typename GeometryTraits::CoordTraits, dimension, dimWorld > Traits;
00041 
00042       typedef typename GeometryTraits::Caching Caching;
00043 
00044       CachedStorage ()
00045       : affine( false ),
00046         jacobianTransposedComputed( false ),
00047         jacobianInverseTransposedComputed( false ),
00048         integrationElementComputed( false )
00049       {}
00050 
00051       typename Traits::JacobianTransposedType jacobianTransposed;
00052       typename Traits::JacobianType jacobianInverseTransposed;
00053       typename Traits::FieldType integrationElement;
00054 
00055       bool affine : 1;
00056 
00057       bool jacobianTransposedComputed : 1;        // = affine, if jacobian transposed was computed
00058       bool jacobianInverseTransposedComputed : 1; // = affine, if jacobian inverse transposed was computed
00059       bool integrationElementComputed : 1;        // = affine, if integration element was computed
00060     };
00061 
00062 
00063 
00064     // CachedJacobianTranposed
00065     // -----------------------
00066 
00067     template< unsigned int dim, class GeometryTraits >
00068     class CachedJacobianTransposed
00069     {
00070       friend class CachedJacobianInverseTransposed< dim, GeometryTraits >;
00071 
00072       typedef CachedStorage< dim, GeometryTraits > Storage;
00073       typedef typename Storage::Traits Traits;
00074 
00075       typedef typename Traits::MatrixHelper MatrixHelper;
00076 
00077     public:
00078       typedef typename Traits::FieldType ctype;
00079 
00080       static const int rows = Traits::dimension;
00081       static const int cols = Traits::dimWorld;
00082 
00083       typedef typename Traits::JacobianTransposedType FieldMatrix;
00084 
00085       operator bool () const
00086       {
00087         return storage().jacobianTransposedComputed;
00088       }
00089 
00090       operator const FieldMatrix & () const
00091       {
00092         return storage().jacobianTransposed;
00093       }
00094 
00095       template< class X, class Y >
00096       void mv ( const X &x, Y &y ) const
00097       {
00098         static_cast< const FieldMatrix & >( *this ).mv( x, y );
00099       }
00100 
00101       template< class X, class Y >
00102       void mtv ( const X &x, Y &y ) const
00103       {
00104         static_cast< const FieldMatrix & >( *this ).mtv( x, y );
00105       }
00106 
00107       template< class X, class Y >
00108       void umv ( const X &x, Y &y ) const
00109       {
00110         static_cast< const FieldMatrix & >( *this ).umv( x, y );
00111       }
00112 
00113       template< class X, class Y >
00114       void umtv ( const X &x, Y &y ) const
00115       {
00116         static_cast< const FieldMatrix & >( *this ).umtv( x, y );
00117       }
00118 
00119       template< class X, class Y >
00120       void mmv ( const X &x, Y &y ) const
00121       {
00122         static_cast< const FieldMatrix & >( *this ).mmv( x, y );
00123       }
00124 
00125       template< class X, class Y >
00126       void mmtv ( const X &x, Y &y ) const
00127       {
00128         static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
00129       }
00130 
00131       ctype det () const
00132       {
00133         if( !storage().integrationElementComputed )
00134         {
00135           storage().integrationElement = MatrixHelper::template sqrtDetAAT< rows, cols >( storage().jacobianTransposed );
00136           storage().integrationElementComputed = storage().affine;
00137         }
00138         return storage().integrationElement;
00139       }
00140 
00141     private:
00142       Storage &storage () const { return storage_; }
00143 
00144       mutable Storage storage_;
00145     };
00146 
00147 
00148 
00149     // CachedJacobianInverseTransposed
00150     // -------------------------------
00151 
00152     template< unsigned int dim, class GeometryTraits >
00153     class CachedJacobianInverseTransposed
00154     {
00155       template< class, class > friend class CachedMapping;
00156 
00157       typedef CachedJacobianTransposed< dim, GeometryTraits > JacobianTransposed;
00158       typedef typename JacobianTransposed::Storage Storage;
00159       typedef typename JacobianTransposed::Traits Traits;
00160 
00161       typedef typename Traits::MatrixHelper MatrixHelper;
00162 
00163     public:
00164       typedef typename Traits::FieldType ctype;
00165 
00166       static const int rows = Traits::dimWorld;
00167       static const int cols = Traits::dimension;
00168 
00169       typedef typename Traits::JacobianType FieldMatrix;
00170 
00171       operator bool () const
00172       {
00173         return storage().jacobianInverseTransposedComputed;
00174       }
00175 
00176       operator const FieldMatrix & () const
00177       {
00178         return storage().jacobianInverseTransposed;
00179       }
00180 
00181       template< class X, class Y >
00182       void mv ( const X &x, Y &y ) const
00183       {
00184         static_cast< const FieldMatrix & >( *this ).mv( x, y );
00185       }
00186 
00187       template< class X, class Y >
00188       void mtv ( const X &x, Y &y ) const
00189       {
00190         static_cast< const FieldMatrix & >( *this ).mtv( x, y );
00191       }
00192 
00193       template< class X, class Y >
00194       void umv ( const X &x, Y &y ) const
00195       {
00196         static_cast< const FieldMatrix & >( *this ).umv( x, y );
00197       }
00198 
00199       template< class X, class Y >
00200       void umtv ( const X &x, Y &y ) const
00201       {
00202         static_cast< const FieldMatrix & >( *this ).umtv( x, y );
00203       }
00204 
00205       template< class X, class Y >
00206       void mmv ( const X &x, Y &y ) const
00207       {
00208         static_cast< const FieldMatrix & >( *this ).mmv( x, y );
00209       }
00210 
00211       template< class X, class Y >
00212       void mmtv ( const X &x, Y &y ) const
00213       {
00214         static_cast< const FieldMatrix & >( *this ).mmtv( x, y );
00215       }
00216 
00217       ctype det () const
00218       {
00219         // integrationElement is always computed with jacobianInverseTransposed
00220         return ctype( 1 ) / storage().integrationElement;
00221       }
00222 
00223     private:
00224       JacobianTransposed &jacobianTransposed () { return jacobianTransposed_; }
00225       const JacobianTransposed &jacobianTransposed () const { return jacobianTransposed_; }
00226 
00227       Storage &storage () const { return jacobianTransposed().storage(); }
00228 
00229       JacobianTransposed jacobianTransposed_;
00230     };
00231 
00232 
00233 
00234     // CachedMapping
00235     // -------------
00236 
00237     template< class Topology, class GeometryTraits >
00238     class CachedMapping
00239     {
00240       typedef CachedMapping< Topology, GeometryTraits > This;
00241 
00242       typedef typename GeometryTraits::template Mapping< Topology >::type MappingImpl;
00243 
00244     public:
00245       typedef MappingTraits< typename GeometryTraits::CoordTraits, Topology::dimension, GeometryTraits::dimWorld > Traits;
00246 
00247       typedef GenericGeometry::Mapping< typename GeometryTraits::CoordTraits, Topology, GeometryTraits::dimWorld, MappingImpl > Mapping;
00248 
00249       static const unsigned int dimension = Traits::dimension;
00250       static const unsigned int dimWorld = Traits::dimWorld;
00251 
00252       typedef typename Traits::FieldType FieldType;           
00253       typedef typename Traits::LocalCoordinate LocalCoordinate;
00254       typedef typename Traits::GlobalCoordinate GlobalCoordinate;
00255 
00256       typedef CachedStorage< dimension, GeometryTraits > Storage;
00257       typedef CachedJacobianTransposed< dimension, GeometryTraits > JacobianTransposed;
00258       typedef CachedJacobianInverseTransposed< dimension, GeometryTraits > JacobianInverseTransposed;
00259 
00260       typedef GenericGeometry::ReferenceElement< Topology, FieldType > ReferenceElement;
00261 
00262       // can we safely assume that this mapping is always affine?
00263       static const bool alwaysAffine = Mapping::alwaysAffine;
00264 
00265 #if 0
00266       template< unsigned int codim >
00267       struct Codim
00268       {
00269         typedef typename TraceProvider< Topology, GeometryTraits, codim, false >::Trace Trace;
00270       };
00271 #endif
00272 
00273       typedef typename GeometryTraits::Caching Caching;
00274   
00275     private:
00276       typedef typename Traits::MatrixHelper MatrixHelper;
00277 
00278     public:
00279       template< class CoordVector >
00280       explicit CachedMapping ( const CoordVector &coords )
00281       : mapping_( coords )
00282       {
00283         if( alwaysAffine )
00284           storage().affine = true;
00285         else
00286           computeJacobianTransposed( baryCenter() );
00287         preCompute();
00288       }
00289 
00290       template< class CoordVector >
00291       explicit CachedMapping ( const std::pair< const CoordVector &, bool > &coords )
00292       : mapping_( coords.first )
00293       {
00294         storage().affine = coords.second;
00295         preCompute();
00296       }
00297 
00298       bool affine () const { return (alwaysAffine || storage().affine); }
00299       Dune::GeometryType type () const { return Dune::GeometryType( Topology() ); }
00300 
00301 #if 0
00302       unsigned int topologyId () const DUNE_DEPRECATED { return type().id(); }
00303 #endif
00304 
00305       int numCorners () const { return ReferenceElement::numCorners; }
00306       GlobalCoordinate corner ( int i ) const { return mapping().corner( i ); }
00307       GlobalCoordinate center () const { return global( ReferenceElement::template baryCenter< 0 >( 0 ) ); }
00308 
00309       static bool checkInside ( const LocalCoordinate &x ) { return ReferenceElement::checkInside( x ); }
00310 
00311       GlobalCoordinate global ( const LocalCoordinate &x ) const
00312       {
00313         GlobalCoordinate y;
00314         if( jacobianTransposed() )
00315         {
00316           y = corner( 0 );
00317           jacobianTransposed().umtv( x, y );
00318           //MatrixHelper::template ATx< dimension, dimWorld >( jacobianTransposed_, x, y );
00319           //y += corner( 0 );
00320         }
00321         else
00322           mapping().global( x, y );
00323         return y;
00324       }
00325 
00326       LocalCoordinate local ( const GlobalCoordinate &y ) const
00327       {
00328         LocalCoordinate x;
00329         if( jacobianInverseTransposed() )
00330         {
00331           GlobalCoordinate z = y - corner( 0 );
00332           jacobianInverseTransposed().mtv( z, x );
00333           // MatrixHelper::template ATx< dimWorld, dimension >( jacobianInverseTransposed(), z, x );
00334         }
00335         else if( affine() )
00336         {
00337           const JacobianTransposed &JT = jacobianTransposed( baryCenter() );
00338           GlobalCoordinate z = y - corner( 0 );
00339           MatrixHelper::template xTRightInvA< dimension, dimWorld >( JT, z, x );
00340         }
00341         else
00342           mapping().local( y, x );
00343         return x;
00344       }
00345 
00346       FieldType integrationElement ( const LocalCoordinate &x ) const
00347       {
00348         const EvaluationType evaluateI = Caching::evaluateIntegrationElement;
00349         const EvaluationType evaluateJ = Caching::evaluateJacobianInverseTransposed;
00350         if( ((evaluateI == PreCompute) || (evaluateJ == PreCompute)) && alwaysAffine )
00351           return storage().integrationElement;
00352         else
00353           return jacobianTransposed( x ).det();
00354       }
00355 
00356       FieldType volume () const
00357       {
00358         // do we need a quadrature of higher order, here?
00359         const FieldType refVolume = ReferenceElement::volume();
00360         return refVolume * integrationElement( baryCenter() );
00361       }
00362 
00363       const JacobianTransposed &jacobianTransposed ( const LocalCoordinate &x ) const
00364       {
00365         const EvaluationType evaluate = Caching::evaluateJacobianTransposed;
00366         if( (evaluate == PreCompute) && alwaysAffine )
00367           return jacobianTransposed();
00368 
00369         if( !jacobianTransposed() )
00370           computeJacobianTransposed( x );
00371         return jacobianTransposed();
00372       }
00373 
00374       const JacobianInverseTransposed &
00375       jacobianInverseTransposed ( const LocalCoordinate &x ) const
00376       {
00377         const EvaluationType evaluate = Caching::evaluateJacobianInverseTransposed;
00378         if( (evaluate == PreCompute) && alwaysAffine )
00379           return jacobianInverseTransposed();
00380 
00381         if( !jacobianInverseTransposed() )
00382           computeJacobianInverseTransposed( x );
00383         return jacobianInverseTransposed();
00384       }
00385 
00386 #if 0
00387       This *clone () const
00388       {
00389         return new This( *this );
00390       }
00391 
00392       This* clone ( char *mappingStorage ) const
00393       {
00394         return new( mappingStorage ) This( *this );
00395       }
00396 
00397       template< unsigned int codim, bool hybrid >
00398       typename TraceProvider< Topology, GeometryTraits, codim, hybrid >::Trace*
00399       trace ( unsigned int i, char *mappingStorage ) const
00400       {
00401         return TraceProvider< Topology, GeometryTraits, codim, hybrid >::construct( mapping(), i, mappingStorage );
00402       }
00403 #endif
00404 
00405       const Mapping &mapping () const { return mapping_; }
00406 
00407     private:
00408       static const LocalCoordinate &baryCenter ()
00409       {
00410         return ReferenceElement::template baryCenter< 0 >( 0 );
00411       }
00412 
00413       Storage &storage () const
00414       {
00415         return jacobianInverseTransposed().storage();
00416       }
00417 
00418       const JacobianTransposed &jacobianTransposed () const
00419       {
00420         return jacobianInverseTransposed().jacobianTransposed();
00421       }
00422 
00423       const JacobianInverseTransposed &jacobianInverseTransposed () const
00424       {
00425         return jacobianInverseTransposed_;
00426       }
00427 
00428       void preCompute ()
00429       {
00430         assert( affine() == mapping().jacobianTransposed( baryCenter(), storage().jacobianTransposed ) );
00431         if( !affine() )
00432           return;
00433 
00434         if( (Caching::evaluateJacobianTransposed == PreCompute) && !jacobianTransposed() )
00435           computeJacobianTransposed( baryCenter() );
00436 
00437         if( Caching::evaluateJacobianInverseTransposed == PreCompute )
00438           computeJacobianInverseTransposed( baryCenter() );
00439         else if( Caching::evaluateIntegrationElement == PreCompute )
00440           jacobianTransposed().det();
00441       }
00442 
00443       void computeJacobianTransposed ( const LocalCoordinate &x ) const
00444       {
00445         storage().affine = mapping().jacobianTransposed( x, storage().jacobianTransposed );
00446         storage().jacobianTransposedComputed = affine();
00447       }
00448 
00449       void computeJacobianInverseTransposed ( const LocalCoordinate &x ) const
00450       {
00451         storage().integrationElement
00452           = MatrixHelper::template rightInvA< dimension, dimWorld >( jacobianTransposed( x ), storage().jacobianInverseTransposed );
00453         storage().integrationElementComputed = affine();
00454         storage().jacobianInverseTransposedComputed = affine();
00455       }
00456 
00457     private:
00458       Mapping mapping_;
00459       JacobianInverseTransposed jacobianInverseTransposed_;
00460     };
00461 
00462   } // namespace GenericGeometry
00463 
00464 } // namespace Dune
00465 
00466 #endif // #ifndef DUNE_GEOMETRY_GENERICGEOMETRY_CACHED_MAPPING_HH