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