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/type.hh> 69 template<
class Topology >
72 template<
class Topology,
class F >
86 static This &instance ()
88 static This _instance;
92 typedef Impl::Point Topology;
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 ];
122 unsigned int maxOrder ()
const 127 void computeSizes (
unsigned int order )
const 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 >
155 static This &instance ()
157 static This _instance;
161 typedef Impl::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 ];
191 unsigned int maxOrder()
const 196 void computeSizes (
unsigned int order )
const 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 >
230 static This &instance ()
232 static This _instance;
236 typedef Impl::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 ];
266 unsigned int maxOrder()
const 271 void computeSizes (
unsigned int order )
const 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 ];
333 const F *
const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
334 for( ; rit != drend ; ++rit, ++wit )
338 for (
unsigned int j=1; j<d; ++j)
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;
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 )
362 template<
class Topology,
class F >
374 static const unsigned int dimDomain = Topology::dimension;
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
390 F *
const end = values + block;
391 for( Field *it = values+1 ; it != end; ++it )
395 void integrate (
const unsigned int order,
396 const unsigned int *
const offsets,
397 Field *
const values )
const 403 template<
class BaseTopology,
class F >
412 static const unsigned int dimDomain = Topology::dimension;
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 436 const BaseSize &size = BaseSize::instance();
438 const Field &z = x[ dimDomain-1 ];
441 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
443 Field *row0 = values;
444 for(
unsigned int k = 1; k <= order; ++k )
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 );
454 void integrate (
const unsigned int order,
455 const unsigned int *
const offsets,
456 Field *
const values )
const 458 const BaseSize &size = BaseSize::instance();
459 const Size &mySize = Size::instance();
461 baseBasis_.integrate( order, offsets, values );
462 const unsigned int *
const baseSizes = size.sizes_;
464 Field *row0 = values;
465 for(
unsigned int k = 1; k <= order; ++k )
467 Field *
const row1begin = values + offsets[ k-1 ];
468 Field *
const row1End = row1begin + mySize.sizes_[ k ];
469 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
471 Field *row1 = row1begin;
472 Field *it = row1begin + baseSizes[ k ];
473 for(
unsigned int j = 1; j <= k; ++j )
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);
480 for( ; it != row1End; ++row0, ++it )
481 *it = (Field( k ) / Field( k+1 )) * (*row0);
487 template<
class BaseTopology,
class F >
496 static const unsigned int dimDomain = Topology::dimension;
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,
520 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
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,
535 FieldVector< Field, dimD > y;
536 for(
unsigned int i = 0; i < dimDomain-1; ++i )
537 y[ i ] = x[ i ] * invomz;
540 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
543 for(
unsigned int k = 1; k <= order; ++k )
545 Field *it = values + block*offsets[ k-1 ];
546 Field *
const end = it + block*size.sizes_[ k ];
547 for( ; it != end; ++it )
556 for(
unsigned int k = 1; k <= order; ++k )
558 Field *it = values + block*offsets[ k-1 ];
559 Field *
const end = it + block*size.sizes_[ k ];
560 for( ; it != end; ++it )
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 573 const BaseSize &size = BaseSize::instance();
576 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
578 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
580 Field *row0 = values;
581 for(
unsigned int k = 1; k <= order; ++k )
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 ] );
590 void integrate (
const unsigned int order,
591 const unsigned int *
const offsets,
592 Field *
const values )
const 594 const BaseSize &size = BaseSize::instance();
597 baseBasis_.integrate( order, offsets, values );
599 const unsigned int *
const baseSizes = size.sizes_;
602 Field *
const col0End = values + baseSizes[ 0 ];
603 for( Field *it = values; it != col0End; ++it )
604 *it *= Field( 1 ) / Field(
int(dimDomain) );
607 Field *row0 = values;
608 for(
unsigned int k = 1; k <= order; ++k )
610 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
612 Field *
const row1 = values+offsets[ k-1 ];
613 Field *
const col0End = row1 + baseSizes[ k ];
615 for( ; it != col0End; ++it )
617 for(
unsigned int i = 1; i <= k; ++i )
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);
634 template<
class Topology,
class F >
642 static const unsigned int dimension = Base::dimDomain;
643 static const unsigned int dimRange = 1;
656 size_(Size::instance())
661 const unsigned int *
sizes (
unsigned int order )
const 663 size_.computeSizes( order );
664 return size_.numBaseFunctions_;
669 return sizes( order_ );
674 size_.computeSizes( order_ );
675 return size_( order_ );
678 unsigned int derivSize (
const unsigned int deriv )
const 680 typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology;
695 void evaluate (
const unsigned int deriv,
const DomainVector &x,
696 Field *
const values )
const 698 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
701 template <
unsigned int deriv>
703 Field *
const values )
const 705 evaluate( deriv, x, values );
708 template<
unsigned int deriv,
class Vector >
710 Vector &values )
const 712 evaluate<deriv>(x,&(values[0]));
714 template<
unsigned int deriv, DerivativeLayout layout >
718 evaluate<deriv>(x,&(values->block()));
720 template<
unsigned int deriv >
724 evaluate(0,x,&(values[0][0]));
727 template<
class Vector >
729 Vector &values )
const 731 evaluate<0>(x,&(values[0]));
734 template<
class DVector,
class RVector >
735 void evaluate (
const DVector &x, RVector &values )
const 737 assert( DVector::dimension == dimension);
739 for(
int d = 0; d < dimension; ++d )
741 evaluate<0>( bx, values );
746 Base::integrate( order_, sizes( order_ ), values );
748 template <
class Vector>
751 integrate( &(values[ 0 ]) );
755 This& operator=(
const This&);
765 template<
int dim,
class F >
767 :
public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F >
773 typedef typename Impl::SimplexTopology< dim >::type
Topology;
774 static const int dimension = dim;
786 template<
int dim,
class F >
788 :
public MonomialBasis< typename Impl::CubeTopology< dim >::type, F >
794 typedef typename Impl::CubeTopology< dim >::type
Topology;
795 static const int dimension = dim;
807 template<
int dim,
class F >
815 static const int dimension = dim;
816 static const unsigned int dimRange = 1;
823 : order_(order), topologyId_(topologyId) {}
827 virtual const unsigned int *sizes ( )
const = 0;
831 return sizes( )[ order_ ];
844 virtual void evaluate (
const unsigned int deriv,
const DomainVector &x,
845 Field *
const values )
const = 0;
846 template <
unsigned int deriv >
848 Field *
const values )
const 850 evaluate( deriv, x, values );
852 template <
unsigned int deriv,
int size >
854 Dune::FieldVector<Field,size> *
const values )
const 856 evaluate( deriv, x, &(values[0][0]) );
858 template<
unsigned int deriv, DerivativeLayout layout >
862 evaluate<deriv>(x,&(values->block()));
864 template <
unsigned int deriv,
class Vector>
866 Vector &values )
const 868 evaluate<deriv>( x, &(values[ 0 ]) );
870 template<
class Vector >
872 Vector &values )
const 874 evaluate<0>(x,values);
876 template<
class DVector,
class RVector >
877 void evaluate (
const DVector &x, RVector &values )
const 879 assert( DVector::dimension == dimension);
881 for(
int d = 0; d < dimension; ++d )
883 evaluate<0>( bx, values );
885 template<
unsigned int deriv,
class DVector,
class RVector >
886 void evaluate (
const DVector &x, RVector &values )
const 888 assert( DVector::dimension == dimension);
890 for(
int d = 0; d < dimension; ++d )
892 evaluate<deriv>( bx, values );
895 virtual void integrate ( Field *
const values )
const = 0;
896 template <
class Vector>
899 integrate( &(values[ 0 ]) );
906 template<
class Topology,
class F >
918 : Base(Topology::id,order), basis_(order)
923 return basis_.sizes(order_);
926 void evaluate (
const unsigned int deriv,
const DomainVector &x,
927 Field *
const values )
const 929 basis_.evaluate(deriv,x,values);
934 basis_.integrate(values);
945 template<
int dim,
class F >
948 template<
int dim,
class F >
951 static const unsigned int dimension = dim;
957 template<
int dim,
class F >
959 public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
961 static const unsigned int dimension = dim;
968 template <
int dd,
class FF >
974 template<
class Topology >
986 template<
int dim,
class SF >
988 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
990 static const unsigned int dimension = dim;
992 template <
int dd,
class FF >
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
Definition: monomialbasis.hh:993
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
Definition: monomialbasis.hh:969
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
F Field
Definition: monomialbasis.hh:494
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
F Field
Definition: monomialbasis.hh:410
Impl::SimplexTopology< dim >::type Topology
Definition: monomialbasis.hh:773