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