dune-geometry  2.2.0
referenceelements.hh
Go to the documentation of this file.
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