dune-localfunctions  2.2.0
monomialbasis.hh
Go to the documentation of this file.
00001 #ifndef DUNE_MONOMIALBASIS_HH
00002 #define DUNE_MONOMIALBASIS_HH
00003 
00004 #include <vector>
00005 
00006 #include <dune/common/fvector.hh>
00007 #include <dune/common/fmatrix.hh>
00008 
00009 #include <dune/geometry/topologyfactory.hh>
00010 #include <dune/geometry/genericgeometry/topologytypes.hh>
00011 
00012 #include <dune/localfunctions/utility/field.hh>
00013 #include <dune/localfunctions/utility/multiindex.hh>
00014 #include <dune/localfunctions/utility/tensor.hh>
00015 
00016 namespace Dune
00017 {
00018 /************************************************
00019  * Classes for evaluating ''Monomials'' on any order
00020  * for all reference element type. 
00021  * For a simplex topology these are the nomral 
00022  * monomials for cube topologies the bimonomials.
00023  * The construction follows the construction of the
00024  * generic geometries using tensor products for
00025  * prism generation and duffy transform for pyramid
00026  * construction.
00027  * A derivative argument can be applied, in which case
00028  * all derivatives up to the desired order are 
00029  * evaluated. Note that in for higher order derivatives
00030  * only the ''lower'' part of the symmetric tensor
00031  * is evaluated, e.g., passing derivative equal to 2
00032  * to the class will provide the vector
00033  *    (d/dxdx p, d/dxydx p, d/dydy p, 
00034  *     d/dx p, d/dy p, p)
00035  * Important:
00036  * So far the computation of the derivatives has not
00037  * been fully implemented for general pyramid
00038  * construction, i.e., in the case where a pyramid is
00039  * build over a non simplex base geometry.
00040  *
00041  * Central classes:
00042  * 1) template< class Topology, class F >
00043  *    class MonomialBasisImpl;
00044  *    Implementation of the monomial evaluation for
00045  *    a given topology and field type.
00046  *    The method evaluate fills a F* vector
00047  * 2) template< class Topology, class F >
00048  *    class MonomialBasis
00049  *    The base class for the static monomial evaluation
00050  *    providing addiional evaluate methods including
00051  *    one taking std::vector<F>.
00052  * 3) template< int dim, class F >
00053  *    class VirtualMonomialBasis
00054  *    Virtualization of the MonomialBasis.
00055  * 4) template< int dim, class F >
00056  *    struct MonomialBasisFactory;
00057  *    A factory class for the VirtualMonomialBasis
00058  * 5) template< int dim, class F >
00059  *    struct MonomialBasisProvider
00060  *    A singleton container for the virtual monomial
00061  *    basis
00062  ************************************************/
00063 
00064   // Internal Forward Declarations
00065   // -----------------------------
00066   
00067   template< class Topology >
00068   class MonomialBasisSize;
00069 
00070   template< class Topology, class F >
00071   class MonomialBasis;
00072 
00073 
00074 
00075   // MonomialBasisSize
00076   // -----------------
00077 
00078   template<>
00079   class MonomialBasisSize< GenericGeometry::Point >
00080   {
00081     typedef MonomialBasisSize< GenericGeometry::Point > This;
00082 
00083   public:
00084     static This &instance ()
00085     {
00086       static This _instance;
00087       return _instance;
00088     }
00089 
00090     typedef GenericGeometry::Point Topology;
00091 
00092     friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
00093     friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
00094 
00095     mutable unsigned int maxOrder_;
00096     // sizes_[ k ]: number of basis functions of exactly order k
00097     mutable unsigned int *sizes_;
00098     // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
00099     mutable unsigned int *numBaseFunctions_;
00100 
00101     MonomialBasisSize ()
00102     : maxOrder_( 0 ),
00103       sizes_( 0 ),
00104       numBaseFunctions_( 0 )
00105     {
00106       computeSizes( 2 );
00107     }
00108 
00109     ~MonomialBasisSize ()
00110     {
00111       delete[] sizes_;
00112       delete[] numBaseFunctions_;
00113     }
00114 
00115     unsigned int operator() ( const unsigned int order ) const
00116     {
00117       return numBaseFunctions_[ order ];
00118     }
00119 
00120     unsigned int maxOrder () const
00121     {
00122       return maxOrder_;
00123     }
00124 
00125     void computeSizes ( unsigned int order ) const
00126     {
00127       if (order <= maxOrder_)
00128         return;
00129 
00130       maxOrder_ = order;
00131 
00132       delete [] sizes_;
00133       delete [] numBaseFunctions_;
00134       sizes_            = new unsigned int [ order+1 ];
00135       numBaseFunctions_ = new unsigned int [ order+1 ];
00136 
00137       sizes_[ 0 ] = 1;
00138       numBaseFunctions_[ 0 ] = 1;
00139       for( unsigned int k = 1; k <= order; ++k )
00140       {
00141         sizes_[ k ]            = 0;
00142         numBaseFunctions_[ k ] = 1;
00143       }
00144     }
00145   };
00146 
00147   template< class BaseTopology >
00148   class MonomialBasisSize< GenericGeometry::Prism< BaseTopology > >
00149   {
00150     typedef MonomialBasisSize< GenericGeometry::Prism< BaseTopology > > This;
00151 
00152   public:
00153     static This &instance ()
00154     {
00155       static This _instance;
00156       return _instance;
00157     }
00158 
00159     typedef GenericGeometry::Prism< BaseTopology > Topology;
00160 
00161     friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
00162     friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
00163 
00164     mutable unsigned int maxOrder_;
00165     // sizes_[ k ]: number of basis functions of exactly order k
00166     mutable unsigned int *sizes_;
00167     // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
00168     mutable unsigned int *numBaseFunctions_;
00169 
00170     MonomialBasisSize ()
00171     : maxOrder_( 0 ),
00172       sizes_( 0 ),
00173       numBaseFunctions_( 0 )
00174     {
00175       computeSizes( 2 );
00176     }
00177 
00178     ~MonomialBasisSize ()
00179     {
00180       delete[] sizes_;
00181       delete[] numBaseFunctions_;
00182     }
00183 
00184     unsigned int operator() ( const unsigned int order ) const
00185     {
00186       return numBaseFunctions_[ order ];
00187     }
00188 
00189     unsigned int maxOrder() const
00190     {
00191       return maxOrder_;
00192     }
00193 
00194     void computeSizes ( unsigned int order ) const
00195     {
00196       if (order <= maxOrder_)
00197         return;
00198 
00199       maxOrder_ = order;
00200 
00201       delete[] sizes_;
00202       delete[] numBaseFunctions_;
00203       sizes_            = new unsigned int[ order+1 ];
00204       numBaseFunctions_ = new unsigned int[ order+1 ];
00205 
00206       MonomialBasisSize<BaseTopology> &baseBasis =
00207             MonomialBasisSize<BaseTopology>::instance();
00208       baseBasis.computeSizes( order );
00209       const unsigned int *const baseSizes = baseBasis.sizes_;
00210       const unsigned int *const baseNBF   = baseBasis.numBaseFunctions_;
00211 
00212       sizes_[ 0 ] = 1;
00213       numBaseFunctions_[ 0 ] = 1;
00214       for( unsigned int k = 1; k <= order; ++k )
00215       {
00216         sizes_[ k ]            = baseNBF[ k ] + k*baseSizes[ k ];
00217         numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
00218       }
00219     }
00220   };
00221 
00222   template< class BaseTopology >
00223   class MonomialBasisSize< GenericGeometry::Pyramid< BaseTopology > >
00224   {
00225     typedef MonomialBasisSize< GenericGeometry::Pyramid< BaseTopology > > This;
00226 
00227   public:
00228     static This &instance ()
00229     {
00230       static This _instance;
00231       return _instance;
00232     }
00233 
00234     typedef GenericGeometry::Pyramid< BaseTopology > Topology;
00235 
00236     friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
00237     friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
00238 
00239     mutable unsigned int maxOrder_;
00240     // sizes_[ k ]: number of basis functions of exactly order k
00241     mutable unsigned int *sizes_;
00242     // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
00243     mutable unsigned int *numBaseFunctions_;
00244 
00245     MonomialBasisSize ()
00246     : maxOrder_( 0 ),
00247       sizes_( 0 ),
00248       numBaseFunctions_( 0 )
00249     {
00250       computeSizes( 2 );
00251     }
00252 
00253     ~MonomialBasisSize ()
00254     {
00255       delete[] sizes_;
00256       delete[] numBaseFunctions_;
00257     }
00258 
00259     unsigned int operator() ( const unsigned int order ) const
00260     {
00261       return numBaseFunctions_[ order ];
00262     }
00263 
00264     unsigned int maxOrder() const
00265     {
00266       return maxOrder_;
00267     }
00268 
00269     void computeSizes ( unsigned int order ) const
00270     {
00271       if (order <= maxOrder_)
00272         return;
00273 
00274       maxOrder_ = order;
00275 
00276       delete[] sizes_;
00277       delete[] numBaseFunctions_;
00278       sizes_            = new unsigned int[ order+1 ];
00279       numBaseFunctions_ = new unsigned int[ order+1 ];
00280 
00281       MonomialBasisSize<BaseTopology> &baseBasis =
00282             MonomialBasisSize<BaseTopology>::instance();
00283 
00284       baseBasis.computeSizes( order );
00285 
00286       const unsigned int *const baseNBF = baseBasis.numBaseFunctions_;
00287       sizes_[ 0 ] = 1;
00288       numBaseFunctions_[ 0 ] = 1;
00289       for( unsigned int k = 1; k <= order; ++k )
00290       {
00291         sizes_[ k ]            = baseNBF[ k ];
00292         numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
00293       }
00294     }
00295   };
00296 
00297 
00298 
00299   // MonomialBasisHelper
00300   // -------------------
00301 
00302 
00303   template< int mydim, int dim, class F >
00304   struct MonomialBasisHelper
00305   {
00306     typedef MonomialBasisSize< typename GenericGeometry::SimplexTopology< mydim >::type > MySize;
00307     typedef MonomialBasisSize< typename GenericGeometry::SimplexTopology< dim >::type > Size;
00308 
00309     static void copy ( const unsigned int deriv, F *&wit, F *&rit,
00310                        const unsigned int numBaseFunctions, const F &z )
00311     {
00312       // n(d,k) = size<k>[d];
00313       MySize &mySize = MySize::instance();
00314       Size &size = Size::instance();
00315 
00316       const F *const rend = rit + size( deriv )*numBaseFunctions;
00317       for( ; rit != rend; )
00318       {
00319         F *prit = rit;
00320 
00321         *wit = z * *rit; 
00322         ++rit, ++wit; 
00323 
00324         for( unsigned d = 1; d <= deriv; ++d )
00325         {
00326           #ifndef NDEBUG
00327           const F *const derivEnd = rit + mySize.sizes_[ d ];
00328           #endif
00329           const F *const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
00330           for( ; rit != drend ; ++rit, ++wit )
00331             *wit = z * *rit; 
00332           for (unsigned int j=1;j<d;++j) 
00333           {
00334             const F *const drend = rit + mySize.sizes_[ d-j ] - mySize.sizes_[ d-j-1 ];
00335             for( ; rit != drend ; ++prit, ++rit, ++wit )
00336               *wit = F(j) * *prit + z * *rit; 
00337           }
00338           *wit = F(d) * *prit + z * *rit; 
00339           ++prit, ++rit, ++wit;
00340           assert(derivEnd == rit);
00341           rit += size.sizes_[d] - mySize.sizes_[d];
00342           prit += size.sizes_[d-1] - mySize.sizes_[d-1];
00343           const F *const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
00344           for ( ; wit != emptyWitEnd; ++wit )
00345             *wit = Zero<F>();
00346         }
00347       }
00348     }
00349   };
00350 
00351 
00352 
00353   // MonomialBasisImpl
00354   // -----------------
00355 
00356   template< class Topology, class F >
00357   class MonomialBasisImpl;
00358 
00359   template< class F >
00360   class MonomialBasisImpl< GenericGeometry::Point, F >
00361   {
00362     typedef MonomialBasisImpl< GenericGeometry::Point, F > This;
00363 
00364   public:
00365     typedef GenericGeometry::Point Topology;
00366     typedef F Field;
00367 
00368     static const unsigned int dimDomain = Topology::dimension;
00369 
00370     typedef FieldVector< Field, dimDomain > DomainVector;
00371 
00372   private:
00373     friend class MonomialBasis< Topology, Field >;
00374     friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
00375     friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
00376 
00377     template< int dimD >
00378     void evaluate ( const unsigned int deriv, const unsigned int order,
00379                     const FieldVector< Field, dimD > &x,
00380                     const unsigned int block, const unsigned int *const offsets,
00381                     Field *const values ) const  
00382     {
00383       *values = Unity< F >();
00384       F *const end = values + block;
00385       for( Field *it = values+1 ; it != end; ++it )
00386         *it = Zero< F >();
00387     }
00388 
00389     void integrate ( const unsigned int order,
00390                      const unsigned int *const offsets,
00391                      Field *const values ) const
00392     {
00393       values[ 0 ] = Unity< Field >();
00394     }
00395   };
00396 
00397   template< class BaseTopology, class F >
00398   class MonomialBasisImpl< GenericGeometry::Prism< BaseTopology >, F >
00399   {
00400     typedef MonomialBasisImpl< GenericGeometry::Prism< BaseTopology >, F > This;
00401 
00402   public:
00403     typedef GenericGeometry::Prism< BaseTopology > Topology;
00404     typedef F Field;
00405 
00406     static const unsigned int dimDomain = Topology::dimension;
00407 
00408     typedef FieldVector< Field, dimDomain > DomainVector;
00409 
00410   private:
00411     friend class MonomialBasis< Topology, Field >;
00412     friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
00413     friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
00414 
00415     typedef MonomialBasisSize< BaseTopology > BaseSize; 
00416     typedef MonomialBasisSize< Topology > Size; 
00417 
00418     MonomialBasisImpl< BaseTopology, Field > baseBasis_;
00419 
00420     MonomialBasisImpl ()
00421     {}
00422 
00423     template< int dimD >
00424     void evaluate ( const unsigned int deriv, const unsigned int order,
00425                     const FieldVector< Field, dimD > &x,
00426                     const unsigned int block, const unsigned int *const offsets,
00427                     Field *const values ) const
00428     {
00429       typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
00430       const BaseSize &size = BaseSize::instance();
00431 
00432       const Field &z = x[ dimDomain-1 ];
00433       
00434       // fill first column
00435       baseBasis_.evaluate( deriv, order, x, block, offsets, values );
00436 
00437       Field *row0 = values;
00438       for( unsigned int k = 1; k <= order; ++k )
00439       {
00440         Field *row1 = values + block*offsets[ k-1 ];
00441         Field *wit = row1 + block*size.sizes_[ k ];
00442         Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
00443         Helper::copy( deriv, wit, row0, size( k-1 ), z );
00444         row0 = row1;
00445       }
00446     }
00447 
00448     void integrate ( const unsigned int order,
00449                      const unsigned int *const offsets,
00450                      Field *const values ) const
00451     {
00452       const BaseSize &size = BaseSize::instance();
00453       const Size &mySize = Size::instance();
00454       // fill first column
00455       baseBasis_.integrate( order, offsets, values );
00456       const unsigned int *const baseSizes = size.sizes_;
00457 
00458       Field *row0 = values;
00459       for( unsigned int k = 1; k <= order; ++k )
00460       {
00461         Field *const row1begin = values + offsets[ k-1 ];
00462         Field *const row1End = row1begin + mySize.sizes_[ k ];
00463         assert( (unsigned int)(row1End - values) <= offsets[ k ] );
00464 
00465         Field *row1 = row1begin;
00466         Field *it = row1begin + baseSizes[ k ];
00467         for( unsigned int j = 1; j <= k; ++j )
00468         {
00469           Field *const end = it + baseSizes[ k ];
00470           assert( (unsigned int)(end - values) <= offsets[ k ] );
00471           for( ; it != end; ++row1, ++it )
00472             *it = (Field( j ) / Field( j+1 )) * (*row1);
00473         }
00474         for( ; it != row1End; ++row0, ++it )
00475           *it = (Field( k ) / Field( k+1 )) * (*row0);
00476         row0 = row1;
00477       }
00478     }
00479   };
00480 
00481   template< class BaseTopology, class F >
00482   class MonomialBasisImpl< GenericGeometry::Pyramid< BaseTopology >, F >
00483   {
00484     typedef MonomialBasisImpl< GenericGeometry::Pyramid< BaseTopology >, F > This;
00485 
00486   public:
00487     typedef GenericGeometry::Pyramid< BaseTopology > Topology;
00488     typedef F Field;
00489 
00490     static const unsigned int dimDomain = Topology::dimension;
00491 
00492     typedef FieldVector< Field, dimDomain > DomainVector;
00493 
00494   private:
00495     friend class MonomialBasis< Topology, Field >;
00496     friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
00497     friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
00498 
00499     typedef MonomialBasisSize< BaseTopology > BaseSize; 
00500     typedef MonomialBasisSize< Topology > Size; 
00501 
00502     MonomialBasisImpl< BaseTopology, Field > baseBasis_;
00503 
00504     MonomialBasisImpl ()
00505     {
00506     }
00507 
00508     template< int dimD >
00509     void evaluateSimplexBase ( const unsigned int deriv, const unsigned int order,
00510                                const FieldVector< Field, dimD > &x,
00511                                const unsigned int block, const unsigned int *const offsets,
00512                                Field *const values, 
00513                                const BaseSize &size ) const
00514     {
00515       baseBasis_.evaluate( deriv, order, x, block, offsets, values );
00516     }
00517 
00518     template< int dimD >
00519     void evaluatePyramidBase ( const unsigned int deriv, const unsigned int order,
00520                                const FieldVector< Field, dimD > &x,
00521                                const unsigned int block, const unsigned int *const offsets,
00522                                Field *const values,
00523                                const BaseSize &size ) const
00524     {
00525       Field omz = Unity< Field >() - x[ dimDomain-1 ];
00526 
00527       if( Zero< Field >() < omz ) 
00528       {
00529         const Field invomz = Unity< Field >() / omz;
00530         FieldVector< Field, dimD > y;
00531         for( unsigned int i = 0; i < dimDomain-1; ++i )
00532           y[ i ] = x[ i ] * invomz;
00533       
00534         // fill first column
00535         baseBasis_.evaluate( deriv, order, y, block, offsets, values );
00536 
00537         Field omzk = omz;
00538         for( unsigned int k = 1; k <= order; ++k )
00539         {
00540           Field *it = values + block*offsets[ k-1 ];
00541           Field *const end = it + block*size.sizes_[ k ];
00542           for( ; it != end; ++it )
00543             *it *= omzk;
00544           omzk *= omz;
00545         }
00546       }
00547       else
00548       {
00549         assert( deriv==0 );
00550         *values = Unity< Field >();
00551         for( unsigned int k = 1; k <= order; ++k )
00552         {
00553           Field *it = values + block*offsets[ k-1 ];
00554           Field *const end = it + block*size.sizes_[ k ];
00555           for( ; it != end; ++it )
00556             *it = Zero< Field >();
00557         }
00558       }
00559     }
00560 
00561     template< int dimD >
00562     void evaluate ( const unsigned int deriv, const unsigned int order,
00563                     const FieldVector< Field, dimD > &x,
00564                     const unsigned int block, const unsigned int *const offsets,
00565                     Field *const values ) const
00566     {
00567       typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
00568       const BaseSize &size = BaseSize::instance();
00569       
00570       if( GenericGeometry::IsSimplex< Topology >::value )
00571         evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
00572       else
00573         evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
00574 
00575       Field *row0 = values;
00576       for( unsigned int k = 1; k <= order; ++k )
00577       {
00578         Field *row1 = values + block*offsets[ k-1 ];
00579         Field *wit = row1 + block*size.sizes_[ k ];
00580         Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
00581         row0 = row1;
00582       }
00583     }
00584 
00585     void integrate ( const unsigned int order,
00586                      const unsigned int *const offsets,
00587                      Field *const values ) const
00588     {
00589       const BaseSize &size = BaseSize::instance();
00590 
00591       // fill first column
00592       baseBasis_.integrate( order, offsets, values );
00593 
00594       const unsigned int *const baseSizes = size.sizes_;
00595 
00596       Field *const col0End = values + baseSizes[ 0 ];
00597       for( Field *it = values; it != col0End; ++it )
00598         *it *= Field( 1 ) /  Field( int(dimDomain) );   
00599       Field *row0 = values;
00600 
00601       for( unsigned int k = 1; k <= order; ++k )
00602       {
00603         const Field factor = (Field( 1 ) / Field( k + dimDomain ));
00604 
00605         Field *const row1 = values+offsets[ k-1 ];
00606         Field *const col0End = row1 + baseSizes[ k ];
00607         Field *it = row1;
00608         for( ; it != col0End; ++it )
00609           *it *= factor;
00610         for( unsigned int i = 1; i <= k; ++i )
00611         {
00612           Field *const end = it + baseSizes[ k-i ];
00613           assert( (unsigned int)(end - values) <= offsets[ k ] );
00614           for( ; it != end; ++row0, ++it )
00615             *it = (*row0) * (Field( i ) * factor);
00616         }
00617         row0 = row1;
00618       }
00619     }
00620   };
00621 
00622 
00623 
00624   // MonomialBasis
00625   // -------------
00626 
00627   template< class Topology, class F >
00628   class MonomialBasis
00629   : public MonomialBasisImpl< Topology, F >
00630   {
00631     typedef MonomialBasis< Topology, F > This;
00632     typedef MonomialBasisImpl< Topology, F > Base;
00633 
00634   public:
00635     static const unsigned int dimension = Base::dimDomain;
00636     static const unsigned int dimRange = 1;
00637 
00638     typedef typename Base::Field Field;
00639 
00640     typedef typename Base::DomainVector DomainVector;
00641 
00642     typedef Dune::FieldVector<Field,dimRange> RangeVector;
00643 
00644     typedef MonomialBasisSize<Topology> Size;
00645 
00646     MonomialBasis (unsigned int order)
00647     : Base(),
00648       order_(order),
00649       size_(Size::instance())
00650     {
00651       assert(order<=1024); // avoid wrapping of unsigned int (0-1) order=1024 is quite hight...)
00652     }
00653 
00654     const unsigned int *sizes ( unsigned int order ) const
00655     {
00656       size_.computeSizes( order );
00657       return size_.numBaseFunctions_;
00658     }
00659 
00660     const unsigned int *sizes () const
00661     {
00662       return sizes( order_ );
00663     }
00664 
00665     const unsigned int size () const
00666     {
00667       size_.computeSizes( order_ );
00668       return size_( order_ );
00669     }
00670 
00671     const unsigned int derivSize ( const unsigned int deriv ) const
00672     {
00673       typedef typename GenericGeometry::SimplexTopology< dimension >::type SimplexTopology;
00674       MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
00675       return MonomialBasisSize< SimplexTopology >::instance()( deriv );
00676     }
00677 
00678     const unsigned int order () const
00679     {
00680       return order_ ;
00681     }
00682 
00683     const unsigned int topologyId ( ) const
00684     {
00685       return Topology::id;
00686     }
00687 
00688     void evaluate ( const unsigned int deriv, const DomainVector &x,
00689                     Field *const values ) const
00690     {
00691       Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
00692     }
00693 
00694     template <unsigned int deriv>
00695     void evaluate ( const DomainVector &x,
00696                     Field *const values ) const
00697     {
00698       evaluate( deriv, x, values );
00699     }
00700 
00701     template<unsigned int deriv, class Vector >
00702     void evaluate ( const DomainVector &x,
00703                     Vector &values ) const
00704     {
00705       evaluate<deriv>(x,&(values[0]));
00706     }
00707     template<unsigned int deriv, DerivativeLayout layout >
00708     void evaluate ( const DomainVector &x,
00709                     Derivatives<Field,dimension,1,deriv,layout> *values ) const
00710     {
00711       evaluate<deriv>(x,&(values->block()));
00712     }
00713     template< unsigned int deriv >
00714     void evaluate ( const DomainVector &x,
00715                     FieldVector<Field,Derivatives<Field,dimension,1,deriv,value>::size> *values ) const
00716     {
00717       evaluate(0,x,&(values[0][0]));
00718     }
00719 
00720     template<class Vector >
00721     void evaluate ( const DomainVector &x,
00722                     Vector &values ) const
00723     {
00724       evaluate<0>(x,&(values[0]));
00725     }
00726 
00727     template< class DVector, class RVector >
00728     void evaluate ( const DVector &x, RVector &values ) const
00729     {
00730       assert( DVector::dimension == dimension);
00731       DomainVector bx;
00732       for( int d = 0; d < dimension; ++d )
00733         field_cast( x[ d ], bx[ d ] );
00734       evaluate<0>( bx, values );
00735     }
00736 
00737     void integrate ( Field *const values ) const
00738     {
00739       Base::integrate( order_, sizes( order_ ), values );
00740     }
00741     template <class Vector>
00742     void integrate ( Vector &values ) const
00743     {
00744       integrate( &(values[ 0 ]) );
00745     }
00746   private:
00747     MonomialBasis(const This&);
00748     This& operator=(const This&);
00749     unsigned int order_;
00750     Size &size_; 
00751   };
00752 
00753 
00754 
00755   // StdMonomialBasis
00756   // ----------------
00757 
00758   template< int dim,class F >
00759   class StandardMonomialBasis
00760   : public MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F >
00761   {
00762     typedef StandardMonomialBasis< dim, F > This;
00763     typedef MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F > Base;
00764 
00765   public:
00766     typedef typename GenericGeometry::SimplexTopology< dim >::type Topology;
00767     static const int dimension = dim;
00768 
00769     StandardMonomialBasis ( unsigned int order )
00770     : Base( order )
00771     {}
00772   };
00773 
00774 
00775 
00776   // StandardBiMonomialBasis
00777   // -----------------------
00778 
00779   template< int dim, class F >
00780   class StandardBiMonomialBasis
00781   : public MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F >
00782   {
00783     typedef StandardBiMonomialBasis< dim, F > This;
00784     typedef MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F > Base;
00785 
00786   public:
00787     typedef typename GenericGeometry::CubeTopology< dim >::type Topology;
00788     static const int dimension = dim;
00789 
00790     StandardBiMonomialBasis ( unsigned int order )
00791     : Base( order )
00792     {}
00793   };
00794 
00795   // -----------------------------------------------------------
00796   // -----------------------------------------------------------
00797   // VirtualMonomialBasis
00798   // -------------------
00799 
00800   template< int dim, class F >
00801   class VirtualMonomialBasis 
00802   {
00803     typedef VirtualMonomialBasis< dim, F > This;
00804 
00805   public:
00806     typedef F Field;
00807     typedef F StorageField;
00808     static const int dimension = dim;
00809     static const unsigned int dimRange = 1;
00810 
00811     typedef FieldVector<Field,dimension> DomainVector;
00812     typedef FieldVector<Field,dimRange> RangeVector;
00813 
00814     explicit VirtualMonomialBasis(unsigned int topologyId,
00815                                   unsigned int order)
00816       : order_(order), topologyId_(topologyId) {}
00817 
00818     virtual ~VirtualMonomialBasis() {}
00819 
00820     virtual const unsigned int *sizes ( ) const = 0;
00821 
00822     const unsigned int size ( ) const
00823     {
00824       return sizes( )[ order_ ];
00825     }
00826 
00827     const unsigned int order () const
00828     {
00829       return order_;
00830     }
00831 
00832     const unsigned int topologyId ( ) const
00833     {
00834       return topologyId_;
00835     }
00836 
00837     virtual void evaluate ( const unsigned int deriv, const DomainVector &x,
00838                             Field *const values ) const = 0;
00839     template < unsigned int deriv >
00840     void evaluate ( const DomainVector &x,
00841                     Field *const values ) const 
00842     {
00843       evaluate( deriv, x, values );
00844     }
00845     template < unsigned int deriv, int size >
00846     void evaluate ( const DomainVector &x,
00847                     Dune::FieldVector<Field,size> *const values ) const
00848     {
00849       evaluate( deriv, x, &(values[0][0]) );
00850     }
00851     template<unsigned int deriv, DerivativeLayout layout >
00852     void evaluate ( const DomainVector &x,
00853                     Derivatives<Field,dimension,1,deriv,layout> *values ) const
00854     {
00855       evaluate<deriv>(x,&(values->block()));
00856     }
00857     template <unsigned int deriv, class Vector>
00858     void evaluate ( const DomainVector &x,
00859                     Vector &values ) const
00860     {
00861       evaluate<deriv>( x, &(values[ 0 ]) );
00862     }
00863     template< class Vector >
00864     void evaluate ( const DomainVector &x,
00865                     Vector &values ) const
00866     {
00867       evaluate<0>(x,values);
00868     }
00869     template< class DVector, class RVector >
00870     void evaluate ( const DVector &x, RVector &values ) const
00871     {
00872       assert( DVector::dimension == dimension);
00873       DomainVector bx;
00874       for( int d = 0; d < dimension; ++d )
00875         field_cast( x[ d ], bx[ d ] );
00876       evaluate<0>( bx, values );
00877     }
00878     template< unsigned int deriv, class DVector, class RVector >
00879     void evaluate ( const DVector &x, RVector &values ) const
00880     {
00881       assert( DVector::dimension == dimension);
00882       DomainVector bx;
00883       for( int d = 0; d < dimension; ++d )
00884         field_cast( x[ d ], bx[ d ] );
00885       evaluate<deriv>( bx, values );
00886     }
00887 
00888     virtual void integrate ( Field *const values ) const = 0;
00889     template <class Vector>
00890     void integrate ( Vector &values ) const
00891     {
00892       integrate( &(values[ 0 ]) );
00893     }
00894   protected:
00895     unsigned int order_;
00896     unsigned int topologyId_;
00897   };
00898 
00899   template< class Topology, class F >
00900   class VirtualMonomialBasisImpl 
00901   : public VirtualMonomialBasis< Topology::dimension, F >
00902   {
00903     typedef VirtualMonomialBasis< Topology::dimension, F > Base;
00904     typedef VirtualMonomialBasisImpl< Topology, F > This;
00905 
00906   public:
00907     typedef typename Base::Field Field;
00908     typedef typename Base::DomainVector DomainVector;
00909 
00910     VirtualMonomialBasisImpl(unsigned int order)
00911     : Base(Topology::id,order), basis_(order) 
00912     {}
00913 
00914     const unsigned int *sizes ( ) const
00915     {
00916       return basis_.sizes(order_);
00917     }
00918 
00919     void evaluate ( const unsigned int deriv, const DomainVector &x,
00920                     Field *const values ) const
00921     {
00922       basis_.evaluate(deriv,x,values);
00923     }
00924 
00925     void integrate ( Field *const values ) const
00926     {
00927       basis_.integrate(values);
00928     }
00929 
00930   private:
00931     MonomialBasis<Topology,Field> basis_;
00932     using Base::order_;
00933   };
00934 
00935   // MonomialBasisFactory
00936   // --------------------
00937   
00938   template< int dim, class F >
00939   struct MonomialBasisFactory;
00940 
00941   template< int dim, class F >
00942   struct MonomialBasisFactoryTraits
00943   {
00944     static const unsigned int dimension = dim;
00945     typedef unsigned int Key;
00946     typedef const VirtualMonomialBasis< dimension, F > Object;
00947     typedef MonomialBasisFactory<dim,F> Factory;
00948   };
00949 
00950   template< int dim, class F >
00951   struct MonomialBasisFactory : 
00952     public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
00953   {
00954     static const unsigned int dimension = dim;
00955     typedef F StorageField;
00956     typedef MonomialBasisFactoryTraits<dim,F> Traits;
00957 
00958     typedef typename Traits::Key Key;
00959     typedef typename Traits::Object Object;
00960 
00961     template < int dd, class FF >
00962     struct EvaluationBasisFactory
00963     {
00964       typedef MonomialBasisFactory<dd,FF> Type;
00965     };
00966 
00967     template< class Topology >
00968     static Object* createObject ( const Key &order )
00969     {
00970       return new VirtualMonomialBasisImpl< Topology, StorageField >( order );
00971     }
00972   };
00973 
00974 
00975 
00976   // MonomialBasisProvider
00977   // ---------------------
00978 
00979   template< int dim, class SF >
00980   struct MonomialBasisProvider 
00981   : public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
00982   {
00983     static const unsigned int dimension = dim;
00984     typedef SF StorageField;
00985     template < int dd, class FF >
00986     struct EvaluationBasisFactory
00987     {
00988       typedef MonomialBasisProvider<dd,FF> Type;
00989     };
00990   };
00991 
00992 }
00993 
00994 #endif