dune-geometry
2.2.0
|
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set et ts=8 sw=2 sts=2: 00003 #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH 00004 #define DUNE_GEOMETRY_REFERENCEELEMENTS_HH 00005 00006 #include <dune/common/forloop.hh> 00007 #include <dune/common/typetraits.hh> 00008 00009 #include <dune/geometry/genericgeometry/subtopologies.hh> 00010 #include <dune/geometry/genericgeometry/referencedomain.hh> 00011 #include <dune/geometry/genericgeometry/conversion.hh> 00012 #include <dune/geometry/genericgeometry/hybridmapping.hh> 00013 #include <dune/geometry/genericgeometry/mappingprovider.hh> 00014 00015 namespace Dune 00016 { 00017 00018 // Internal Forward Declarations 00019 // ----------------------------- 00020 00021 template< class ctype, int dim > 00022 class GenericReferenceElementContainer; 00023 00024 00025 00026 // GenericReferenceElement 00027 // ----------------------- 00028 00045 template< class ctype, int dim > 00046 class GenericReferenceElement 00047 { 00048 typedef GenericReferenceElement< ctype, dim > This; 00049 00050 friend class GenericReferenceElementContainer< ctype, dim >; 00051 00052 // make copy constructor private 00053 GenericReferenceElement(const GenericReferenceElement &); 00054 00055 GenericReferenceElement () {} 00056 00057 ~GenericReferenceElement () 00058 { 00059 ForLoop< Destroy, 0, dim >::apply( mappings_ ); 00060 integral_constant< int, 0 > codim0Variable; 00061 if(mappings_[ codim0Variable ].size()) 00062 delete mappings_[ codim0Variable ][ 0 ]; 00063 } 00064 00065 class SubEntityInfo; 00066 template< class Topology > class CornerStorage; 00067 template< class Topology > struct Initialize; 00068 template< int codim > struct Destroy; 00069 00070 struct GeometryTraits 00071 : public GenericGeometry::DefaultGeometryTraits< ctype, dim, dim > 00072 { 00073 typedef GenericGeometry::DefaultGeometryTraits< ctype, dim, dim > Base; 00074 00075 typedef typename Base::CoordTraits CoordTraits; 00076 00077 template< class Topology > 00078 struct Mapping 00079 { 00080 typedef GenericGeometry::CornerMapping< CoordTraits, Topology, dim, CornerStorage< Topology >, true > type; 00081 }; 00082 00083 struct Caching 00084 { 00085 static const GenericGeometry::EvaluationType evaluateJacobianTransposed = GenericGeometry::PreCompute; 00086 static const GenericGeometry::EvaluationType evaluateJacobianInverseTransposed = GenericGeometry::PreCompute; 00087 static const GenericGeometry::EvaluationType evaluateIntegrationElement = GenericGeometry::PreCompute; 00088 static const GenericGeometry::EvaluationType evaluateNormal = GenericGeometry::PreCompute; 00089 }; 00090 00091 }; 00092 00093 public: 00095 template< int codim > 00096 struct Codim 00097 { 00099 typedef GenericGeometry::HybridMapping< dim-codim, GeometryTraits > Mapping; 00100 }; 00101 00102 private: 00104 template< int codim > 00105 struct MappingArray 00106 : public std::vector< typename Codim< codim >::Mapping * > 00107 {}; 00108 00110 typedef GenericGeometry::CodimTable< MappingArray, dim > MappingsTable; 00111 00112 std::vector< SubEntityInfo > info_[ dim+1 ]; 00113 00115 ctype volume_; 00116 std::vector< FieldVector< ctype, dim > > volumeNormals_; 00117 00119 MappingsTable mappings_; 00120 00121 public: 00126 int size ( int c ) const 00127 { 00128 assert( (c >= 0) && (c <= dim) ); 00129 return info_[ c ].size(); 00130 } 00131 00143 int size ( int i, int c, int cc ) const 00144 { 00145 assert( (c >= 0) && (c <= dim) ); 00146 return info_[ c ][ i ].size( cc ); 00147 } 00148 00162 int subEntity ( int i, int c, int ii, int cc ) const 00163 { 00164 assert( (c >= 0) && (c <= dim) ); 00165 return info_[ c ][ i ].number( ii, cc ); 00166 } 00167 00177 const FieldVector< ctype, dim > &position( int i, int c ) const 00178 { 00179 assert( (c >= 0) && (c <= dim) ); 00180 return info_[ c ][ i ].position(); 00181 } 00182 00190 bool checkInside ( const FieldVector< ctype, dim > &local ) const 00191 { 00192 return checkInside< 0 >( local, 0 ); 00193 } 00194 00209 template< int codim > 00210 bool checkInside ( const FieldVector< ctype, dim-codim > &local, int i ) const 00211 { 00212 return mapping< codim >( i ).checkInside( local ); 00213 } 00214 00236 template< int codim > 00237 FieldVector< ctype, dim > 00238 global( const FieldVector< ctype, dim-codim > &local, int i, int c ) const 00239 { 00240 if( c != codim ) 00241 DUNE_THROW( Exception, "Local Coordinate Type does not correspond to codimension c." ); 00242 assert( c == codim ); 00243 return mapping< codim >( i ).global( local ); 00244 } 00245 00264 template< int codim > 00265 FieldVector< ctype, dim > 00266 global( const FieldVector< ctype, dim-codim > &local, int i ) const 00267 { 00268 return mapping< codim >( i ).global( local ); 00269 } 00270 00286 template< int codim > 00287 typename Codim< codim >::Mapping &mapping( int i ) const 00288 { 00289 integral_constant< int, codim > codimVariable; 00290 return *(mappings_[ codimVariable ][ i ]); 00291 } 00292 00301 const GeometryType &type ( int i, int c ) const 00302 { 00303 assert( (c >= 0) && (c <= dim) ); 00304 return info_[ c ][ i ].type(); 00305 } 00306 00308 const GeometryType &type () const { return type( 0, 0 ); } 00309 00310 unsigned int topologyId ( int i, int c ) const DUNE_DEPRECATED 00311 { 00312 assert( (c >= 0) && (c <= dim) ); 00313 return info_[ c ][ i ].topologyId(); 00314 } 00315 00317 ctype volume () const 00318 { 00319 return volume_; 00320 } 00321 00329 const FieldVector< ctype, dim > &volumeOuterNormal ( int face ) const 00330 { 00331 assert( (face >= 0) && (face < int( volumeNormals_.size())) ); 00332 return volumeNormals_[ face ]; 00333 } 00334 00341 template< class Topology > 00342 void initializeTopology () 00343 { 00344 dune_static_assert( (Topology::dimension == dim), 00345 "Cannot initialize reference element for different dimension." ); 00346 typedef Initialize< Topology > Init; 00347 typedef GenericGeometry::VirtualMapping< Topology, GeometryTraits > VirtualMapping; 00348 00349 // set up subentities 00350 integral_constant< int, 0 > codim0Variable; 00351 mappings_[ codim0Variable ].resize( 1 ); 00352 mappings_[ codim0Variable ][ 0 ] = new VirtualMapping( codim0Variable ); 00353 00354 Dune::ForLoop< Init::template Codim, 0, dim >::apply( info_, mappings_ ); 00355 00356 // compute reference element volume 00357 typedef GenericGeometry::ReferenceDomain< Topology > ReferenceDomain; 00358 volume_ = ReferenceDomain::template volume< ctype >(); 00359 00360 // compute normals 00361 volumeNormals_.resize( ReferenceDomain::numNormals ); 00362 for( unsigned int i = 0; i < ReferenceDomain::numNormals; ++i ) 00363 ReferenceDomain::integrationOuterNormal( i ,volumeNormals_[ i ] ); 00364 } 00365 }; 00366 00367 00371 template< class ctype, int dim > 00372 class GenericReferenceElement< ctype, dim >::SubEntityInfo 00373 { 00374 template< class Topology, int codim > struct Initialize 00375 { 00376 template< int subcodim > struct SubCodim; 00377 }; 00378 00379 int codim_; 00380 std::vector< int > numbering_[ dim+1 ]; 00381 FieldVector< ctype, dim > baryCenter_; 00382 GeometryType type_; 00383 00384 public: 00385 int size ( int cc ) const 00386 { 00387 assert( (cc >= codim_) && (cc <= dim) ); 00388 return numbering_[ cc ].size(); 00389 } 00390 00391 int number ( int ii, int cc ) const 00392 { 00393 assert( (cc >= codim_) && (cc <= dim) ); 00394 return numbering_[ cc ][ ii ]; 00395 } 00396 00397 const FieldVector< ctype, dim > &position () const 00398 { 00399 return baryCenter_; 00400 } 00401 00402 const GeometryType &type () const 00403 { 00404 return type_; 00405 } 00406 00407 unsigned int topologyId () const DUNE_DEPRECATED 00408 { 00409 return type_.id(); 00410 } 00411 00412 template< class Topology, unsigned int codim, unsigned int i > 00413 void initialize () 00414 { 00415 typedef Initialize< Topology, codim > Init; 00416 typedef GenericGeometry::ReferenceDomain< Topology > RefDomain; 00417 00418 codim_ = codim; 00419 00420 const unsigned int iVariable = i; 00421 Dune::ForLoop< Init::template SubCodim, 0, dim-codim >::apply( iVariable, numbering_ ); 00422 00423 baryCenter_ = ctype( 0 ); 00424 static const unsigned int numCorners = size( dim ); 00425 for( unsigned int j = 0; j < numCorners; ++j ) 00426 { 00427 FieldVector< ctype, dim > corner; 00428 RefDomain::corner( number( j, dim ), corner ); 00429 baryCenter_ += corner; 00430 } 00431 baryCenter_ *= ctype( 1 ) / ctype( numCorners ); 00432 00433 typedef typename GenericGeometry::SubTopology< Topology, codim, i >::type SubTopology; 00434 type_ = GeometryType( SubTopology::id, SubTopology::dimension ); 00435 // type_ = GenericGeometry::DuneGeometryType< SubTopology, GeometryType::simplex >::type(); 00436 } 00437 }; 00438 00439 00440 template< class ctype, int dim > 00441 template< class Topology > 00442 class GenericReferenceElement< ctype, dim >::CornerStorage 00443 { 00444 typedef GenericGeometry::ReferenceDomain< Topology > RefDomain; 00445 00446 public: 00447 static const unsigned int size = Topology::numCorners; 00448 00449 template< class SubTopology > 00450 struct SubStorage 00451 { 00452 typedef CornerStorage< SubTopology > type; 00453 }; 00454 00455 explicit CornerStorage ( const integral_constant< int, 0 > & ) 00456 { 00457 for( unsigned int i = 0; i < size; ++i ) 00458 RefDomain::corner( i, coords_[ i ] ); 00459 } 00460 00461 template< class Mapping, unsigned int codim > 00462 explicit 00463 CornerStorage ( const GenericGeometry::SubMappingCoords< Mapping, codim > &coords ) 00464 { 00465 for( unsigned int i = 0; i < size; ++i ) 00466 coords_[ i ] = coords[ i ]; 00467 } 00468 00469 const FieldVector< ctype, dim > &operator[] ( unsigned int i ) const 00470 { 00471 return coords_[ i ]; 00472 } 00473 00474 private: 00475 FieldVector< ctype, dim > coords_[ size ]; 00476 }; 00477 00478 00479 template< class ctype, int dim > 00480 template< class Topology, int codim > 00481 template< int subcodim > 00482 struct GenericReferenceElement< ctype, dim >::SubEntityInfo::Initialize< Topology, codim >::SubCodim 00483 { 00484 typedef GenericGeometry::SubTopologySize< Topology, codim, subcodim > SubSize; 00485 typedef GenericGeometry::GenericSubTopologyNumbering< Topology, codim, subcodim > SubNumbering; 00486 00487 static void apply ( unsigned int i, std::vector< int > (&numbering)[ dim+1 ] ) 00488 { 00489 const unsigned int size = SubSize::size( i ); 00490 numbering[ codim+subcodim ].resize( size ); 00491 for( unsigned int j = 0; j < size; ++j ) 00492 numbering[ codim+subcodim ][ j ] = SubNumbering::number( i, j ); 00493 } 00494 }; 00495 00496 00497 template< class ctype, int dim > 00498 template< class Topology > 00499 struct GenericReferenceElement< ctype, dim >::Initialize 00500 { 00501 typedef Dune::GenericReferenceElement< ctype, dim > GenericReferenceElement; 00502 00503 typedef typename GenericReferenceElement::template Codim< 0 >::Mapping ReferenceMapping; 00504 00505 template< int codim > 00506 struct Codim 00507 { 00508 template< int i > 00509 struct SubTopology 00510 { 00511 static void apply ( std::vector< SubEntityInfo > &info ) 00512 { 00513 info[ i ].template initialize< Topology, codim, i >(); 00514 } 00515 }; 00516 00517 static void 00518 apply ( std::vector< SubEntityInfo > (&info)[ dim+1 ], 00519 MappingsTable &mappings ) 00520 { 00521 const unsigned int size = GenericGeometry::Size< Topology, codim >::value; 00522 info[ codim ].resize( size ); 00523 Dune::ForLoop< SubTopology, 0, size-1 >::apply( info[ codim ] ); 00524 00525 if( codim > 0 ) 00526 { 00527 integral_constant< int, 0 > codim0Variable; 00528 const ReferenceMapping &refMapping = *(mappings[ codim0Variable ][ 0 ]); 00529 00530 typedef typename GenericGeometry::MappingProvider< ReferenceMapping, codim > MappingProvider; 00531 00532 integral_constant< int, codim > codimVariable; 00533 mappings[ codimVariable ].resize( size ); 00534 for( unsigned int i = 0; i < size; ++i ) { 00535 char* storage = new char[MappingProvider::maxMappingSize]; 00536 mappings[ codimVariable ][ i ] = refMapping.template trace< codim >( i, storage ); 00537 } 00538 } 00539 } 00540 }; 00541 }; 00542 00543 00544 00545 template< class ctype, int dim > 00546 template< int codim > 00547 struct GenericReferenceElement< ctype, dim >::Destroy 00548 { 00549 static void apply ( MappingsTable &mappings ) 00550 { 00551 if (codim > 0 ) 00552 { 00553 integral_constant< int, codim > codimVariable; 00554 for( size_t i = 0; i < mappings[ codimVariable ].size(); ++i ) { 00555 typedef typename Codim<codim>::Mapping Mapping; 00556 mappings[ codimVariable ][ i ]->~Mapping(); 00557 char* storage = (char*)mappings[ codimVariable ][ i ]; 00558 delete[](storage); 00559 } 00560 } 00561 } 00562 }; 00563 00564 00565 00566 // GenericReferenceElementContainer 00567 // -------------------------------- 00568 00569 template< class ctype, int dim > 00570 class GenericReferenceElementContainer 00571 { 00572 static const unsigned int numTopologies = (1u << dim); 00573 00574 public: 00575 typedef GenericReferenceElement< ctype, dim > value_type; 00576 typedef const value_type *const_iterator; 00577 00578 GenericReferenceElementContainer () 00579 { 00580 ForLoop< Builder, 0, numTopologies-1 >::apply( values_ ); 00581 } 00582 00583 const value_type &operator() ( const unsigned int topologyId ) const DUNE_DEPRECATED 00584 { 00585 return values_[ topologyId ]; 00586 } 00587 00588 const value_type &operator() ( const GeometryType &type ) const 00589 { 00590 assert( type.dim() == dim ); 00591 return values_[ type.id() ]; 00592 } 00593 00594 const value_type &simplex () const 00595 { 00596 return values_[ GenericGeometry::SimplexTopology< dim >::type::id ]; 00597 } 00598 00599 const value_type &cube () const 00600 { 00601 return values_[ GenericGeometry::CubeTopology< dim >::type::id ]; 00602 } 00603 00604 const value_type &pyramid () const 00605 { 00606 return values_[ GenericGeometry::PyramidTopology< dim >::type::id ]; 00607 } 00608 00609 const value_type &prism () const 00610 { 00611 return values_[ GenericGeometry::PrismTopology< dim >::type::id ]; 00612 } 00613 00614 const_iterator begin () const { return values_; } 00615 const_iterator end () const { return values_ + numTopologies; } 00616 00617 static const GenericReferenceElementContainer &instance () DUNE_DEPRECATED 00618 { 00619 static GenericReferenceElementContainer inst; 00620 return inst; 00621 } 00622 00623 private: 00624 template< int topologyId > 00625 struct Builder 00626 { 00627 static void apply ( value_type (&values)[ numTopologies ] ) 00628 { 00629 typedef typename GenericGeometry::Topology< topologyId, dim >::type Topology; 00630 values[ topologyId ].template initializeTopology< Topology >(); 00631 } 00632 }; 00633 00634 value_type values_[ numTopologies ]; 00635 }; 00636 00637 00638 00639 // GenericReferenceElements 00640 // ------------------------ 00641 00650 template< class ctype, int dim > 00651 struct GenericReferenceElements 00652 { 00653 typedef typename GenericReferenceElementContainer< ctype, dim >::const_iterator Iterator; 00654 00656 static const GenericReferenceElement< ctype, dim > & 00657 general ( const GeometryType &type ) 00658 { 00659 return container()( type ); 00660 } 00661 00663 static const GenericReferenceElement< ctype, dim > &simplex () 00664 { 00665 return container().simplex(); 00666 } 00667 00669 static const GenericReferenceElement< ctype, dim > &cube () 00670 { 00671 return container().cube(); 00672 } 00673 00674 static Iterator begin () { return container().begin(); } 00675 static Iterator end () { return container().end(); } 00676 00677 private: 00678 static const GenericReferenceElementContainer< ctype, dim > &container () 00679 { 00680 static GenericReferenceElementContainer< ctype, dim > container; 00681 return container; 00682 } 00683 }; 00684 00685 } // namespace Dune 00686 00687 #endif // #ifndef DUNE_GEOMETRY_REFERENCEELEMENTS_HH