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 >
70 class MonomialBasisSize;
72 template<
class Topology,
class F >
80 template<
class TopologyType >
88 static This _instance;
133 sizes_ =
new unsigned int[ order+1 ];
136 constexpr
auto dim = TopologyType::dimension;
139 for(
unsigned int k = 1; k <= order; ++k )
144 for(
int codim=dim-1; codim>=0; codim--)
146 if (Impl::isPrism(TopologyType::id,dim,codim))
148 for(
unsigned int k = 1; k <= order; ++k )
156 for(
unsigned int k = 1; k <= order; ++k )
172 template<
int mydim,
int dim,
class F >
178 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
179 const unsigned int numBaseFunctions,
const F &z )
185 const F *
const rend = rit + size( deriv )*numBaseFunctions;
186 for( ; rit != rend; )
193 for(
unsigned d = 1; d <= deriv; ++d )
196 const F *
const derivEnd = rit + mySize.
sizes_[ d ];
200 const F *
const drend = rit + mySize.
sizes_[ d ] - mySize.
sizes_[ d-1 ];
201 for( ; rit != drend ; ++rit, ++wit )
205 for (
unsigned int j=1; j<d; ++j)
207 const F *
const drend = rit + mySize.
sizes_[ d-j ] - mySize.
sizes_[ d-j-1 ];
208 for( ; rit != drend ; ++prit, ++rit, ++wit )
209 *wit = F(j) * *prit + z * *rit;
211 *wit = F(d) * *prit + z * *rit;
212 ++prit, ++rit, ++wit;
213 assert(derivEnd == rit);
216 const F *
const emptyWitEnd = wit + size.
sizes_[d] - mySize.
sizes_[d];
217 for ( ; wit != emptyWitEnd; ++wit )
229 template<
class Topology,
class F >
241 static const unsigned int dimDomain = Topology::dimension;
251 void evaluate ( const unsigned int deriv, const unsigned int order,
252 const FieldVector< Field, dimD > &x,
253 const unsigned int block, const unsigned int *const offsets,
254 Field *const values ) const
257 F *
const end = values + block;
258 for(
Field *it = values+1 ; it != end; ++it )
262 void integrate (
const unsigned int order,
263 const unsigned int *
const offsets,
264 Field *
const values )
const
270 template<
class BaseTopology,
class F >
279 static const unsigned int dimDomain = Topology::dimension;
297 void evaluate (
const unsigned int deriv,
const unsigned int order,
298 const FieldVector< Field, dimD > &x,
299 const unsigned int block,
const unsigned int *
const offsets,
300 Field *
const values )
const
303 const BaseSize &size = BaseSize::instance();
304 const_cast<BaseSize&
>(size).computeSizes(order);
306 const Field &z = x[ dimDomain-1 ];
309 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
311 Field *row0 = values;
312 for(
unsigned int k = 1; k <= order; ++k )
314 Field *row1 = values + block*offsets[ k-1 ];
315 Field *wit = row1 + block*size.sizes_[ k ];
316 Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
317 Helper::copy( deriv, wit, row0, size( k-1 ), z );
322 void integrate (
const unsigned int order,
323 const unsigned int *
const offsets,
324 Field *
const values )
const
326 const BaseSize &size = BaseSize::instance();
327 const Size &mySize = Size::instance();
329 baseBasis_.integrate( order, offsets, values );
330 const unsigned int *
const baseSizes = size.sizes_;
332 Field *row0 = values;
333 for(
unsigned int k = 1; k <= order; ++k )
335 Field *
const row1begin = values + offsets[ k-1 ];
336 Field *
const row1End = row1begin + mySize.sizes_[ k ];
337 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
339 Field *row1 = row1begin;
340 Field *it = row1begin + baseSizes[ k ];
341 for(
unsigned int j = 1; j <= k; ++j )
343 Field *
const end = it + baseSizes[ k ];
344 assert( (
unsigned int)(end - values) <= offsets[ k ] );
345 for( ; it != end; ++row1, ++it )
346 *it = (Field( j ) / Field( j+1 )) * (*row1);
348 for( ; it != row1End; ++row0, ++it )
349 *it = (Field( k ) / Field( k+1 )) * (*row0);
355 template<
class BaseTopology,
class F >
364 static const unsigned int dimDomain = Topology::dimension;
382 void evaluateSimplexBase (
const unsigned int deriv,
const unsigned int order,
383 const FieldVector< Field, dimD > &x,
384 const unsigned int block,
const unsigned int *
const offsets,
386 const BaseSize &size )
const
388 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
392 void evaluatePyramidBase (
const unsigned int deriv,
const unsigned int order,
393 const FieldVector< Field, dimD > &x,
394 const unsigned int block,
const unsigned int *
const offsets,
396 const BaseSize &size )
const
398 Field omz = Unity< Field >() - x[ dimDomain-1 ];
400 if( Zero< Field >() < omz )
402 const Field invomz = Unity< Field >() / omz;
403 FieldVector< Field, dimD > y;
404 for(
unsigned int i = 0; i < dimDomain-1; ++i )
405 y[ i ] = x[ i ] * invomz;
408 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
411 for(
unsigned int k = 1; k <= order; ++k )
413 Field *it = values + block*offsets[ k-1 ];
414 Field *
const end = it + block*size.sizes_[ k ];
415 for( ; it != end; ++it )
423 *values = Unity< Field >();
424 for(
unsigned int k = 1; k <= order; ++k )
426 Field *it = values + block*offsets[ k-1 ];
427 Field *
const end = it + block*size.sizes_[ k ];
428 for( ; it != end; ++it )
429 *it = Zero< Field >();
435 void evaluate (
const unsigned int deriv,
const unsigned int order,
436 const FieldVector< Field, dimD > &x,
437 const unsigned int block,
const unsigned int *
const offsets,
438 Field *
const values )
const
440 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
441 const BaseSize &size = BaseSize::instance();
442 const_cast<BaseSize&
>(size).computeSizes(order);
445 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
447 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
449 Field *row0 = values;
450 for(
unsigned int k = 1; k <= order; ++k )
452 Field *row1 = values + block*offsets[ k-1 ];
453 Field *wit = row1 + block*size.sizes_[ k ];
454 Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
459 void integrate (
const unsigned int order,
460 const unsigned int *
const offsets,
461 Field *
const values )
const
463 const BaseSize &size = BaseSize::instance();
466 baseBasis_.integrate( order, offsets, values );
468 const unsigned int *
const baseSizes = size.sizes_;
471 Field *
const col0End = values + baseSizes[ 0 ];
472 for( Field *it = values; it != col0End; ++it )
473 *it *= Field( 1 ) / Field(
int(dimDomain) );
476 Field *row0 = values;
477 for(
unsigned int k = 1; k <= order; ++k )
479 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
481 Field *
const row1 = values+offsets[ k-1 ];
482 Field *
const col0End = row1 + baseSizes[ k ];
484 for( ; it != col0End; ++it )
486 for(
unsigned int i = 1; i <= k; ++i )
488 Field *
const end = it + baseSizes[ k-i ];
489 assert( (
unsigned int)(end - values) <= offsets[ k ] );
490 for( ; it != end; ++row0, ++it )
491 *it = (*row0) * (Field( i ) * factor);
503 template<
class Topology,
class F >
525 size_(
Size::instance())
538 return sizes( order_ );
544 return size_( order_ );
547 unsigned int derivSize (
const unsigned int deriv )
const
549 typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology;
565 Field *
const values )
const
567 Base::evaluate( deriv, order_, x,
derivSize( deriv ),
sizes( order_ ), values );
570 template <
unsigned int deriv>
572 Field *
const values )
const
577 template<
unsigned int deriv,
class Vector >
579 Vector &values )
const
581 evaluate<deriv>(x,&(values[0]));
583 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
587 evaluate<deriv>(x,&(values->block()));
589 template<
unsigned int deriv >
596 template<
class Vector >
598 Vector &values )
const
600 evaluate<0>(x,&(values[0]));
603 template<
class DVector,
class RVector >
604 void evaluate (
const DVector &x, RVector &values )
const
606 assert( DVector::dimension ==
dimension);
610 evaluate<0>( bx, values );
615 Base::integrate( order_,
sizes( order_ ), values );
617 template <
class Vector>
624 This& operator=(
const This&);
634 template<
int dim,
class F >
636 :
public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F >
642 typedef typename Impl::SimplexTopology< dim >::type
Topology;
655 template<
int dim,
class F >
657 :
public MonomialBasis< typename Impl::CubeTopology< dim >::type, F >
663 typedef typename Impl::CubeTopology< dim >::type
Topology;
676 template<
int dim,
class F >
696 virtual const unsigned int *
sizes ( )
const = 0;
714 Field *
const values )
const = 0;
715 template <
unsigned int deriv >
717 Field *
const values )
const
721 template <
unsigned int deriv,
int size >
723 Dune::FieldVector<Field,size> *
const values )
const
725 evaluate( deriv, x, &(values[0][0]) );
727 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
731 evaluate<deriv>(x,&(values->block()));
733 template <
unsigned int deriv,
class Vector>
735 Vector &values )
const
737 evaluate<deriv>( x, &(values[ 0 ]) );
739 template<
class Vector >
741 Vector &values )
const
743 evaluate<0>(x,values);
745 template<
class DVector,
class RVector >
746 void evaluate (
const DVector &x, RVector &values )
const
748 assert( DVector::dimension ==
dimension);
752 evaluate<0>( bx, values );
754 template<
unsigned int deriv,
class DVector,
class RVector >
755 void evaluate (
const DVector &x, RVector &values )
const
757 assert( DVector::dimension ==
dimension);
761 evaluate<deriv>( bx, values );
765 template <
class Vector>
775 template<
class Topology,
class F >
796 Field *
const values )
const
814 template<
int dim,
class F >
823 template <
int dd,
class FF >
829 template<
class Topology >
842 template<
int dim,
class SF >
844 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
848 template <
int dd,
class FF >
Definition: bdfmcube.hh:16
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
@ value
Definition: tensor.hh:166
A class representing the unit of a given Field.
Definition: field.hh:28
A class representing the zero of a given Field.
Definition: field.hh:77
Definition: monomialbasis.hh:82
unsigned int * sizes_
Definition: monomialbasis.hh:95
~MonomialBasisSize()
Definition: monomialbasis.hh:108
unsigned int * numBaseFunctions_
Definition: monomialbasis.hh:98
MonomialBasisSize()
Definition: monomialbasis.hh:100
static This & instance()
Definition: monomialbasis.hh:86
void computeSizes(unsigned int order)
Definition: monomialbasis.hh:124
unsigned int operator()(const unsigned int order) const
Definition: monomialbasis.hh:114
unsigned int maxOrder() const
Definition: monomialbasis.hh:119
unsigned int maxOrder_
Definition: monomialbasis.hh:92
Definition: monomialbasis.hh:506
Base::Field Field
Definition: monomialbasis.hh:514
static const unsigned int dimension
Definition: monomialbasis.hh:511
unsigned int size() const
Definition: monomialbasis.hh:541
unsigned int derivSize(const unsigned int deriv) const
Definition: monomialbasis.hh:547
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:584
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:571
void integrate(Vector &values) const
Definition: monomialbasis.hh:618
unsigned int topologyId() const
Definition: monomialbasis.hh:559
Dune::FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:518
void evaluate(const DomainVector &x, FieldVector< Field, Derivatives< Field, dimension, 1, deriv, DerivativeLayoutNS::value >::size > *values) const
Definition: monomialbasis.hh:590
unsigned int order() const
Definition: monomialbasis.hh:554
MonomialBasis(unsigned int order)
Definition: monomialbasis.hh:522
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:597
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:564
Base::DomainVector DomainVector
Definition: monomialbasis.hh:516
void integrate(Field *const values) const
Definition: monomialbasis.hh:613
static const unsigned int dimRange
Definition: monomialbasis.hh:512
const unsigned int * sizes() const
Definition: monomialbasis.hh:536
MonomialBasisSize< Topology > Size
Definition: monomialbasis.hh:520
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:604
const unsigned int * sizes(unsigned int order) const
Definition: monomialbasis.hh:530
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:578
Definition: monomialbasis.hh:174
MonomialBasisSize< typename Impl::SimplexTopology< mydim >::type > MySize
Definition: monomialbasis.hh:175
MonomialBasisSize< typename Impl::SimplexTopology< dim >::type > Size
Definition: monomialbasis.hh:176
static void copy(const unsigned int deriv, F *&wit, F *&rit, const unsigned int numBaseFunctions, const F &z)
Definition: monomialbasis.hh:178
Definition: monomialbasis.hh:230
Definition: monomialbasis.hh:234
F Field
Definition: monomialbasis.hh:239
Impl::Point Topology
Definition: monomialbasis.hh:238
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:243
Definition: monomialbasis.hh:272
F Field
Definition: monomialbasis.hh:277
Impl::Prism< BaseTopology > Topology
Definition: monomialbasis.hh:276
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:281
Definition: monomialbasis.hh:357
F Field
Definition: monomialbasis.hh:362
Impl::Pyramid< BaseTopology > Topology
Definition: monomialbasis.hh:361
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:366
Definition: monomialbasis.hh:637
StandardMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:645
Impl::SimplexTopology< dim >::type Topology
Definition: monomialbasis.hh:642
static const int dimension
Definition: monomialbasis.hh:643
Definition: monomialbasis.hh:658
static const int dimension
Definition: monomialbasis.hh:664
Impl::CubeTopology< dim >::type Topology
Definition: monomialbasis.hh:663
StandardBiMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:666
Definition: monomialbasis.hh:678
unsigned int topologyId() const
Definition: monomialbasis.hh:708
FieldVector< Field, dimension > DomainVector
Definition: monomialbasis.hh:687
unsigned int order_
Definition: monomialbasis.hh:771
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:740
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:716
F Field
Definition: monomialbasis.hh:682
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:734
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:746
unsigned int order() const
Definition: monomialbasis.hh:703
virtual const unsigned int * sizes() const =0
unsigned int topologyId_
Definition: monomialbasis.hh:772
static const unsigned int dimRange
Definition: monomialbasis.hh:685
F StorageField
Definition: monomialbasis.hh:683
static const int dimension
Definition: monomialbasis.hh:684
VirtualMonomialBasis(unsigned int topologyId, unsigned int order)
Definition: monomialbasis.hh:690
unsigned int size() const
Definition: monomialbasis.hh:698
FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:688
virtual ~VirtualMonomialBasis()
Definition: monomialbasis.hh:694
virtual void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const =0
virtual void integrate(Field *const values) const =0
void evaluate(const DomainVector &x, Dune::FieldVector< Field, size > *const values) const
Definition: monomialbasis.hh:722
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:755
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:728
void integrate(Vector &values) const
Definition: monomialbasis.hh:766
Definition: monomialbasis.hh:778
VirtualMonomialBasisImpl(unsigned int order)
Definition: monomialbasis.hh:786
Base::DomainVector DomainVector
Definition: monomialbasis.hh:784
void integrate(Field *const values) const
Definition: monomialbasis.hh:801
const unsigned int * sizes() const
Definition: monomialbasis.hh:790
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:795
Base::Field Field
Definition: monomialbasis.hh:783
Definition: monomialbasis.hh:816
static void release(Object *object)
Definition: monomialbasis.hh:834
const VirtualMonomialBasis< dimension, F > Object
Definition: monomialbasis.hh:821
F StorageField
Definition: monomialbasis.hh:818
static const unsigned int dimension
Definition: monomialbasis.hh:817
static Object * create(const Key &order)
Definition: monomialbasis.hh:830
unsigned int Key
Definition: monomialbasis.hh:820
Definition: monomialbasis.hh:825
MonomialBasisFactory< dd, FF > Type
Definition: monomialbasis.hh:826
Definition: monomialbasis.hh:845
static const unsigned int dimension
Definition: monomialbasis.hh:846
SF StorageField
Definition: monomialbasis.hh:847
Definition: monomialbasis.hh:850
MonomialBasisProvider< dd, FF > Type
Definition: monomialbasis.hh:851
Definition: tensor.hh:170