dune-localfunctions  2.2.0
polynomialbasis.hh
Go to the documentation of this file.
00001 #ifndef DUNE_POLYNOMIALBASIS_HH
00002 #define DUNE_POLYNOMIALBASIS_HH
00003 
00004 #include <fstream>
00005 
00006 #include <dune/common/fmatrix.hh>
00007 
00008 #include <dune/localfunctions/common/localbasis.hh>
00009 
00010 #include <dune/localfunctions/utility/coeffmatrix.hh>
00011 #include <dune/localfunctions/utility/monomialbasis.hh>
00012 #include <dune/localfunctions/utility/multiindex.hh>
00013 #include <dune/localfunctions/utility/basisevaluator.hh>
00014 
00015 namespace Dune
00016 {
00017 
00018   // PolynomialBasis
00019   // ---------------
00020 
00058   template< class Eval, class CM, class D=double, class R=double >
00059   class PolynomialBasis 
00060   {
00061     typedef PolynomialBasis< Eval, CM > This;
00062     typedef Eval Evaluator;
00063 
00064   public:
00065     typedef CM CoefficientMatrix;
00066 
00067     typedef typename CoefficientMatrix::Field StorageField;
00068 
00069     static const unsigned int dimension = Evaluator::dimension;
00070     static const unsigned int dimRange = Evaluator::dimRange*CoefficientMatrix::blockSize;
00071     typedef LocalBasisTraits<D,dimension,FieldVector<D,dimension>,
00072                                R,dimRange,FieldVector<R,dimRange>,
00073                                FieldMatrix<R,dimRange,dimension> > Traits;
00074     typedef typename Evaluator::Basis Basis;
00075     typedef typename Evaluator::DomainVector DomainVector;
00076 
00077     PolynomialBasis (const Basis &basis, 
00078                      const CoefficientMatrix &coeffMatrix,
00079                      unsigned int size)
00080     : basis_(basis),
00081       coeffMatrix_(&coeffMatrix),
00082       eval_(basis), 
00083       order_(basis.order()),
00084       size_(size)
00085     {
00086       // assert(coeffMatrix_);
00087       // assert(size_ <= coeffMatrix.size()); // !!!
00088     }
00089 
00090     const Basis &basis () const
00091     {
00092       return basis_;
00093     }
00094 
00095     const CoefficientMatrix &matrix () const
00096     {
00097       return *coeffMatrix_;
00098     }
00099 
00100     const unsigned int order () const 
00101     {
00102       return order_;
00103     }
00104 
00105     const unsigned int size () const
00106     {
00107       return size_;
00108     }
00109 
00111     void evaluateFunction (const typename Traits::DomainType& x,
00112                            std::vector<typename Traits::RangeType>& out) const
00113     { 
00114       out.resize(size());
00115       evaluate(x,out);
00116     }
00117 
00119     void evaluateJacobian (const typename Traits::DomainType& x,         // position
00120                            std::vector<typename Traits::JacobianType>& out) const      // return value
00121     {  
00122       out.resize(size());
00123       jacobian(x,out);
00124     }
00125 
00126     template< unsigned int deriv, class F >
00127     void evaluate ( const DomainVector &x, F *values ) const
00128     {
00129       coeffMatrix_->mult( eval_.template evaluate<deriv>( x ), size(), values);
00130     }
00131     template< unsigned int deriv, class DVector, class F >
00132     void evaluate ( const DVector &x, F *values ) const
00133     {
00134       assert( DVector::dimension == dimension);
00135       DomainVector bx;
00136       for( int d = 0; d < dimension; ++d )
00137         field_cast( x[ d ], bx[ d ] );
00138       evaluate<deriv>( bx, values );
00139     }
00140 
00141     template <bool dummy,class DVector>
00142     struct Convert
00143     {
00144       static DomainVector apply( const DVector &x )
00145       {
00146         assert( DVector::dimension == dimension);
00147         DomainVector bx;
00148         for( unsigned int d = 0; d < dimension; ++d )
00149           field_cast( x[ d ], bx[ d ] );
00150         return bx;
00151       }
00152     };
00153     template <bool dummy>
00154     struct Convert<dummy,DomainVector>
00155     {
00156       static const DomainVector &apply( const DomainVector &x )
00157       {
00158         return x;
00159       }
00160     };
00161     template< unsigned int deriv, class DVector, class RVector >
00162     void evaluate ( const DVector &x, RVector &values ) const
00163     {
00164       assert(values.size()>=size());
00165       const DomainVector &bx = Convert<true,DVector>::apply(x);
00166       coeffMatrix_->mult( eval_.template evaluate<deriv>( bx ), values );
00167     }
00168 
00169     template <class Fy>
00170     void evaluate ( const DomainVector &x, std::vector<FieldVector<Fy,dimRange> > &values ) const
00171     {
00172       evaluate<0>(x,values);
00173     }
00174     template< class DVector, class RVector >
00175     void evaluate ( const DVector &x, RVector &values ) const
00176     {
00177       assert( DVector::dimension == dimension);
00178       DomainVector bx;
00179       for( unsigned int d = 0; d < dimension; ++d )
00180         field_cast( x[ d ], bx[ d ] );
00181       evaluate<0>( bx, values );
00182     }
00183 
00184     template< unsigned int deriv, class Vector >
00185     void evaluateSingle ( const DomainVector &x, Vector &values ) const
00186     {
00187       assert(values.size()>=size());
00188       coeffMatrix_->template mult<deriv>( eval_.template evaluate<deriv>( x ), values );
00189     }
00190     template< unsigned int deriv, class Fy >
00191     void evaluateSingle ( const DomainVector &x, 
00192                           std::vector< FieldVector<FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size>,dimRange> > &values) const
00193     {
00194       evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
00195     }
00196     template< unsigned int deriv, class Fy >
00197     void evaluateSingle ( const DomainVector &x, 
00198                           std::vector< FieldVector<LFETensor<Fy,dimension,deriv>,dimRange> > &values) const
00199     {
00200       evaluateSingle<deriv>(x,reinterpret_cast<std::vector< FieldVector<Fy,LFETensor<Fy,dimension,deriv>::size*dimRange> >&>(values));
00201     }
00202 
00203     template <class Fy>
00204     void jacobian ( const DomainVector &x, std::vector<FieldMatrix<Fy,dimRange,dimension> > &values ) const
00205     {
00206       assert(values.size()>=size());
00207       evaluateSingle<1>(x,reinterpret_cast<std::vector<FieldVector<Fy,dimRange*dimension> >&>(values));
00208     }
00209     template< class DVector, class RVector >
00210     void jacobian ( const DVector &x, RVector &values ) const
00211     {
00212       assert( DVector::dimension == dimension);
00213       DomainVector bx;
00214       for( unsigned int d = 0; d < dimension; ++d )
00215         field_cast( x[ d ], bx[ d ] );
00216       jacobian( bx, values );
00217     }
00218 
00219     template <class Fy>
00220     void integrate ( std::vector<Fy> &values ) const
00221     {
00222       assert(values.size()>=size());
00223       coeffMatrix_->mult( eval_.template integrate(), values );
00224     }
00225 
00226   protected:
00227     PolynomialBasis(const PolynomialBasis &other) 
00228       : basis_(other.basis_),
00229         coeffMatrix_(other.coeffMatrix_),
00230         eval_(basis_),
00231         order_(basis_.order()),
00232         size_(other.size_)
00233     {
00234     }
00235     PolynomialBasis &operator=(const PolynomialBasis&);
00236     const Basis &basis_;
00237     const CoefficientMatrix* coeffMatrix_;
00238     mutable Evaluator eval_;
00239     unsigned int order_,size_;
00240   };
00241 
00248   template< class Eval, class CM = SparseCoeffMatrix<typename Eval::Field,Eval::dimRange>,
00249             class D=double, class R=double>
00250   class PolynomialBasisWithMatrix 
00251   : public PolynomialBasis< Eval, CM, D, R >
00252   {
00253   public:
00254     typedef CM CoefficientMatrix;
00255 
00256   private:
00257     typedef Eval Evaluator;
00258 
00259     typedef PolynomialBasisWithMatrix< Evaluator, CM > This;
00260     typedef PolynomialBasis<Evaluator,CM> Base;
00261 
00262   public:
00263     typedef typename Base::Basis Basis;
00264 
00265     PolynomialBasisWithMatrix (const Basis &basis)
00266     : Base(basis,coeffMatrix_,0)
00267     {}
00268 
00269     template <class Matrix>
00270     void fill(const Matrix& matrix)
00271     {
00272       coeffMatrix_.fill(matrix);
00273       this->size_ = coeffMatrix_.size();
00274     }
00275     template <class Matrix>
00276     void fill(const Matrix& matrix,int size)
00277     {
00278       coeffMatrix_.fill(matrix);
00279       assert(size<=coeffMatrix_.size());
00280       this->size_ = size;
00281     }
00282 
00283   private:
00284     PolynomialBasisWithMatrix(const PolynomialBasisWithMatrix &);
00285     PolynomialBasisWithMatrix &operator=(const PolynomialBasisWithMatrix &);
00286     CoefficientMatrix coeffMatrix_;
00287   };
00288 }
00289 #endif // DUNE_POLYNOMIALBASIS_HH
00290