dune-localfunctions  2.2.0
basisevaluator.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 
00004 #ifndef DUNE_BASISEVALUATOR_HH
00005 #define DUNE_BASISEVALUATOR_HH
00006 
00007 #include <vector>
00008 
00009 #include <dune/common/fmatrix.hh>
00010 #include <dune/common/fvector.hh>
00011 #include <dune/common/typetraits.hh>
00012 
00013 #include <dune/geometry/genericgeometry/topologytypes.hh>
00014 
00015 #include <dune/localfunctions/utility/field.hh>
00016 #include <dune/localfunctions/utility/multiindex.hh>
00017 #include <dune/localfunctions/utility/tensor.hh>
00018 
00019 namespace Dune
00020 {
00021   /*******************************************
00022    * Should be removed as soon as the Tensor
00023    * classes have been revisited. See remarks
00024    * in tensor.hh (also hold true here).
00025    *******************************************/
00026  
00027 
00028   template <class B>
00029   struct MonomialEvaluator
00030   {
00031     typedef B Basis;
00032     typedef typename Basis::Field Field;
00033     typedef typename Basis::DomainVector DomainVector;
00034     static const int dimension = Basis::dimension;
00035     static const int dimRange = Basis::dimRange;
00036 
00037     typedef std::vector<Field> Container;
00038 
00039     template< class Deriv >
00040     struct BaseIterator;
00041 
00042     template <unsigned int deriv>
00043     struct Iterator 
00044     {
00045       typedef BaseIterator<Derivatives<Field,dimension,dimRange,deriv,derivative> > All;
00046       typedef BaseIterator<Derivatives<Field,dimension,1,0,value> > Integrate;
00047     };
00048 
00049     unsigned int size() const 
00050     {
00051       return size_;
00052     }
00053 
00054     protected:
00055     MonomialEvaluator(const Basis &basis,unsigned int order,unsigned int size) 
00056     : basis_(basis),
00057       order_(order),
00058       size_(size),
00059       container_(0)
00060     {
00061     }
00062     template <int deriv>
00063     void resize() 
00064     {
00065       const int totalSize = Derivatives<Field,dimension,dimRange,deriv,derivative>::size*size_;
00066       container_.resize(totalSize);
00067     }
00068     MonomialEvaluator(const MonomialEvaluator&);
00069     const Basis &basis_;
00070     unsigned int order_,size_;
00071     Container container_;
00072   };
00073 
00074 
00075   template< class B >
00076   template< class Deriv >
00077   struct MonomialEvaluator< B >::BaseIterator
00078   {
00079     typedef Deriv Derivatives;
00080     typedef typename Deriv::Field Field;
00081     static const unsigned int blockSize = Deriv::size;
00082     typedef Dune::FieldVector<Field,blockSize> Block;
00083     static const DerivativeLayout layout = Deriv::layout;
00084     static const unsigned int dimDomain = Deriv::dimDomain;
00085     static const unsigned int dimRange = Deriv::dimRange;
00086 
00087     typedef std::vector<Field> Container;
00088     typedef typename Container::iterator CIter;
00089 
00090     explicit BaseIterator ( Container &container ) 
00091     : pos_( container.begin() ),
00092       end_( container.end() )
00093     {}
00094 
00095     const Deriv &operator*() const 
00096     {
00097       assert(!done());
00098       return reinterpret_cast<const Deriv&>(*pos_);
00099     }
00100 
00101     const Deriv *operator->() const 
00102     {
00103       return &(operator*());
00104     }
00105 
00106     bool done () const
00107     {
00108       return pos_ == end_;
00109     }
00110 
00111     BaseIterator &operator++ () 
00112     {
00113       pos_ += blockSize;
00114       return *this;
00115     }
00116 
00117     BaseIterator &operator+= ( unsigned int skip )
00118     {
00119       pos_ += skip*blockSize;
00120       return *this;
00121     }
00122 
00123   private:
00124     CIter pos_;
00125     const CIter end_;
00126   };
00127 
00128   template< class B >
00129   struct StandardEvaluator
00130   : public MonomialEvaluator< B >
00131   {
00132     typedef B Basis;
00133     typedef typename Basis::Field Field;
00134     typedef typename Basis::DomainVector DomainVector;
00135     typedef std::vector<Field> Container;
00136     static const int dimension = Basis::dimension;
00137     static const int dimRange = Basis::dimRange;
00138     typedef MonomialEvaluator<B> Base;
00139 
00140     template <unsigned int deriv>
00141     struct Iterator : public Base::template Iterator<deriv>
00142     {
00143     };
00144 
00145     StandardEvaluator(const Basis &basis) 
00146     : Base(basis,basis.order(),basis.size())
00147     {
00148     }
00149     template <unsigned int deriv,class DVector>
00150     typename Iterator<deriv>::All evaluate(const DVector &x) 
00151     {
00152       Base::template resize<deriv>();
00153       basis_.template evaluate<deriv>(x,&(container_[0]));
00154       return typename Iterator<deriv>::All(container_);
00155     }
00156     typename Iterator<0>::Integrate integrate() 
00157     {
00158       Base::template resize<0>();
00159       basis_.integrate(&(container_[0]));
00160       return typename Iterator<0>::Integrate(container_);
00161     }
00162 
00163   protected:
00164     StandardEvaluator ( const Basis &basis, unsigned int size )
00165     : Base( basis, basis.order(), size )
00166     {}
00167 
00168   private:
00169     StandardEvaluator(const StandardEvaluator&);
00170     using Base::basis_;
00171     using Base::container_;
00172   };
00173 
00174 #if 0 // OLD OLD
00175   template< class B, class Fill >
00176   struct VecEvaluator
00177   : public StandardEvaluator< B >
00178   {
00179     typedef B Basis;
00180     typedef typename Basis::Field Field;
00181     static const int dimension = Basis::dimension;
00182     static const int dimRange = Basis::dimRange*Fill::dimRange;
00183     typedef typename Basis::DomainVector DomainVector;
00184     typedef std::vector<Field> Container;
00185     typedef StandardEvaluator<B> Base;
00186 
00187     template <unsigned int deriv>
00188     struct Iterator 
00189     {
00190       typedef typename Base::template BaseIterator<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> > All;
00191     };
00192 
00193     VecEvaluator ( const Basis &basis, const Fill &fill )
00194     : Base( basis, basis.size() ),
00195       fill_( fill ),
00196       size_( basis.size()*dimRange )
00197     {
00198     }
00199     template <unsigned int deriv>
00200     typename Iterator<deriv>::All evaluate(const DomainVector &x) 
00201     {
00202       resize< deriv >();
00203       fill_.template apply<deriv>( x,Base::template evaluate<deriv>(x), vecContainer_ );
00204       std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >& derivContainer =
00205          reinterpret_cast<std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >&>(vecContainer_);
00206       return typename Iterator<deriv>::All(derivContainer);
00207     }
00208     template <unsigned int deriv,class DVector>
00209     typename Iterator<deriv>::All evaluate(const DVector &x) 
00210     {
00211       resize< deriv >();
00212       fill_.template apply<deriv>( x,Base::template evaluate<deriv>(x), vecContainer_ );
00213       std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >& derivContainer =
00214          reinterpret_cast<std::vector<Derivatives<Field,dimension,dimRange,deriv,Fill::layout> >&>(vecContainer_);
00215       return typename Iterator<deriv>::All(derivContainer);
00216     }
00217     unsigned int size() const
00218     {
00219       return size_;
00220     }
00221 
00222   protected:
00223     VecEvaluator ( const Basis &basis, const Fill &fill, unsigned int size )
00224     : Base( basis, basis.size() ),
00225       fill_( fill ),
00226       size_( size )
00227     {
00228       resize< 2 >();
00229     }
00230 
00231     template <int deriv>
00232     void resize() 
00233     {
00234       const int totalSize = Derivatives<Field,dimension,dimRange,deriv,derivative>::size*size_;
00235       vecContainer_.resize(totalSize);
00236     }
00237 
00238     VecEvaluator(const VecEvaluator&);
00239 
00240     Container vecContainer_;
00241     const Fill &fill_;
00242     unsigned int size_;
00243   };
00244 
00245   template <int dimR,DerivativeLayout layout>
00246   struct DiagonalFill;
00247 
00248   template <int dimR>
00249   struct DiagonalFill<dimR,derivative>
00250   {
00251     static const DerivativeLayout layout = derivative;
00252     static const int dimRange = dimR;
00253     template <int deriv, class Domain, class Iter,class Field>
00254     void apply(const Domain &x,
00255                Iter iter,std::vector<Field> &vecContainer) const
00256     {
00257       typedef std::vector<Field> Container;
00258       typename Container::iterator vecIter = vecContainer.begin();
00259       for ( ; !iter.done(); ++iter)
00260       {
00261         const typename Iter::Block &block = iter->block();
00262         for (int r1=0;r1<dimR;++r1) 
00263         {
00264           unsigned int b = 0;
00265           apply<Field>(r1,x,block,b,vecIter);
00266         }
00267       }
00268     }
00269     template <class Field, class Domain, class Block,class VecIter>
00270     void apply(int r1, const Domain &x,
00271                const Block &block,unsigned int &b, 
00272                VecIter &vecIter) const
00273     {
00274       unsigned int bStart = b;
00275       unsigned int bEnd = b+Block::size;
00276       apply<Field>(r1,x,block,bStart,bEnd,vecIter);
00277       b=bEnd;
00278     }
00279     template <class Field, class Domain, class Block,class VecIter>
00280     void apply(int r1, const Domain &x,const Block &block,
00281                unsigned int bStart, unsigned int bEnd,
00282                VecIter &vecIter) const
00283     {
00284       for (int r2=0;r2<dimR;++r2) 
00285       {
00286         for (unsigned int bb=bStart;bb<bEnd;++bb)
00287         {
00288           *vecIter = (r1==r2?block[bb]:Field(0));
00289           ++vecIter;
00290         }
00291       }
00292     }
00293   };
00294   template <int dimR>
00295   struct DiagonalFill<dimR,value>
00296   {
00297     static const DerivativeLayout layout = value;
00298     static const int dimRange = dimR;
00299     template <int deriv, class Domain, class Iter,class Field>
00300     void apply(const Domain &x,
00301                Iter iter,std::vector<Field> &vecContainer) const
00302     {
00303       typedef std::vector<Field> Container;
00304       typename Container::iterator vecIter = vecContainer.begin();
00305       for ( ; !iter.done(); ++iter)
00306       {
00307         const typename Iter::Block &block = iter->block();
00308         for (int r1=0;r1<dimR;++r1) 
00309         {
00310           unsigned int b = 0;
00311           apply<Field>(integral_constant<int,deriv>(),r1,x,block,b,vecIter);
00312         }
00313       }
00314     }
00315     template <class Field, class Domain, class Block,class VecIter,int deriv>
00316     void apply(const integral_constat<int,deriv>&, int r1, const Domain &x,
00317                const Block &block,unsigned int &b, 
00318                VecIter &vecIter) const
00319     {
00320       apply<Field>(integral_constant<int,deriv-1>(),r1,x,block,b,vecIter);
00321       unsigned int bStart = b;
00322       unsigned int bEnd = b+LFETensor<Field,Domain::dimension,deriv>::size;
00323       apply<Field>(r1,x,block,bStart,bEnd,vecIter);
00324       b=bEnd;
00325     }
00326     template <class Field, class Domain, class Block,class VecIter>
00327     void apply(const integral_constant<int,0>&, int r1, const Domain &x,
00328                const Block &block,unsigned int &b, 
00329                VecIter &vecIter) const
00330     {
00331       apply<Field>(r1,x,block,b,b+1,vecIter);
00332       ++b;
00333     }
00334     template <class Field, class Domain, class Block,class VecIter>
00335     void apply(int r1, const Domain &x,const Block &block,
00336                unsigned int bStart, unsigned int bEnd,
00337                VecIter &vecIter) const
00338     {
00339       for (int r2=0;r2<dimR;++r2) 
00340       {
00341         for (unsigned int bb=bStart;bb<bEnd;++bb)
00342         {
00343           *vecIter = (r1==r2?block[bb]:Field(0));
00344           ++vecIter;
00345         }
00346       }
00347     }
00348   };
00349 
00350   template <class B,int dimR,DerivativeLayout layout>
00351   struct VectorialEvaluator 
00352   : public VecEvaluator<B,DiagonalFill<dimR,layout> >
00353   {
00354     typedef DiagonalFill<dimR,layout> Fill;
00355     typedef VecEvaluator< B,Fill > Base;
00356     VectorialEvaluator(const B &basis)
00357     : Base(basis,fill_,basis.size()*dimR)
00358     {
00359     }
00360     private:
00361     Fill fill_;
00362   };
00363 #endif // OLD OLD
00364 
00365 }
00366 
00367 #endif