3 #ifndef DUNE_MONOMIALBASIS_HH 4 #define DUNE_MONOMIALBASIS_HH 8 #include <dune/common/fvector.hh> 9 #include <dune/common/fmatrix.hh> 11 #include <dune/geometry/topologyfactory.hh> 12 #include <dune/geometry/genericgeometry/topologytypes.hh> 69 template<
class Topology >
72 template<
class Topology,
class F >
88 static This _instance;
97 mutable unsigned int maxOrder_;
99 mutable unsigned int *sizes_;
101 mutable unsigned int *numBaseFunctions_;
106 numBaseFunctions_( 0 )
114 delete[] numBaseFunctions_;
117 unsigned int operator() (
const unsigned int order )
const 119 return numBaseFunctions_[ order ];
129 if (order <= maxOrder_)
135 delete [] numBaseFunctions_;
136 sizes_ =
new unsigned int [ order+1 ];
137 numBaseFunctions_ =
new unsigned int [ order+1 ];
140 numBaseFunctions_[ 0 ] = 1;
141 for(
unsigned int k = 1; k <= order; ++k )
144 numBaseFunctions_[ k ] = 1;
149 template<
class BaseTopology >
157 static This _instance;
161 typedef GenericGeometry::Prism< BaseTopology >
Topology;
166 mutable unsigned int maxOrder_;
168 mutable unsigned int *sizes_;
170 mutable unsigned int *numBaseFunctions_;
175 numBaseFunctions_( 0 )
183 delete[] numBaseFunctions_;
186 unsigned int operator() (
const unsigned int order )
const 188 return numBaseFunctions_[ order ];
198 if (order <= maxOrder_)
204 delete[] numBaseFunctions_;
205 sizes_ =
new unsigned int[ order+1 ];
206 numBaseFunctions_ =
new unsigned int[ order+1 ];
210 baseBasis.computeSizes( order );
211 const unsigned int *
const baseSizes = baseBasis.sizes_;
212 const unsigned int *
const baseNBF = baseBasis.numBaseFunctions_;
215 numBaseFunctions_[ 0 ] = 1;
216 for(
unsigned int k = 1; k <= order; ++k )
218 sizes_[ k ] = baseNBF[ k ] + k*baseSizes[ k ];
219 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
224 template<
class BaseTopology >
232 static This _instance;
236 typedef GenericGeometry::Pyramid< BaseTopology >
Topology;
241 mutable unsigned int maxOrder_;
243 mutable unsigned int *sizes_;
245 mutable unsigned int *numBaseFunctions_;
250 numBaseFunctions_( 0 )
258 delete[] numBaseFunctions_;
261 unsigned int operator() (
const unsigned int order )
const 263 return numBaseFunctions_[ order ];
273 if (order <= maxOrder_)
279 delete[] numBaseFunctions_;
280 sizes_ =
new unsigned int[ order+1 ];
281 numBaseFunctions_ =
new unsigned int[ order+1 ];
286 baseBasis.computeSizes( order );
288 const unsigned int *
const baseNBF = baseBasis.numBaseFunctions_;
290 numBaseFunctions_[ 0 ] = 1;
291 for(
unsigned int k = 1; k <= order; ++k )
293 sizes_[ k ] = baseNBF[ k ];
294 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
305 template<
int mydim,
int dim,
class F >
311 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
312 const unsigned int numBaseFunctions,
const F &z )
315 MySize &mySize = MySize::instance();
316 Size &size = Size::instance();
318 const F *
const rend = rit + size( deriv )*numBaseFunctions;
319 for( ; rit != rend; )
326 for(
unsigned d = 1; d <= deriv; ++d )
329 const F *
const derivEnd = rit + mySize.sizes_[ d ];
331 const F *
const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
332 for( ; rit != drend ; ++rit, ++wit )
334 for (
unsigned int j=1; j<d; ++j)
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;
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 )
358 template<
class Topology,
class F >
370 static const unsigned int dimDomain = Topology::dimension;
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
386 F *
const end = values + block;
387 for( Field *it = values+1 ; it != end; ++it )
391 void integrate (
const unsigned int order,
392 const unsigned int *
const offsets,
393 Field *
const values )
const 399 template<
class BaseTopology,
class F >
405 typedef GenericGeometry::Prism< BaseTopology >
Topology;
408 static const unsigned int dimDomain = Topology::dimension;
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 432 const BaseSize &size = BaseSize::instance();
434 const Field &z = x[ dimDomain-1 ];
437 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
439 Field *row0 = values;
440 for(
unsigned int k = 1; k <= order; ++k )
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 );
450 void integrate (
const unsigned int order,
451 const unsigned int *
const offsets,
452 Field *
const values )
const 454 const BaseSize &size = BaseSize::instance();
455 const Size &mySize = Size::instance();
457 baseBasis_.integrate( order, offsets, values );
458 const unsigned int *
const baseSizes = size.sizes_;
460 Field *row0 = values;
461 for(
unsigned int k = 1; k <= order; ++k )
463 Field *
const row1begin = values + offsets[ k-1 ];
464 Field *
const row1End = row1begin + mySize.sizes_[ k ];
465 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
467 Field *row1 = row1begin;
468 Field *it = row1begin + baseSizes[ k ];
469 for(
unsigned int j = 1; j <= k; ++j )
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);
476 for( ; it != row1End; ++row0, ++it )
477 *it = (Field( k ) / Field( k+1 )) * (*row0);
483 template<
class BaseTopology,
class F >
489 typedef GenericGeometry::Pyramid< BaseTopology >
Topology;
492 static const unsigned int dimDomain = Topology::dimension;
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,
516 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
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,
531 FieldVector< Field, dimD > y;
532 for(
unsigned int i = 0; i < dimDomain-1; ++i )
533 y[ i ] = x[ i ] * invomz;
536 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
539 for(
unsigned int k = 1; k <= order; ++k )
541 Field *it = values + block*offsets[ k-1 ];
542 Field *
const end = it + block*size.sizes_[ k ];
543 for( ; it != end; ++it )
552 for(
unsigned int k = 1; k <= order; ++k )
554 Field *it = values + block*offsets[ k-1 ];
555 Field *
const end = it + block*size.sizes_[ k ];
556 for( ; it != end; ++it )
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 569 const BaseSize &size = BaseSize::instance();
572 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
574 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
576 Field *row0 = values;
577 for(
unsigned int k = 1; k <= order; ++k )
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 ] );
586 void integrate (
const unsigned int order,
587 const unsigned int *
const offsets,
588 Field *
const values )
const 590 const BaseSize &size = BaseSize::instance();
593 baseBasis_.integrate( order, offsets, values );
595 const unsigned int *
const baseSizes = size.sizes_;
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;
602 for(
unsigned int k = 1; k <= order; ++k )
604 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
606 Field *
const row1 = values+offsets[ k-1 ];
607 Field *
const col0End = row1 + baseSizes[ k ];
609 for( ; it != col0End; ++it )
611 for(
unsigned int i = 1; i <= k; ++i )
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);
628 template<
class Topology,
class F >
636 static const unsigned int dimension = Base::dimDomain;
637 static const unsigned int dimRange = 1;
650 size_(Size::instance())
655 const unsigned int *
sizes (
unsigned int order )
const 657 size_.computeSizes( order );
658 return size_.numBaseFunctions_;
663 return sizes( order_ );
666 const unsigned int size ()
const 668 size_.computeSizes( order_ );
669 return size_( order_ );
672 const unsigned int derivSize (
const unsigned int deriv )
const 674 typedef typename GenericGeometry::SimplexTopology< dimension >::type SimplexTopology;
689 void evaluate (
const unsigned int deriv,
const DomainVector &x,
690 Field *
const values )
const 692 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
695 template <
unsigned int deriv>
697 Field *
const values )
const 699 evaluate( deriv, x, values );
702 template<
unsigned int deriv,
class Vector >
704 Vector &values )
const 706 evaluate<deriv>(x,&(values[0]));
708 template<
unsigned int deriv, DerivativeLayout layout >
712 evaluate<deriv>(x,&(values->block()));
714 template<
unsigned int deriv >
718 evaluate(0,x,&(values[0][0]));
721 template<
class Vector >
723 Vector &values )
const 725 evaluate<0>(x,&(values[0]));
728 template<
class DVector,
class RVector >
729 void evaluate (
const DVector &x, RVector &values )
const 731 assert( DVector::dimension == dimension);
733 for(
int d = 0; d < dimension; ++d )
735 evaluate<0>( bx, values );
740 Base::integrate( order_, sizes( order_ ), values );
742 template <
class Vector>
745 integrate( &(values[ 0 ]) );
749 This& operator=(
const This&);
759 template<
int dim,
class F >
761 :
public MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F >
767 typedef typename GenericGeometry::SimplexTopology< dim >::type
Topology;
768 static const int dimension = dim;
780 template<
int dim,
class F >
782 :
public MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F >
788 typedef typename GenericGeometry::CubeTopology< dim >::type
Topology;
789 static const int dimension = dim;
801 template<
int dim,
class F >
809 static const int dimension = dim;
810 static const unsigned int dimRange = 1;
817 : order_(order), topologyId_(topologyId) {}
821 virtual const unsigned int *sizes ( )
const = 0;
823 const unsigned int size ( )
const 825 return sizes( )[ order_ ];
838 virtual void evaluate (
const unsigned int deriv,
const DomainVector &x,
839 Field *
const values )
const = 0;
840 template <
unsigned int deriv >
842 Field *
const values )
const 844 evaluate( deriv, x, values );
846 template <
unsigned int deriv,
int size >
848 Dune::FieldVector<Field,size> *
const values )
const 850 evaluate( deriv, x, &(values[0][0]) );
852 template<
unsigned int deriv, DerivativeLayout layout >
856 evaluate<deriv>(x,&(values->block()));
858 template <
unsigned int deriv,
class Vector>
860 Vector &values )
const 862 evaluate<deriv>( x, &(values[ 0 ]) );
864 template<
class Vector >
866 Vector &values )
const 868 evaluate<0>(x,values);
870 template<
class DVector,
class RVector >
871 void evaluate (
const DVector &x, RVector &values )
const 873 assert( DVector::dimension == dimension);
875 for(
int d = 0; d < dimension; ++d )
877 evaluate<0>( bx, values );
879 template<
unsigned int deriv,
class DVector,
class RVector >
880 void evaluate (
const DVector &x, RVector &values )
const 882 assert( DVector::dimension == dimension);
884 for(
int d = 0; d < dimension; ++d )
886 evaluate<deriv>( bx, values );
889 virtual void integrate ( Field *
const values )
const = 0;
890 template <
class Vector>
893 integrate( &(values[ 0 ]) );
900 template<
class Topology,
class F >
912 : Base(Topology::id,order), basis_(order)
917 return basis_.sizes(order_);
920 void evaluate (
const unsigned int deriv,
const DomainVector &x,
921 Field *
const values )
const 923 basis_.evaluate(deriv,x,values);
928 basis_.integrate(values);
939 template<
int dim,
class F >
942 template<
int dim,
class F >
945 static const unsigned int dimension = dim;
951 template<
int dim,
class F >
953 public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
955 static const unsigned int dimension = dim;
962 template <
int dd,
class FF >
968 template<
class Topology >
980 template<
int dim,
class SF >
982 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
984 static const unsigned int dimension = dim;
986 template <
int dd,
class FF >
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
Definition: monomialbasis.hh:150
F Field
Definition: monomialbasis.hh:490
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
~MonomialBasisSize()
Definition: monomialbasis.hh:180
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
Definition: monomialbasis.hh:484
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
~MonomialBasisSize()
Definition: monomialbasis.hh:255
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
F Field
Definition: monomialbasis.hh:406
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
Definition: monomialbasis.hh:987
unsigned int order_
Definition: monomialbasis.hh:896
MonomialBasisProvider< dd, FF > Type
Definition: monomialbasis.hh:989
Definition: monomialbasis.hh:400
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
F Field
Definition: monomialbasis.hh:368
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
Definition: monomialbasis.hh:225
Definition: monomialbasis.hh:963
Definition: monomialbasis.hh:362
MonomialBasis(unsigned int order)
Definition: monomialbasis.hh:647
Definition: monomialbasis.hh:802
Definition: monomialbasis.hh:81
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