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