dune-localfunctions
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 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