dune-localfunctions  2.2.1
monomialbasis.hh
Go to the documentation of this file.
1 #ifndef DUNE_MONOMIALBASIS_HH
2 #define DUNE_MONOMIALBASIS_HH
3 
4 #include <vector>
5 
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
8 
9 #include <dune/geometry/topologyfactory.hh>
10 #include <dune/geometry/genericgeometry/topologytypes.hh>
11 
15 
16 namespace Dune
17 {
18 /************************************************
19  * Classes for evaluating ''Monomials'' on any order
20  * for all reference element type.
21  * For a simplex topology these are the nomral
22  * monomials for cube topologies the bimonomials.
23  * The construction follows the construction of the
24  * generic geometries using tensor products for
25  * prism generation and duffy transform for pyramid
26  * construction.
27  * A derivative argument can be applied, in which case
28  * all derivatives up to the desired order are
29  * evaluated. Note that in for higher order derivatives
30  * only the ''lower'' part of the symmetric tensor
31  * is evaluated, e.g., passing derivative equal to 2
32  * to the class will provide the vector
33  * (d/dxdx p, d/dxydx p, d/dydy p,
34  * d/dx p, d/dy p, p)
35  * Important:
36  * So far the computation of the derivatives has not
37  * been fully implemented for general pyramid
38  * construction, i.e., in the case where a pyramid is
39  * build over a non simplex base geometry.
40  *
41  * Central classes:
42  * 1) template< class Topology, class F >
43  * class MonomialBasisImpl;
44  * Implementation of the monomial evaluation for
45  * a given topology and field type.
46  * The method evaluate fills a F* vector
47  * 2) template< class Topology, class F >
48  * class MonomialBasis
49  * The base class for the static monomial evaluation
50  * providing addiional evaluate methods including
51  * one taking std::vector<F>.
52  * 3) template< int dim, class F >
53  * class VirtualMonomialBasis
54  * Virtualization of the MonomialBasis.
55  * 4) template< int dim, class F >
56  * struct MonomialBasisFactory;
57  * A factory class for the VirtualMonomialBasis
58  * 5) template< int dim, class F >
59  * struct MonomialBasisProvider
60  * A singleton container for the virtual monomial
61  * basis
62  ************************************************/
63 
64  // Internal Forward Declarations
65  // -----------------------------
66 
67  template< class Topology >
69 
70  template< class Topology, class F >
72 
73 
74 
75  // MonomialBasisSize
76  // -----------------
77 
78  template<>
79  class MonomialBasisSize< GenericGeometry::Point >
80  {
82 
83  public:
84  static This &instance ()
85  {
86  static This _instance;
87  return _instance;
88  }
89 
90  typedef GenericGeometry::Point Topology;
91 
92  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
93  friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
94 
95  mutable unsigned int maxOrder_;
96  // sizes_[ k ]: number of basis functions of exactly order k
97  mutable unsigned int *sizes_;
98  // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
99  mutable unsigned int *numBaseFunctions_;
100 
102  : maxOrder_( 0 ),
103  sizes_( 0 ),
104  numBaseFunctions_( 0 )
105  {
106  computeSizes( 2 );
107  }
108 
110  {
111  delete[] sizes_;
112  delete[] numBaseFunctions_;
113  }
114 
115  unsigned int operator() ( const unsigned int order ) const
116  {
117  return numBaseFunctions_[ order ];
118  }
119 
120  unsigned int maxOrder () const
121  {
122  return maxOrder_;
123  }
124 
125  void computeSizes ( unsigned int order ) const
126  {
127  if (order <= maxOrder_)
128  return;
129 
130  maxOrder_ = order;
131 
132  delete [] sizes_;
133  delete [] numBaseFunctions_;
134  sizes_ = new unsigned int [ order+1 ];
135  numBaseFunctions_ = new unsigned int [ order+1 ];
136 
137  sizes_[ 0 ] = 1;
138  numBaseFunctions_[ 0 ] = 1;
139  for( unsigned int k = 1; k <= order; ++k )
140  {
141  sizes_[ k ] = 0;
142  numBaseFunctions_[ k ] = 1;
143  }
144  }
145  };
146 
147  template< class BaseTopology >
148  class MonomialBasisSize< GenericGeometry::Prism< BaseTopology > >
149  {
151 
152  public:
153  static This &instance ()
154  {
155  static This _instance;
156  return _instance;
157  }
158 
159  typedef GenericGeometry::Prism< BaseTopology > Topology;
160 
161  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
162  friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
163 
164  mutable unsigned int maxOrder_;
165  // sizes_[ k ]: number of basis functions of exactly order k
166  mutable unsigned int *sizes_;
167  // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
168  mutable unsigned int *numBaseFunctions_;
169 
171  : maxOrder_( 0 ),
172  sizes_( 0 ),
173  numBaseFunctions_( 0 )
174  {
175  computeSizes( 2 );
176  }
177 
179  {
180  delete[] sizes_;
181  delete[] numBaseFunctions_;
182  }
183 
184  unsigned int operator() ( const unsigned int order ) const
185  {
186  return numBaseFunctions_[ order ];
187  }
188 
189  unsigned int maxOrder() const
190  {
191  return maxOrder_;
192  }
193 
194  void computeSizes ( unsigned int order ) const
195  {
196  if (order <= maxOrder_)
197  return;
198 
199  maxOrder_ = order;
200 
201  delete[] sizes_;
202  delete[] numBaseFunctions_;
203  sizes_ = new unsigned int[ order+1 ];
204  numBaseFunctions_ = new unsigned int[ order+1 ];
205 
208  baseBasis.computeSizes( order );
209  const unsigned int *const baseSizes = baseBasis.sizes_;
210  const unsigned int *const baseNBF = baseBasis.numBaseFunctions_;
211 
212  sizes_[ 0 ] = 1;
213  numBaseFunctions_[ 0 ] = 1;
214  for( unsigned int k = 1; k <= order; ++k )
215  {
216  sizes_[ k ] = baseNBF[ k ] + k*baseSizes[ k ];
217  numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
218  }
219  }
220  };
221 
222  template< class BaseTopology >
223  class MonomialBasisSize< GenericGeometry::Pyramid< BaseTopology > >
224  {
226 
227  public:
228  static This &instance ()
229  {
230  static This _instance;
231  return _instance;
232  }
233 
234  typedef GenericGeometry::Pyramid< BaseTopology > Topology;
235 
236  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
237  friend class MonomialBasisSize< GenericGeometry::Pyramid< Topology > >;
238 
239  mutable unsigned int maxOrder_;
240  // sizes_[ k ]: number of basis functions of exactly order k
241  mutable unsigned int *sizes_;
242  // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
243  mutable unsigned int *numBaseFunctions_;
244 
246  : maxOrder_( 0 ),
247  sizes_( 0 ),
248  numBaseFunctions_( 0 )
249  {
250  computeSizes( 2 );
251  }
252 
254  {
255  delete[] sizes_;
256  delete[] numBaseFunctions_;
257  }
258 
259  unsigned int operator() ( const unsigned int order ) const
260  {
261  return numBaseFunctions_[ order ];
262  }
263 
264  unsigned int maxOrder() const
265  {
266  return maxOrder_;
267  }
268 
269  void computeSizes ( unsigned int order ) const
270  {
271  if (order <= maxOrder_)
272  return;
273 
274  maxOrder_ = order;
275 
276  delete[] sizes_;
277  delete[] numBaseFunctions_;
278  sizes_ = new unsigned int[ order+1 ];
279  numBaseFunctions_ = new unsigned int[ order+1 ];
280 
283 
284  baseBasis.computeSizes( order );
285 
286  const unsigned int *const baseNBF = baseBasis.numBaseFunctions_;
287  sizes_[ 0 ] = 1;
288  numBaseFunctions_[ 0 ] = 1;
289  for( unsigned int k = 1; k <= order; ++k )
290  {
291  sizes_[ k ] = baseNBF[ k ];
292  numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
293  }
294  }
295  };
296 
297 
298 
299  // MonomialBasisHelper
300  // -------------------
301 
302 
303  template< int mydim, int dim, class F >
305  {
308 
309  static void copy ( const unsigned int deriv, F *&wit, F *&rit,
310  const unsigned int numBaseFunctions, const F &z )
311  {
312  // n(d,k) = size<k>[d];
313  MySize &mySize = MySize::instance();
314  Size &size = Size::instance();
315 
316  const F *const rend = rit + size( deriv )*numBaseFunctions;
317  for( ; rit != rend; )
318  {
319  F *prit = rit;
320 
321  *wit = z * *rit;
322  ++rit, ++wit;
323 
324  for( unsigned d = 1; d <= deriv; ++d )
325  {
326  #ifndef NDEBUG
327  const F *const derivEnd = rit + mySize.sizes_[ d ];
328  #endif
329  const F *const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
330  for( ; rit != drend ; ++rit, ++wit )
331  *wit = z * *rit;
332  for (unsigned int j=1;j<d;++j)
333  {
334  const F *const drend = rit + mySize.sizes_[ d-j ] - mySize.sizes_[ d-j-1 ];
335  for( ; rit != drend ; ++prit, ++rit, ++wit )
336  *wit = F(j) * *prit + z * *rit;
337  }
338  *wit = F(d) * *prit + z * *rit;
339  ++prit, ++rit, ++wit;
340  assert(derivEnd == rit);
341  rit += size.sizes_[d] - mySize.sizes_[d];
342  prit += size.sizes_[d-1] - mySize.sizes_[d-1];
343  const F *const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
344  for ( ; wit != emptyWitEnd; ++wit )
345  *wit = Zero<F>();
346  }
347  }
348  }
349  };
350 
351 
352 
353  // MonomialBasisImpl
354  // -----------------
355 
356  template< class Topology, class F >
358 
359  template< class F >
360  class MonomialBasisImpl< GenericGeometry::Point, F >
361  {
363 
364  public:
365  typedef GenericGeometry::Point Topology;
366  typedef F Field;
367 
368  static const unsigned int dimDomain = Topology::dimension;
369 
370  typedef FieldVector< Field, dimDomain > DomainVector;
371 
372  private:
373  friend class MonomialBasis< Topology, Field >;
374  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
375  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
376 
377  template< int dimD >
378  void evaluate ( const unsigned int deriv, const unsigned int order,
379  const FieldVector< Field, dimD > &x,
380  const unsigned int block, const unsigned int *const offsets,
381  Field *const values ) const
382  {
383  *values = Unity< F >();
384  F *const end = values + block;
385  for( Field *it = values+1 ; it != end; ++it )
386  *it = Zero< F >();
387  }
388 
389  void integrate ( const unsigned int order,
390  const unsigned int *const offsets,
391  Field *const values ) const
392  {
393  values[ 0 ] = Unity< Field >();
394  }
395  };
396 
397  template< class BaseTopology, class F >
398  class MonomialBasisImpl< GenericGeometry::Prism< BaseTopology >, F >
399  {
401 
402  public:
403  typedef GenericGeometry::Prism< BaseTopology > Topology;
404  typedef F Field;
405 
406  static const unsigned int dimDomain = Topology::dimension;
407 
408  typedef FieldVector< Field, dimDomain > DomainVector;
409 
410  private:
411  friend class MonomialBasis< Topology, Field >;
412  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
413  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
414 
415  typedef MonomialBasisSize< BaseTopology > BaseSize;
416  typedef MonomialBasisSize< Topology > Size;
417 
418  MonomialBasisImpl< BaseTopology, Field > baseBasis_;
419 
421  {}
422 
423  template< int dimD >
424  void evaluate ( const unsigned int deriv, const unsigned int order,
425  const FieldVector< Field, dimD > &x,
426  const unsigned int block, const unsigned int *const offsets,
427  Field *const values ) const
428  {
430  const BaseSize &size = BaseSize::instance();
431 
432  const Field &z = x[ dimDomain-1 ];
433 
434  // fill first column
435  baseBasis_.evaluate( deriv, order, x, block, offsets, values );
436 
437  Field *row0 = values;
438  for( unsigned int k = 1; k <= order; ++k )
439  {
440  Field *row1 = values + block*offsets[ k-1 ];
441  Field *wit = row1 + block*size.sizes_[ k ];
442  Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
443  Helper::copy( deriv, wit, row0, size( k-1 ), z );
444  row0 = row1;
445  }
446  }
447 
448  void integrate ( const unsigned int order,
449  const unsigned int *const offsets,
450  Field *const values ) const
451  {
452  const BaseSize &size = BaseSize::instance();
453  const Size &mySize = Size::instance();
454  // fill first column
455  baseBasis_.integrate( order, offsets, values );
456  const unsigned int *const baseSizes = size.sizes_;
457 
458  Field *row0 = values;
459  for( unsigned int k = 1; k <= order; ++k )
460  {
461  Field *const row1begin = values + offsets[ k-1 ];
462  Field *const row1End = row1begin + mySize.sizes_[ k ];
463  assert( (unsigned int)(row1End - values) <= offsets[ k ] );
464 
465  Field *row1 = row1begin;
466  Field *it = row1begin + baseSizes[ k ];
467  for( unsigned int j = 1; j <= k; ++j )
468  {
469  Field *const end = it + baseSizes[ k ];
470  assert( (unsigned int)(end - values) <= offsets[ k ] );
471  for( ; it != end; ++row1, ++it )
472  *it = (Field( j ) / Field( j+1 )) * (*row1);
473  }
474  for( ; it != row1End; ++row0, ++it )
475  *it = (Field( k ) / Field( k+1 )) * (*row0);
476  row0 = row1;
477  }
478  }
479  };
480 
481  template< class BaseTopology, class F >
482  class MonomialBasisImpl< GenericGeometry::Pyramid< BaseTopology >, F >
483  {
485 
486  public:
487  typedef GenericGeometry::Pyramid< BaseTopology > Topology;
488  typedef F Field;
489 
490  static const unsigned int dimDomain = Topology::dimension;
491 
492  typedef FieldVector< Field, dimDomain > DomainVector;
493 
494  private:
495  friend class MonomialBasis< Topology, Field >;
496  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
497  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
498 
499  typedef MonomialBasisSize< BaseTopology > BaseSize;
500  typedef MonomialBasisSize< Topology > Size;
501 
502  MonomialBasisImpl< BaseTopology, Field > baseBasis_;
503 
505  {
506  }
507 
508  template< int dimD >
509  void evaluateSimplexBase ( const unsigned int deriv, const unsigned int order,
510  const FieldVector< Field, dimD > &x,
511  const unsigned int block, const unsigned int *const offsets,
512  Field *const values,
513  const BaseSize &size ) const
514  {
515  baseBasis_.evaluate( deriv, order, x, block, offsets, values );
516  }
517 
518  template< int dimD >
519  void evaluatePyramidBase ( const unsigned int deriv, const unsigned int order,
520  const FieldVector< Field, dimD > &x,
521  const unsigned int block, const unsigned int *const offsets,
522  Field *const values,
523  const BaseSize &size ) const
524  {
525  Field omz = Unity< Field >() - x[ dimDomain-1 ];
526 
527  if( Zero< Field >() < omz )
528  {
529  const Field invomz = Unity< Field >() / omz;
530  FieldVector< Field, dimD > y;
531  for( unsigned int i = 0; i < dimDomain-1; ++i )
532  y[ i ] = x[ i ] * invomz;
533 
534  // fill first column
535  baseBasis_.evaluate( deriv, order, y, block, offsets, values );
536 
537  Field omzk = omz;
538  for( unsigned int k = 1; k <= order; ++k )
539  {
540  Field *it = values + block*offsets[ k-1 ];
541  Field *const end = it + block*size.sizes_[ k ];
542  for( ; it != end; ++it )
543  *it *= omzk;
544  omzk *= omz;
545  }
546  }
547  else
548  {
549  assert( deriv==0 );
550  *values = Unity< Field >();
551  for( unsigned int k = 1; k <= order; ++k )
552  {
553  Field *it = values + block*offsets[ k-1 ];
554  Field *const end = it + block*size.sizes_[ k ];
555  for( ; it != end; ++it )
556  *it = Zero< Field >();
557  }
558  }
559  }
560 
561  template< int dimD >
562  void evaluate ( const unsigned int deriv, const unsigned int order,
563  const FieldVector< Field, dimD > &x,
564  const unsigned int block, const unsigned int *const offsets,
565  Field *const values ) const
566  {
567  typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
568  const BaseSize &size = BaseSize::instance();
569 
571  evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
572  else
573  evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
574 
575  Field *row0 = values;
576  for( unsigned int k = 1; k <= order; ++k )
577  {
578  Field *row1 = values + block*offsets[ k-1 ];
579  Field *wit = row1 + block*size.sizes_[ k ];
580  Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
581  row0 = row1;
582  }
583  }
584 
585  void integrate ( const unsigned int order,
586  const unsigned int *const offsets,
587  Field *const values ) const
588  {
589  const BaseSize &size = BaseSize::instance();
590 
591  // fill first column
592  baseBasis_.integrate( order, offsets, values );
593 
594  const unsigned int *const baseSizes = size.sizes_;
595 
596  Field *const col0End = values + baseSizes[ 0 ];
597  for( Field *it = values; it != col0End; ++it )
598  *it *= Field( 1 ) / Field( int(dimDomain) );
599  Field *row0 = values;
600 
601  for( unsigned int k = 1; k <= order; ++k )
602  {
603  const Field factor = (Field( 1 ) / Field( k + dimDomain ));
604 
605  Field *const row1 = values+offsets[ k-1 ];
606  Field *const col0End = row1 + baseSizes[ k ];
607  Field *it = row1;
608  for( ; it != col0End; ++it )
609  *it *= factor;
610  for( unsigned int i = 1; i <= k; ++i )
611  {
612  Field *const end = it + baseSizes[ k-i ];
613  assert( (unsigned int)(end - values) <= offsets[ k ] );
614  for( ; it != end; ++row0, ++it )
615  *it = (*row0) * (Field( i ) * factor);
616  }
617  row0 = row1;
618  }
619  }
620  };
621 
622 
623 
624  // MonomialBasis
625  // -------------
626 
627  template< class Topology, class F >
628  class MonomialBasis
629  : public MonomialBasisImpl< Topology, F >
630  {
631  typedef MonomialBasis< Topology, F > This;
632  typedef MonomialBasisImpl< Topology, F > Base;
633 
634  public:
635  static const unsigned int dimension = Base::dimDomain;
636  static const unsigned int dimRange = 1;
637 
638  typedef typename Base::Field Field;
639 
640  typedef typename Base::DomainVector DomainVector;
641 
642  typedef Dune::FieldVector<Field,dimRange> RangeVector;
643 
645 
646  MonomialBasis (unsigned int order)
647  : Base(),
648  order_(order),
649  size_(Size::instance())
650  {
651  assert(order<=1024); // avoid wrapping of unsigned int (0-1) order=1024 is quite hight...)
652  }
653 
654  const unsigned int *sizes ( unsigned int order ) const
655  {
656  size_.computeSizes( order );
657  return size_.numBaseFunctions_;
658  }
659 
660  const unsigned int *sizes () const
661  {
662  return sizes( order_ );
663  }
664 
665  const unsigned int size () const
666  {
667  size_.computeSizes( order_ );
668  return size_( order_ );
669  }
670 
671  const unsigned int derivSize ( const unsigned int deriv ) const
672  {
673  typedef typename GenericGeometry::SimplexTopology< dimension >::type SimplexTopology;
674  MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
676  }
677 
678  const unsigned int order () const
679  {
680  return order_ ;
681  }
682 
683  const unsigned int topologyId ( ) const
684  {
685  return Topology::id;
686  }
687 
688  void evaluate ( const unsigned int deriv, const DomainVector &x,
689  Field *const values ) const
690  {
691  Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
692  }
693 
694  template <unsigned int deriv>
695  void evaluate ( const DomainVector &x,
696  Field *const values ) const
697  {
698  evaluate( deriv, x, values );
699  }
700 
701  template<unsigned int deriv, class Vector >
702  void evaluate ( const DomainVector &x,
703  Vector &values ) const
704  {
705  evaluate<deriv>(x,&(values[0]));
706  }
707  template<unsigned int deriv, DerivativeLayout layout >
708  void evaluate ( const DomainVector &x,
710  {
711  evaluate<deriv>(x,&(values->block()));
712  }
713  template< unsigned int deriv >
714  void evaluate ( const DomainVector &x,
715  FieldVector<Field,Derivatives<Field,dimension,1,deriv,value>::size> *values ) const
716  {
717  evaluate(0,x,&(values[0][0]));
718  }
719 
720  template<class Vector >
721  void evaluate ( const DomainVector &x,
722  Vector &values ) const
723  {
724  evaluate<0>(x,&(values[0]));
725  }
726 
727  template< class DVector, class RVector >
728  void evaluate ( const DVector &x, RVector &values ) const
729  {
730  assert( DVector::dimension == dimension);
731  DomainVector bx;
732  for( int d = 0; d < dimension; ++d )
733  field_cast( x[ d ], bx[ d ] );
734  evaluate<0>( bx, values );
735  }
736 
737  void integrate ( Field *const values ) const
738  {
739  Base::integrate( order_, sizes( order_ ), values );
740  }
741  template <class Vector>
742  void integrate ( Vector &values ) const
743  {
744  integrate( &(values[ 0 ]) );
745  }
746  private:
747  MonomialBasis(const This&);
748  This& operator=(const This&);
749  unsigned int order_;
750  Size &size_;
751  };
752 
753 
754 
755  // StdMonomialBasis
756  // ----------------
757 
758  template< int dim,class F >
760  : public MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F >
761  {
764 
765  public:
766  typedef typename GenericGeometry::SimplexTopology< dim >::type Topology;
767  static const int dimension = dim;
768 
769  StandardMonomialBasis ( unsigned int order )
770  : Base( order )
771  {}
772  };
773 
774 
775 
776  // StandardBiMonomialBasis
777  // -----------------------
778 
779  template< int dim, class F >
781  : public MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F >
782  {
785 
786  public:
787  typedef typename GenericGeometry::CubeTopology< dim >::type Topology;
788  static const int dimension = dim;
789 
790  StandardBiMonomialBasis ( unsigned int order )
791  : Base( order )
792  {}
793  };
794 
795  // -----------------------------------------------------------
796  // -----------------------------------------------------------
797  // VirtualMonomialBasis
798  // -------------------
799 
800  template< int dim, class F >
802  {
804 
805  public:
806  typedef F Field;
807  typedef F StorageField;
808  static const int dimension = dim;
809  static const unsigned int dimRange = 1;
810 
811  typedef FieldVector<Field,dimension> DomainVector;
812  typedef FieldVector<Field,dimRange> RangeVector;
813 
814  explicit VirtualMonomialBasis(unsigned int topologyId,
815  unsigned int order)
816  : order_(order), topologyId_(topologyId) {}
817 
819 
820  virtual const unsigned int *sizes ( ) const = 0;
821 
822  const unsigned int size ( ) const
823  {
824  return sizes( )[ order_ ];
825  }
826 
827  const unsigned int order () const
828  {
829  return order_;
830  }
831 
832  const unsigned int topologyId ( ) const
833  {
834  return topologyId_;
835  }
836 
837  virtual void evaluate ( const unsigned int deriv, const DomainVector &x,
838  Field *const values ) const = 0;
839  template < unsigned int deriv >
840  void evaluate ( const DomainVector &x,
841  Field *const values ) const
842  {
843  evaluate( deriv, x, values );
844  }
845  template < unsigned int deriv, int size >
846  void evaluate ( const DomainVector &x,
847  Dune::FieldVector<Field,size> *const values ) const
848  {
849  evaluate( deriv, x, &(values[0][0]) );
850  }
851  template<unsigned int deriv, DerivativeLayout layout >
852  void evaluate ( const DomainVector &x,
854  {
855  evaluate<deriv>(x,&(values->block()));
856  }
857  template <unsigned int deriv, class Vector>
858  void evaluate ( const DomainVector &x,
859  Vector &values ) const
860  {
861  evaluate<deriv>( x, &(values[ 0 ]) );
862  }
863  template< class Vector >
864  void evaluate ( const DomainVector &x,
865  Vector &values ) const
866  {
867  evaluate<0>(x,values);
868  }
869  template< class DVector, class RVector >
870  void evaluate ( const DVector &x, RVector &values ) const
871  {
872  assert( DVector::dimension == dimension);
873  DomainVector bx;
874  for( int d = 0; d < dimension; ++d )
875  field_cast( x[ d ], bx[ d ] );
876  evaluate<0>( bx, values );
877  }
878  template< unsigned int deriv, class DVector, class RVector >
879  void evaluate ( const DVector &x, RVector &values ) const
880  {
881  assert( DVector::dimension == dimension);
882  DomainVector bx;
883  for( int d = 0; d < dimension; ++d )
884  field_cast( x[ d ], bx[ d ] );
885  evaluate<deriv>( bx, values );
886  }
887 
888  virtual void integrate ( Field *const values ) const = 0;
889  template <class Vector>
890  void integrate ( Vector &values ) const
891  {
892  integrate( &(values[ 0 ]) );
893  }
894  protected:
895  unsigned int order_;
896  unsigned int topologyId_;
897  };
898 
899  template< class Topology, class F >
901  : public VirtualMonomialBasis< Topology::dimension, F >
902  {
905 
906  public:
907  typedef typename Base::Field Field;
909 
910  VirtualMonomialBasisImpl(unsigned int order)
911  : Base(Topology::id,order), basis_(order)
912  {}
913 
914  const unsigned int *sizes ( ) const
915  {
916  return basis_.sizes(order_);
917  }
918 
919  void evaluate ( const unsigned int deriv, const DomainVector &x,
920  Field *const values ) const
921  {
922  basis_.evaluate(deriv,x,values);
923  }
924 
925  void integrate ( Field *const values ) const
926  {
927  basis_.integrate(values);
928  }
929 
930  private:
932  using Base::order_;
933  };
934 
935  // MonomialBasisFactory
936  // --------------------
937 
938  template< int dim, class F >
940 
941  template< int dim, class F >
943  {
944  static const unsigned int dimension = dim;
945  typedef unsigned int Key;
948  };
949 
950  template< int dim, class F >
951  struct MonomialBasisFactory :
952  public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
953  {
954  static const unsigned int dimension = dim;
955  typedef F StorageField;
957 
958  typedef typename Traits::Key Key;
959  typedef typename Traits::Object Object;
960 
961  template < int dd, class FF >
963  {
965  };
966 
967  template< class Topology >
968  static Object* createObject ( const Key &order )
969  {
971  }
972  };
973 
974 
975 
976  // MonomialBasisProvider
977  // ---------------------
978 
979  template< int dim, class SF >
981  : public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
982  {
983  static const unsigned int dimension = dim;
984  typedef SF StorageField;
985  template < int dd, class FF >
987  {
989  };
990  };
991 
992 }
993 
994 #endif