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/type.hh>
12 #include <dune/geometry/topologyfactory.hh>
69 template< GeometryType::Id geometryId >
70 class MonomialBasisSize;
72 template< GeometryType::Id geometryId,
class F >
80 template< GeometryType::Id geometryId >
88 static This _instance;
133 sizes_ =
new unsigned int[ order+1 ];
136 constexpr GeometryType geometry = geometryId;
137 constexpr
auto dim = geometry.dim();
140 for(
unsigned int k = 1; k <= order; ++k )
145 for(
int codim=dim-1; codim>=0; codim--)
147 if (Impl::isPrism(geometry.id(),dim,codim))
149 for(
unsigned int k = 1; k <= order; ++k )
157 for(
unsigned int k = 1; k <= order; ++k )
173 template<
int mydim,
int dim,
class F >
179 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
180 const unsigned int numBaseFunctions,
const F &z )
186 const F *
const rend = rit + size( deriv )*numBaseFunctions;
187 for( ; rit != rend; )
194 for(
unsigned d = 1; d <= deriv; ++d )
197 const F *
const derivEnd = rit + mySize.
sizes_[ d ];
201 const F *
const drend = rit + mySize.
sizes_[ d ] - mySize.
sizes_[ d-1 ];
202 for( ; rit != drend ; ++rit, ++wit )
206 for (
unsigned int j=1; j<d; ++j)
208 const F *
const drend = rit + mySize.
sizes_[ d-j ] - mySize.
sizes_[ d-j-1 ];
209 for( ; rit != drend ; ++prit, ++rit, ++wit )
210 *wit = F(j) * *prit + z * *rit;
212 *wit = F(d) * *prit + z * *rit;
213 ++prit, ++rit, ++wit;
214 assert(derivEnd == rit);
217 const F *
const emptyWitEnd = wit + size.
sizes_[d] - mySize.
sizes_[d];
218 for ( ; wit != emptyWitEnd; ++wit )
230 template< GeometryType::Id geometryId,
class F>
236 static constexpr GeometryType
geometry = geometryId;
249 void evaluate (
const unsigned int deriv,
const unsigned int order,
250 const FieldVector< Field, dimD > &x,
251 const unsigned int block,
const unsigned int *
const offsets,
252 Field *
const values )
const
256 F *
const end = values + block;
257 for(
Field *it = values+1 ; it != end; ++it )
260 constexpr GeometryType gt = GeometryTypes::vertex;
265 evaluate<gt,dimD>(deriv, order, x, block, offsets, values );
268 template<GeometryType::Id baseGeometryId,
int dimD >
269 void evaluate (
const unsigned int deriv,
const unsigned int order,
270 const FieldVector< Field, dimD > &x,
271 const unsigned int block,
const unsigned int *
const offsets,
272 Field *
const values )
const
275 static constexpr GeometryType baseGeometry = baseGeometryId;
277 auto constexpr isPrismatic =
geometry.isPrismatic(baseGeometry.dim());
280 typedef MonomialBasisHelper< baseGeometry.dim() + 1, dimD,
Field > Helper;
281 typedef MonomialBasisSize<baseGeometryId> BaseSize;
283 const BaseSize &size = BaseSize::instance();
284 const_cast<BaseSize&
>(size).computeSizes(order);
286 const Field &z = x[ baseGeometry.dim() ];
288 Field *row0 = values;
289 for(
unsigned int k = 1; k <= order; ++k )
291 Field *row1 = values + block*offsets[ k-1 ];
292 Field *wit = row1 + block*size.sizes_[ k ];
293 if constexpr ( isPrismatic )
294 Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
295 Helper::copy( deriv, wit, row0, size( k-1 ), z );
300 if constexpr( baseGeometry.dim() ==
dimDomain-1)
304 constexpr GeometryType nextGeometry = isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry)
305 : GeometryTypes::conicalExtension(baseGeometry);
307 evaluate<nextGeometry.toId(),dimD>(deriv, order, x, block, offsets, values );
311 void integrate (
const unsigned int order,
312 const unsigned int *
const offsets,
313 Field *
const values )
const
316 values[ 0 ] = Unity< Field >();
317 static constexpr GeometryType gt = GeometryTypes::vertex;
322 integrate<gt>(order, offsets, values);
325 template<GeometryType::Id baseGeometryId>
326 void integrate (
const unsigned int order,
327 const unsigned int *
const offsets,
328 Field *
const values)
const
330 static constexpr GeometryType baseGeometry = baseGeometryId;
332 auto constexpr isPrismatic =
geometry.isPrismatic(baseGeometry.dim());
335 if constexpr ( isPrismatic )
336 integratePrismatic<baseGeometry>(order, offsets, values);
338 integrateConical<baseGeometry>(order, offsets, values);
341 if constexpr( baseGeometry.dim() ==
dimDomain-1)
345 static constexpr GeometryType nextGeometry = (isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry)
346 : GeometryTypes::conicalExtension(baseGeometry));
348 integrate<nextGeometry.toId()>(order, offsets, values);
353 template<GeometryType::Id baseGeometryId>
354 void integratePrismatic (
const unsigned int order,
355 const unsigned int *
const offsets,
356 Field *
const values )
const
358 typedef MonomialBasisSize<baseGeometryId> BaseSize;
359 static const BaseSize &size = BaseSize::instance();
360 const unsigned int *
const baseSizes = size.sizes_;
362 static constexpr GeometryType baseGeometry = baseGeometryId;
363 static constexpr GeometryType nextGeometry = GeometryTypes::prismaticExtension(baseGeometry);
365 typedef MonomialBasisSize<nextGeometry.toId()> Size;
366 static const Size &mySize = Size::instance();
368 Field *row0 = values;
369 for(
unsigned int k = 1; k <= order; ++k )
371 Field *
const row1begin = values + offsets[ k-1 ];
372 Field *
const row1End = row1begin + mySize.sizes_[ k ];
373 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
375 Field *row1 = row1begin;
376 Field *it = row1begin + baseSizes[ k ];
377 for(
unsigned int j = 1; j <= k; ++j )
379 Field *
const end = it + baseSizes[ k ];
380 assert( (
unsigned int)(end - values) <= offsets[ k ] );
381 for( ; it != end; ++row1, ++it )
384 for( ; it != row1End; ++row0, ++it )
391 template<GeometryType::Id baseGeometryId>
392 void integrateConical (
const unsigned int order,
393 const unsigned int *
const offsets,
394 Field *
const values)
const
396 typedef MonomialBasisSize<baseGeometryId> BaseSize;
397 static const BaseSize &size = BaseSize::instance();
398 const unsigned int *
const baseSizes = size.sizes_;
400 static constexpr GeometryType baseGeometry = baseGeometryId;
403 Field *
const col0End = values + baseSizes[ 0 ];
404 for(
Field *it = values; it != col0End; ++it )
405 *it *=
Field( 1 ) /
Field(
int(baseGeometry.dim()+1) );
408 Field *row0 = values;
409 for(
unsigned int k = 1; k <= order; ++k )
413 Field *
const row1 = values+offsets[ k-1 ];
414 Field *
const col0End = row1 + baseSizes[ k ];
416 for( ; it != col0End; ++it )
418 for(
unsigned int i = 1; i <= k; ++i )
420 Field *
const end = it + baseSizes[ k-i ];
421 assert( (
unsigned int)(end - values) <= offsets[ k ] );
422 for( ; it != end; ++row0, ++it )
423 *it = (*row0) * (
Field( i ) * factor);
435 template< GeometryType::Id geometryId,
class F >
439 static constexpr GeometryType geometry = geometryId;
458 size_(
Size::instance())
471 return sizes( order_ );
477 return size_( order_ );
480 unsigned int derivSize (
const unsigned int deriv )
const
493 return geometry.id();
497 Field *
const values )
const
499 Base::evaluate( deriv, order_, x,
derivSize( deriv ),
sizes( order_ ), values );
502 template <
unsigned int deriv>
504 Field *
const values )
const
509 template<
unsigned int deriv,
class Vector >
511 Vector &values )
const
513 evaluate<deriv>(x,&(values[0]));
515 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
519 evaluate<deriv>(x,&(values->block()));
521 template<
unsigned int deriv >
528 template<
class Vector >
530 Vector &values )
const
532 evaluate<0>(x,&(values[0]));
535 template<
class DVector,
class RVector >
536 void evaluate (
const DVector &x, RVector &values )
const
538 assert( DVector::dimension ==
dimension);
542 evaluate<0>( bx, values );
547 Base::integrate( order_,
sizes( order_ ), values );
549 template <
class Vector>
556 This& operator=(
const This&);
566 template<
int dim,
class F >
568 :
public MonomialBasis< GeometryTypes::simplex(dim).toId() , F >
574 static constexpr GeometryType
geometry = GeometryTypes::simplex(dim);
587 template<
int dim,
class F >
589 :
public MonomialBasis< GeometryTypes::cube(dim).toId() , F >
595 static constexpr GeometryType
geometry = GeometryTypes::cube(dim);
608 template<
int dim,
class F >
622 [[deprecated(
"Use VirtualMonomialBasis(GeometryType gt, unsigned int order) instead.")]]
633 virtual const unsigned int *
sizes ( )
const = 0;
645 [[deprecated(
"Use type().id() instead.")]]
657 Field *
const values )
const = 0;
658 template <
unsigned int deriv >
660 Field *
const values )
const
664 template <
unsigned int deriv,
int size >
666 Dune::FieldVector<Field,size> *
const values )
const
668 evaluate( deriv, x, &(values[0][0]) );
670 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
674 evaluate<deriv>(x,&(values->block()));
676 template <
unsigned int deriv,
class Vector>
678 Vector &values )
const
680 evaluate<deriv>( x, &(values[ 0 ]) );
682 template<
class Vector >
684 Vector &values )
const
686 evaluate<0>(x,values);
688 template<
class DVector,
class RVector >
689 void evaluate (
const DVector &x, RVector &values )
const
691 assert( DVector::dimension ==
dimension);
695 evaluate<0>( bx, values );
697 template<
unsigned int deriv,
class DVector,
class RVector >
698 void evaluate (
const DVector &x, RVector &values )
const
700 assert( DVector::dimension ==
dimension);
704 evaluate<deriv>( bx, values );
708 template <
class Vector>
718 template< GeometryType::Id geometryId,
class F >
722 static constexpr GeometryType geometry = geometryId;
736 return basis_.
sizes(order_);
740 Field *
const values )
const
758 template<
int dim,
class F >
767 template <
int dd,
class FF >
773 template< GeometryType::Id geometryId >
786 template<
int dim,
class SF >
788 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
792 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
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 * numBaseFunctions_
Definition: monomialbasis.hh:98
void computeSizes(unsigned int order)
Definition: monomialbasis.hh:124
unsigned int * sizes_
Definition: monomialbasis.hh:95
MonomialBasisSize()
Definition: monomialbasis.hh:100
static This & instance()
Definition: monomialbasis.hh:86
~MonomialBasisSize()
Definition: monomialbasis.hh:108
unsigned int operator()(const unsigned int order) const
Definition: monomialbasis.hh:114
unsigned int maxOrder_
Definition: monomialbasis.hh:92
unsigned int maxOrder() const
Definition: monomialbasis.hh:119
Definition: monomialbasis.hh:438
unsigned int size() const
Definition: monomialbasis.hh:474
static const unsigned int dimension
Definition: monomialbasis.hh:444
Dune::FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:451
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:496
const unsigned int * sizes() const
Definition: monomialbasis.hh:469
unsigned int topologyId() const
Definition: monomialbasis.hh:491
void integrate(Vector &values) const
Definition: monomialbasis.hh:550
Base::Field Field
Definition: monomialbasis.hh:447
void integrate(Field *const values) const
Definition: monomialbasis.hh:545
static const unsigned int dimRange
Definition: monomialbasis.hh:445
Base::DomainVector DomainVector
Definition: monomialbasis.hh:449
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:503
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:529
void evaluate(const DomainVector &x, FieldVector< Field, Derivatives< Field, dimension, 1, deriv, DerivativeLayoutNS::value >::size > *values) const
Definition: monomialbasis.hh:522
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:516
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:536
MonomialBasisSize< geometryId > Size
Definition: monomialbasis.hh:453
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:510
MonomialBasis(unsigned int order)
Definition: monomialbasis.hh:455
unsigned int derivSize(const unsigned int deriv) const
Definition: monomialbasis.hh:480
const unsigned int * sizes(unsigned int order) const
Definition: monomialbasis.hh:463
unsigned int order() const
Definition: monomialbasis.hh:486
Definition: monomialbasis.hh:175
MonomialBasisSize< GeometryTypes::simplex(dim).toId() > Size
Definition: monomialbasis.hh:177
static void copy(const unsigned int deriv, F *&wit, F *&rit, const unsigned int numBaseFunctions, const F &z)
Definition: monomialbasis.hh:179
MonomialBasisSize< GeometryTypes::simplex(mydim).toId() > MySize
Definition: monomialbasis.hh:176
Definition: monomialbasis.hh:232
FieldVector< Field, dimDomain > DomainVector
Definition: monomialbasis.hh:240
static constexpr GeometryType geometry
Definition: monomialbasis.hh:236
F Field
Definition: monomialbasis.hh:234
static const unsigned int dimDomain
Definition: monomialbasis.hh:238
Definition: monomialbasis.hh:569
static constexpr GeometryType geometry
Definition: monomialbasis.hh:574
StandardMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:577
static const int dimension
Definition: monomialbasis.hh:575
Definition: monomialbasis.hh:590
static const int dimension
Definition: monomialbasis.hh:596
static constexpr GeometryType geometry
Definition: monomialbasis.hh:595
StandardBiMonomialBasis(unsigned int order)
Definition: monomialbasis.hh:598
Definition: monomialbasis.hh:610
unsigned int topologyId() const
Definition: monomialbasis.hh:646
GeometryType geometry_
Definition: monomialbasis.hh:715
FieldVector< Field, dimension > DomainVector
Definition: monomialbasis.hh:619
unsigned int order_
Definition: monomialbasis.hh:714
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:683
void evaluate(const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:659
F Field
Definition: monomialbasis.hh:614
void evaluate(const DomainVector &x, Vector &values) const
Definition: monomialbasis.hh:677
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:689
unsigned int order() const
Definition: monomialbasis.hh:640
virtual const unsigned int * sizes() const =0
static const unsigned int dimRange
Definition: monomialbasis.hh:617
F StorageField
Definition: monomialbasis.hh:615
static const int dimension
Definition: monomialbasis.hh:616
VirtualMonomialBasis(unsigned int topologyId, unsigned int order)
Definition: monomialbasis.hh:623
unsigned int size() const
Definition: monomialbasis.hh:635
FieldVector< Field, dimRange > RangeVector
Definition: monomialbasis.hh:620
virtual ~VirtualMonomialBasis()
Definition: monomialbasis.hh:631
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:665
void evaluate(const DVector &x, RVector &values) const
Definition: monomialbasis.hh:698
GeometryType type() const
Definition: monomialbasis.hh:651
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition: monomialbasis.hh:671
VirtualMonomialBasis(const GeometryType >, unsigned int order)
Definition: monomialbasis.hh:627
void integrate(Vector &values) const
Definition: monomialbasis.hh:709
Definition: monomialbasis.hh:721
void integrate(Field *const values) const
Definition: monomialbasis.hh:745
const unsigned int * sizes() const
Definition: monomialbasis.hh:734
Base::DomainVector DomainVector
Definition: monomialbasis.hh:728
Base::Field Field
Definition: monomialbasis.hh:727
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition: monomialbasis.hh:739
VirtualMonomialBasisImpl(unsigned int order)
Definition: monomialbasis.hh:730
Definition: monomialbasis.hh:760
static void release(Object *object)
Definition: monomialbasis.hh:778
const VirtualMonomialBasis< dimension, F > Object
Definition: monomialbasis.hh:765
F StorageField
Definition: monomialbasis.hh:762
static const unsigned int dimension
Definition: monomialbasis.hh:761
static Object * create(const Key &order)
Definition: monomialbasis.hh:774
unsigned int Key
Definition: monomialbasis.hh:764
Definition: monomialbasis.hh:769
MonomialBasisFactory< dd, FF > Type
Definition: monomialbasis.hh:770
Definition: monomialbasis.hh:789
static const unsigned int dimension
Definition: monomialbasis.hh:790
SF StorageField
Definition: monomialbasis.hh:791
Definition: monomialbasis.hh:794
MonomialBasisProvider< dd, FF > Type
Definition: monomialbasis.hh:795
Definition: tensor.hh:170