dune-localfunctions  2.4.1
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/genericgeometry/topologytypes.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 nomral
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 in 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< GenericGeometry::Point >
82  {
84 
85  public:
86  static This &instance ()
87  {
88  static This _instance;
89  return _instance;
90  }
91 
92  typedef GenericGeometry::Point Topology;
93 
94  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
95  friend class MonomialBasisSize< GenericGeometry::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< GenericGeometry::Prism< BaseTopology > >
151  {
153 
154  public:
155  static This &instance ()
156  {
157  static This _instance;
158  return _instance;
159  }
160 
161  typedef GenericGeometry::Prism< BaseTopology > Topology;
162 
163  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
164  friend class MonomialBasisSize< GenericGeometry::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< GenericGeometry::Pyramid< BaseTopology > >
226  {
228 
229  public:
230  static This &instance ()
231  {
232  static This _instance;
233  return _instance;
234  }
235 
236  typedef GenericGeometry::Pyramid< BaseTopology > Topology;
237 
238  friend class MonomialBasisSize< GenericGeometry::Prism< Topology > >;
239  friend class MonomialBasisSize< GenericGeometry::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  const F *const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
332  for( ; rit != drend ; ++rit, ++wit )
333  *wit = z * *rit;
334  for (unsigned int j=1; j<d; ++j)
335  {
336  const F *const drend = rit + mySize.sizes_[ d-j ] - mySize.sizes_[ d-j-1 ];
337  for( ; rit != drend ; ++prit, ++rit, ++wit )
338  *wit = F(j) * *prit + z * *rit;
339  }
340  *wit = F(d) * *prit + z * *rit;
341  ++prit, ++rit, ++wit;
342  assert(derivEnd == rit);
343  rit += size.sizes_[d] - mySize.sizes_[d];
344  prit += size.sizes_[d-1] - mySize.sizes_[d-1];
345  const F *const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
346  for ( ; wit != emptyWitEnd; ++wit )
347  *wit = Zero<F>();
348  }
349  }
350  }
351  };
352 
353 
354 
355  // MonomialBasisImpl
356  // -----------------
357 
358  template< class Topology, class F >
360 
361  template< class F >
362  class MonomialBasisImpl< GenericGeometry::Point, F >
363  {
365 
366  public:
367  typedef GenericGeometry::Point Topology;
368  typedef F Field;
369 
370  static const unsigned int dimDomain = Topology::dimension;
371 
372  typedef FieldVector< Field, dimDomain > DomainVector;
373 
374  private:
375  friend class MonomialBasis< Topology, Field >;
376  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
377  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
378 
379  template< int dimD >
380  void evaluate ( const unsigned int deriv, const unsigned int order,
381  const FieldVector< Field, dimD > &x,
382  const unsigned int block, const unsigned int *const offsets,
383  Field *const values ) const
384  {
385  *values = Unity< F >();
386  F *const end = values + block;
387  for( Field *it = values+1 ; it != end; ++it )
388  *it = Zero< F >();
389  }
390 
391  void integrate ( const unsigned int order,
392  const unsigned int *const offsets,
393  Field *const values ) const
394  {
395  values[ 0 ] = Unity< Field >();
396  }
397  };
398 
399  template< class BaseTopology, class F >
400  class MonomialBasisImpl< GenericGeometry::Prism< BaseTopology >, F >
401  {
403 
404  public:
405  typedef GenericGeometry::Prism< BaseTopology > Topology;
406  typedef F Field;
407 
408  static const unsigned int dimDomain = Topology::dimension;
409 
410  typedef FieldVector< Field, dimDomain > DomainVector;
411 
412  private:
413  friend class MonomialBasis< Topology, Field >;
414  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
415  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
416 
417  typedef MonomialBasisSize< BaseTopology > BaseSize;
418  typedef MonomialBasisSize< Topology > Size;
419 
420  MonomialBasisImpl< BaseTopology, Field > baseBasis_;
421 
423  {}
424 
425  template< int dimD >
426  void evaluate ( const unsigned int deriv, const unsigned int order,
427  const FieldVector< Field, dimD > &x,
428  const unsigned int block, const unsigned int *const offsets,
429  Field *const values ) const
430  {
432  const BaseSize &size = BaseSize::instance();
433 
434  const Field &z = x[ dimDomain-1 ];
435 
436  // fill first column
437  baseBasis_.evaluate( deriv, order, x, block, offsets, values );
438 
439  Field *row0 = values;
440  for( unsigned int k = 1; k <= order; ++k )
441  {
442  Field *row1 = values + block*offsets[ k-1 ];
443  Field *wit = row1 + block*size.sizes_[ k ];
444  Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
445  Helper::copy( deriv, wit, row0, size( k-1 ), z );
446  row0 = row1;
447  }
448  }
449 
450  void integrate ( const unsigned int order,
451  const unsigned int *const offsets,
452  Field *const values ) const
453  {
454  const BaseSize &size = BaseSize::instance();
455  const Size &mySize = Size::instance();
456  // fill first column
457  baseBasis_.integrate( order, offsets, values );
458  const unsigned int *const baseSizes = size.sizes_;
459 
460  Field *row0 = values;
461  for( unsigned int k = 1; k <= order; ++k )
462  {
463  Field *const row1begin = values + offsets[ k-1 ];
464  Field *const row1End = row1begin + mySize.sizes_[ k ];
465  assert( (unsigned int)(row1End - values) <= offsets[ k ] );
466 
467  Field *row1 = row1begin;
468  Field *it = row1begin + baseSizes[ k ];
469  for( unsigned int j = 1; j <= k; ++j )
470  {
471  Field *const end = it + baseSizes[ k ];
472  assert( (unsigned int)(end - values) <= offsets[ k ] );
473  for( ; it != end; ++row1, ++it )
474  *it = (Field( j ) / Field( j+1 )) * (*row1);
475  }
476  for( ; it != row1End; ++row0, ++it )
477  *it = (Field( k ) / Field( k+1 )) * (*row0);
478  row0 = row1;
479  }
480  }
481  };
482 
483  template< class BaseTopology, class F >
484  class MonomialBasisImpl< GenericGeometry::Pyramid< BaseTopology >, F >
485  {
487 
488  public:
489  typedef GenericGeometry::Pyramid< BaseTopology > Topology;
490  typedef F Field;
491 
492  static const unsigned int dimDomain = Topology::dimension;
493 
494  typedef FieldVector< Field, dimDomain > DomainVector;
495 
496  private:
497  friend class MonomialBasis< Topology, Field >;
498  friend class MonomialBasisImpl< GenericGeometry::Prism< Topology >, Field >;
499  friend class MonomialBasisImpl< GenericGeometry::Pyramid< Topology >, Field >;
500 
501  typedef MonomialBasisSize< BaseTopology > BaseSize;
502  typedef MonomialBasisSize< Topology > Size;
503 
504  MonomialBasisImpl< BaseTopology, Field > baseBasis_;
505 
507  {}
508 
509  template< int dimD >
510  void evaluateSimplexBase ( const unsigned int deriv, const unsigned int order,
511  const FieldVector< Field, dimD > &x,
512  const unsigned int block, const unsigned int *const offsets,
513  Field *const values,
514  const BaseSize &size ) const
515  {
516  baseBasis_.evaluate( deriv, order, x, block, offsets, values );
517  }
518 
519  template< int dimD >
520  void evaluatePyramidBase ( const unsigned int deriv, const unsigned int order,
521  const FieldVector< Field, dimD > &x,
522  const unsigned int block, const unsigned int *const offsets,
523  Field *const values,
524  const BaseSize &size ) const
525  {
526  Field omz = Unity< Field >() - x[ dimDomain-1 ];
527 
528  if( Zero< Field >() < omz )
529  {
530  const Field invomz = Unity< Field >() / omz;
531  FieldVector< Field, dimD > y;
532  for( unsigned int i = 0; i < dimDomain-1; ++i )
533  y[ i ] = x[ i ] * invomz;
534 
535  // fill first column
536  baseBasis_.evaluate( deriv, order, y, block, offsets, values );
537 
538  Field omzk = omz;
539  for( unsigned int k = 1; k <= order; ++k )
540  {
541  Field *it = values + block*offsets[ k-1 ];
542  Field *const end = it + block*size.sizes_[ k ];
543  for( ; it != end; ++it )
544  *it *= omzk;
545  omzk *= omz;
546  }
547  }
548  else
549  {
550  assert( deriv==0 );
551  *values = Unity< Field >();
552  for( unsigned int k = 1; k <= order; ++k )
553  {
554  Field *it = values + block*offsets[ k-1 ];
555  Field *const end = it + block*size.sizes_[ k ];
556  for( ; it != end; ++it )
557  *it = Zero< Field >();
558  }
559  }
560  }
561 
562  template< int dimD >
563  void evaluate ( const unsigned int deriv, const unsigned int order,
564  const FieldVector< Field, dimD > &x,
565  const unsigned int block, const unsigned int *const offsets,
566  Field *const values ) const
567  {
569  const BaseSize &size = BaseSize::instance();
570 
572  evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
573  else
574  evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
575 
576  Field *row0 = values;
577  for( unsigned int k = 1; k <= order; ++k )
578  {
579  Field *row1 = values + block*offsets[ k-1 ];
580  Field *wit = row1 + block*size.sizes_[ k ];
581  Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
582  row0 = row1;
583  }
584  }
585 
586  void integrate ( const unsigned int order,
587  const unsigned int *const offsets,
588  Field *const values ) const
589  {
590  const BaseSize &size = BaseSize::instance();
591 
592  // fill first column
593  baseBasis_.integrate( order, offsets, values );
594 
595  const unsigned int *const baseSizes = size.sizes_;
596 
597  Field *const col0End = values + baseSizes[ 0 ];
598  for( Field *it = values; it != col0End; ++it )
599  *it *= Field( 1 ) / Field( int(dimDomain) );
600  Field *row0 = values;
601 
602  for( unsigned int k = 1; k <= order; ++k )
603  {
604  const Field factor = (Field( 1 ) / Field( k + dimDomain ));
605 
606  Field *const row1 = values+offsets[ k-1 ];
607  Field *const col0End = row1 + baseSizes[ k ];
608  Field *it = row1;
609  for( ; it != col0End; ++it )
610  *it *= factor;
611  for( unsigned int i = 1; i <= k; ++i )
612  {
613  Field *const end = it + baseSizes[ k-i ];
614  assert( (unsigned int)(end - values) <= offsets[ k ] );
615  for( ; it != end; ++row0, ++it )
616  *it = (*row0) * (Field( i ) * factor);
617  }
618  row0 = row1;
619  }
620  }
621  };
622 
623 
624 
625  // MonomialBasis
626  // -------------
627 
628  template< class Topology, class F >
629  class MonomialBasis
630  : public MonomialBasisImpl< Topology, F >
631  {
634 
635  public:
636  static const unsigned int dimension = Base::dimDomain;
637  static const unsigned int dimRange = 1;
638 
639  typedef typename Base::Field Field;
640 
641  typedef typename Base::DomainVector DomainVector;
642 
643  typedef Dune::FieldVector<Field,dimRange> RangeVector;
644 
646 
647  MonomialBasis (unsigned int order)
648  : Base(),
649  order_(order),
650  size_(Size::instance())
651  {
652  assert(order<=1024); // avoid wrapping of unsigned int (0-1) order=1024 is quite hight...)
653  }
654 
655  const unsigned int *sizes ( unsigned int order ) const
656  {
657  size_.computeSizes( order );
658  return size_.numBaseFunctions_;
659  }
660 
661  const unsigned int *sizes () const
662  {
663  return sizes( order_ );
664  }
665 
666  const unsigned int size () const
667  {
668  size_.computeSizes( order_ );
669  return size_( order_ );
670  }
671 
672  const unsigned int derivSize ( const unsigned int deriv ) const
673  {
674  typedef typename GenericGeometry::SimplexTopology< dimension >::type SimplexTopology;
675  MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
677  }
678 
679  const unsigned int order () const
680  {
681  return order_ ;
682  }
683 
684  const unsigned int topologyId ( ) const
685  {
686  return Topology::id;
687  }
688 
689  void evaluate ( const unsigned int deriv, const DomainVector &x,
690  Field *const values ) const
691  {
692  Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
693  }
694 
695  template <unsigned int deriv>
696  void evaluate ( const DomainVector &x,
697  Field *const values ) const
698  {
699  evaluate( deriv, x, values );
700  }
701 
702  template<unsigned int deriv, class Vector >
703  void evaluate ( const DomainVector &x,
704  Vector &values ) const
705  {
706  evaluate<deriv>(x,&(values[0]));
707  }
708  template<unsigned int deriv, DerivativeLayout layout >
709  void evaluate ( const DomainVector &x,
711  {
712  evaluate<deriv>(x,&(values->block()));
713  }
714  template< unsigned int deriv >
715  void evaluate ( const DomainVector &x,
716  FieldVector<Field,Derivatives<Field,dimension,1,deriv,value>::size> *values ) const
717  {
718  evaluate(0,x,&(values[0][0]));
719  }
720 
721  template<class Vector >
722  void evaluate ( const DomainVector &x,
723  Vector &values ) const
724  {
725  evaluate<0>(x,&(values[0]));
726  }
727 
728  template< class DVector, class RVector >
729  void evaluate ( const DVector &x, RVector &values ) const
730  {
731  assert( DVector::dimension == dimension);
732  DomainVector bx;
733  for( int d = 0; d < dimension; ++d )
734  field_cast( x[ d ], bx[ d ] );
735  evaluate<0>( bx, values );
736  }
737 
738  void integrate ( Field *const values ) const
739  {
740  Base::integrate( order_, sizes( order_ ), values );
741  }
742  template <class Vector>
743  void integrate ( Vector &values ) const
744  {
745  integrate( &(values[ 0 ]) );
746  }
747  private:
748  MonomialBasis(const This&);
749  This& operator=(const This&);
750  unsigned int order_;
751  Size &size_;
752  };
753 
754 
755 
756  // StdMonomialBasis
757  // ----------------
758 
759  template< int dim,class F >
761  : public MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F >
762  {
765 
766  public:
767  typedef typename GenericGeometry::SimplexTopology< dim >::type Topology;
768  static const int dimension = dim;
769 
770  StandardMonomialBasis ( unsigned int order )
771  : Base( order )
772  {}
773  };
774 
775 
776 
777  // StandardBiMonomialBasis
778  // -----------------------
779 
780  template< int dim, class F >
782  : public MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F >
783  {
786 
787  public:
788  typedef typename GenericGeometry::CubeTopology< dim >::type Topology;
789  static const int dimension = dim;
790 
791  StandardBiMonomialBasis ( unsigned int order )
792  : Base( order )
793  {}
794  };
795 
796  // -----------------------------------------------------------
797  // -----------------------------------------------------------
798  // VirtualMonomialBasis
799  // -------------------
800 
801  template< int dim, class F >
803  {
805 
806  public:
807  typedef F Field;
808  typedef F StorageField;
809  static const int dimension = dim;
810  static const unsigned int dimRange = 1;
811 
812  typedef FieldVector<Field,dimension> DomainVector;
813  typedef FieldVector<Field,dimRange> RangeVector;
814 
815  explicit VirtualMonomialBasis(unsigned int topologyId,
816  unsigned int order)
817  : order_(order), topologyId_(topologyId) {}
818 
820 
821  virtual const unsigned int *sizes ( ) const = 0;
822 
823  const unsigned int size ( ) const
824  {
825  return sizes( )[ order_ ];
826  }
827 
828  const unsigned int order () const
829  {
830  return order_;
831  }
832 
833  const unsigned int topologyId ( ) const
834  {
835  return topologyId_;
836  }
837 
838  virtual void evaluate ( const unsigned int deriv, const DomainVector &x,
839  Field *const values ) const = 0;
840  template < unsigned int deriv >
841  void evaluate ( const DomainVector &x,
842  Field *const values ) const
843  {
844  evaluate( deriv, x, values );
845  }
846  template < unsigned int deriv, int size >
847  void evaluate ( const DomainVector &x,
848  Dune::FieldVector<Field,size> *const values ) const
849  {
850  evaluate( deriv, x, &(values[0][0]) );
851  }
852  template<unsigned int deriv, DerivativeLayout layout >
853  void evaluate ( const DomainVector &x,
855  {
856  evaluate<deriv>(x,&(values->block()));
857  }
858  template <unsigned int deriv, class Vector>
859  void evaluate ( const DomainVector &x,
860  Vector &values ) const
861  {
862  evaluate<deriv>( x, &(values[ 0 ]) );
863  }
864  template< class Vector >
865  void evaluate ( const DomainVector &x,
866  Vector &values ) const
867  {
868  evaluate<0>(x,values);
869  }
870  template< class DVector, class RVector >
871  void evaluate ( const DVector &x, RVector &values ) const
872  {
873  assert( DVector::dimension == dimension);
874  DomainVector bx;
875  for( int d = 0; d < dimension; ++d )
876  field_cast( x[ d ], bx[ d ] );
877  evaluate<0>( bx, values );
878  }
879  template< unsigned int deriv, class DVector, class RVector >
880  void evaluate ( const DVector &x, RVector &values ) const
881  {
882  assert( DVector::dimension == dimension);
883  DomainVector bx;
884  for( int d = 0; d < dimension; ++d )
885  field_cast( x[ d ], bx[ d ] );
886  evaluate<deriv>( bx, values );
887  }
888 
889  virtual void integrate ( Field *const values ) const = 0;
890  template <class Vector>
891  void integrate ( Vector &values ) const
892  {
893  integrate( &(values[ 0 ]) );
894  }
895  protected:
896  unsigned int order_;
897  unsigned int topologyId_;
898  };
899 
900  template< class Topology, class F >
902  : public VirtualMonomialBasis< Topology::dimension, F >
903  {
906 
907  public:
908  typedef typename Base::Field Field;
910 
911  VirtualMonomialBasisImpl(unsigned int order)
912  : Base(Topology::id,order), basis_(order)
913  {}
914 
915  const unsigned int *sizes ( ) const
916  {
917  return basis_.sizes(order_);
918  }
919 
920  void evaluate ( const unsigned int deriv, const DomainVector &x,
921  Field *const values ) const
922  {
923  basis_.evaluate(deriv,x,values);
924  }
925 
926  void integrate ( Field *const values ) const
927  {
928  basis_.integrate(values);
929  }
930 
931  private:
933  using Base::order_;
934  };
935 
936  // MonomialBasisFactory
937  // --------------------
938 
939  template< int dim, class F >
941 
942  template< int dim, class F >
944  {
945  static const unsigned int dimension = dim;
946  typedef unsigned int Key;
949  };
950 
951  template< int dim, class F >
952  struct MonomialBasisFactory :
953  public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
954  {
955  static const unsigned int dimension = dim;
956  typedef F StorageField;
958 
959  typedef typename Traits::Key Key;
960  typedef typename Traits::Object Object;
961 
962  template < int dd, class FF >
964  {
966  };
967 
968  template< class Topology >
969  static Object* createObject ( const Key &order )
970  {
972  }
973  };
974 
975 
976 
977  // MonomialBasisProvider
978  // ---------------------
979 
980  template< int dim, class SF >
982  : public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
983  {
984  static const unsigned int dimension = dim;
985  typedef SF StorageField;
986  template < int dd, class FF >
988  {
990  };
991  };
992 
993 }
994 
995 #endif
GenericGeometry::Prism< BaseTopology > Topology
Definition: monomialbasis.hh:405
FieldVector< Field, dimension > DomainVector
Definition: monomialbasis.hh:812
Definition: tensor.hh:165
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:703
void evaluate(const DomainVector &x, FieldVector< Field, Derivatives< Field, dimension, 1, deriv, value >::size > *values) const
Definition: monomialbasis.hh:715
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:853
Definition: monomialbasis.hh:940
static void copy(const unsigned int deriv, F *&wit, F *&rit, const unsigned int numBaseFunctions, const F &z)
Definition: monomialbasis.hh:311
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
SF StorageField
Definition: monomialbasis.hh:985
const unsigned int size() const
Definition: monomialbasis.hh:823
MonomialBasisSize< Topology > Size
Definition: monomialbasis.hh:645
Dune::FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:643
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:722
GenericGeometry::Point Topology
Definition: monomialbasis.hh:92
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:841
static This & instance()
Definition: monomialbasis.hh:86
Definition: monomialbasis.hh:901
GenericGeometry::Pyramid< BaseTopology > Topology
Definition: monomialbasis.hh:489
MonomialBasisFactory< dim, F > Factory
Definition: monomialbasis.hh:948
Definition: monomialbasis.hh:70
const unsigned int * sizes() const
Definition: monomialbasis.hh:915
Definition: monomialbasis.hh:781
virtual ~VirtualMonomialBasis()
Definition: monomialbasis.hh:819
void computeSizes(unsigned int order) const
Definition: monomialbasis.hh:271
GenericGeometry::CubeTopology< dim >::type Topology
Definition: monomialbasis.hh:788
GenericGeometry::Point Topology
Definition: monomialbasis.hh:367
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:865
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:871
MonomialBasisFactory< dd, FF > Type
Definition: monomialbasis.hh:965
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:709
Definition: monomialbasis.hh:359
void integrate(Field *const values) const
Definition: monomialbasis.hh:738
unsigned int maxOrder() const
Definition: monomialbasis.hh:122
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:696
StandardBiMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:791
VirtualMonomialBasisImpl(unsigned int order)
Definition: monomialbasis.hh:911
unsigned int Key
Definition: monomialbasis.hh:946
MonomialBasisSize< typename GenericGeometry::SimplexTopology< dim >::type > Size
Definition: monomialbasis.hh:309
void integrate(Vector &values) const
Definition: monomialbasis.hh:743
VirtualMonomialBasis(unsigned int topologyId, unsigned int order)
Definition: monomialbasis.hh:815
Definition: monomialbasis.hh:760
const unsigned int * sizes() const
Definition: monomialbasis.hh:661
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:920
Base::DomainVector DomainVector
Definition: monomialbasis.hh:909
Definition: monomialbasis.hh:306
A class representing the unit of a given Field.
Definition: field.hh:27
Definition: monomialbasis.hh:981
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:880
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:859
Definition: brezzidouglasmarini1cube2dlocalbasis.hh:14
F Field
Definition: monomialbasis.hh:807
GenericGeometry::Pyramid< BaseTopology > Topology
Definition: monomialbasis.hh:236
static This & instance()
Definition: monomialbasis.hh:230
const unsigned int topologyId() const
Definition: monomialbasis.hh:833
MonomialBasisSize< typename GenericGeometry::SimplexTopology< mydim >::type > MySize
Definition: monomialbasis.hh:308
const VirtualMonomialBasis< dimension, F > Object
Definition: monomialbasis.hh:947
MonomialBasisFactoryTraits< dim, F > Traits
Definition: monomialbasis.hh:957
const unsigned int order() const
Definition: monomialbasis.hh:679
FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:813
void integrate(Vector &values) const
Definition: monomialbasis.hh:891
const unsigned int size() const
Definition: monomialbasis.hh:666
Base::DomainVector DomainVector
Definition: monomialbasis.hh:641
const unsigned int derivSize(const unsigned int deriv) const
Definition: monomialbasis.hh:672
~MonomialBasisSize()
Definition: monomialbasis.hh:111
void evaluate(const DomainVector &x, Dune::FieldVector< Field, size > *const values) const
Definition: monomialbasis.hh:847
void computeSizes(unsigned int order) const
Definition: monomialbasis.hh:196
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:729
unsigned int order_
Definition: monomialbasis.hh:896
MonomialBasisProvider< dd, FF > Type
Definition: monomialbasis.hh:989
const unsigned int order() const
Definition: monomialbasis.hh:828
F StorageField
Definition: monomialbasis.hh:808
void integrate(Field *const values) const
Definition: monomialbasis.hh:926
static This & instance()
Definition: monomialbasis.hh:155
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:494
void computeSizes(unsigned int order) const
Definition: monomialbasis.hh:127
Traits::Object Object
Definition: monomialbasis.hh:960
Definition: tensor.hh:168
Base::Field Field
Definition: monomialbasis.hh:908
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:372
unsigned int topologyId_
Definition: monomialbasis.hh:897
Definition: monomialbasis.hh:943
GenericGeometry::Prism< BaseTopology > Topology
Definition: monomialbasis.hh:161
Definition: monomialbasis.hh:73
Traits::Key Key
Definition: monomialbasis.hh:959
F StorageField
Definition: monomialbasis.hh:956
static Object * createObject(const Key &order)
Definition: monomialbasis.hh:969
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:689
A class representing the zero of a given Field.
Definition: field.hh:76
unsigned int maxOrder() const
Definition: monomialbasis.hh:266
MonomialBasis(unsigned int order)
Definition: monomialbasis.hh:647
Definition: monomialbasis.hh:802
const unsigned int * sizes(unsigned int order) const
Definition: monomialbasis.hh:655
const unsigned int topologyId() const
Definition: monomialbasis.hh:684
GenericGeometry::SimplexTopology< dim >::type Topology
Definition: monomialbasis.hh:767
StandardMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:770
unsigned int maxOrder() const
Definition: monomialbasis.hh:191
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:410
Base::Field Field
Definition: monomialbasis.hh:639