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