5#ifndef DUNE_MONOMIALBASIS_HH
6#define DUNE_MONOMIALBASIS_HH
10#include <dune/common/fvector.hh>
11#include <dune/common/fmatrix.hh>
13#include <dune/geometry/type.hh>
14#include <dune/geometry/topologyfactory.hh>
71 template< GeometryType::Id geometryId >
72 class MonomialBasisSize;
74 template< GeometryType::Id geometryId,
class F >
82 template< GeometryType::Id geometryId >
90 static This _instance;
135 sizes_ =
new unsigned int[ order+1 ];
138 constexpr GeometryType geometry = geometryId;
139 constexpr auto dim = geometry.dim();
142 for(
unsigned int k = 1; k <= order; ++k )
147 for(
int codim=dim-1; codim>=0; codim--)
149 if (Impl::isPrism(geometry.id(),dim,codim))
151 for(
unsigned int k = 1; k <= order; ++k )
159 for(
unsigned int k = 1; k <= order; ++k )
175 template<
int mydim,
int dim,
class F >
181 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
182 const unsigned int numBaseFunctions,
const F &z )
188 const F *
const rend = rit + size( deriv )*numBaseFunctions;
189 for( ; rit != rend; )
196 for(
unsigned d = 1; d <= deriv; ++d )
199 const F *
const derivEnd = rit + mySize.
sizes_[ d ];
203 const F *
const drend = rit + mySize.
sizes_[ d ] - mySize.
sizes_[ d-1 ];
204 for( ; rit != drend ; ++rit, ++wit )
208 for (
unsigned int j=1; j<d; ++j)
210 const F *
const drend = rit + mySize.
sizes_[ d-j ] - mySize.
sizes_[ d-j-1 ];
211 for( ; rit != drend ; ++prit, ++rit, ++wit )
212 *wit = F(j) * *prit + z * *rit;
214 *wit = F(d) * *prit + z * *rit;
215 ++prit, ++rit, ++wit;
216 assert(derivEnd == rit);
219 const F *
const emptyWitEnd = wit + size.
sizes_[d] - mySize.
sizes_[d];
220 for ( ; wit != emptyWitEnd; ++wit )
232 template< GeometryType::Id geometryId,
class F>
238 static constexpr GeometryType
geometry = geometryId;
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
258 F *
const end = values + block;
259 for(
Field *it = values+1 ; it != end; ++it )
262 constexpr GeometryType gt = GeometryTypes::vertex;
267 evaluate<gt,dimD>(deriv, order, x, block, offsets, values );
270 template<GeometryType::Id baseGeometryId,
int dimD >
271 void evaluate (
const unsigned int deriv,
const unsigned int order,
272 const FieldVector< Field, dimD > &x,
273 const unsigned int block,
const unsigned int *
const offsets,
274 Field *
const values )
const
277 static constexpr GeometryType baseGeometry = baseGeometryId;
279 auto constexpr isPrismatic =
geometry.isPrismatic(baseGeometry.dim());
282 typedef MonomialBasisHelper< baseGeometry.dim() + 1, dimD,
Field > Helper;
283 typedef MonomialBasisSize<baseGeometryId> BaseSize;
285 const BaseSize &size = BaseSize::instance();
286 const_cast<BaseSize&
>(size).computeSizes(order);
288 const Field &z = x[ baseGeometry.dim() ];
290 Field *row0 = values;
291 for(
unsigned int k = 1; k <= order; ++k )
293 Field *row1 = values + block*offsets[ k-1 ];
294 Field *wit = row1 + block*size.sizes_[ k ];
295 if constexpr ( isPrismatic )
296 Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
297 Helper::copy( deriv, wit, row0, size( k-1 ), z );
302 if constexpr( baseGeometry.dim() ==
dimDomain-1)
306 constexpr GeometryType nextGeometry = isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry)
307 : GeometryTypes::conicalExtension(baseGeometry);
309 evaluate<nextGeometry.toId(),dimD>(deriv, order, x, block, offsets, values );
313 void integrate (
const unsigned int order,
314 const unsigned int *
const offsets,
315 Field *
const values )
const
318 values[ 0 ] = Unity< Field >();
319 static constexpr GeometryType gt = GeometryTypes::vertex;
324 integrate<gt>(order, offsets, values);
327 template<GeometryType::Id baseGeometryId>
328 void integrate (
const unsigned int order,
329 const unsigned int *
const offsets,
330 Field *
const values)
const
332 static constexpr GeometryType baseGeometry = baseGeometryId;
334 auto constexpr isPrismatic =
geometry.isPrismatic(baseGeometry.dim());
337 if constexpr ( isPrismatic )
338 integratePrismatic<baseGeometry>(order, offsets, values);
340 integrateConical<baseGeometry>(order, offsets, values);
343 if constexpr( baseGeometry.dim() ==
dimDomain-1)
347 static constexpr GeometryType nextGeometry = (isPrismatic ? GeometryTypes::prismaticExtension(baseGeometry)
348 : GeometryTypes::conicalExtension(baseGeometry));
350 integrate<nextGeometry.toId()>(order, offsets, values);
355 template<GeometryType::Id baseGeometryId>
356 void integratePrismatic (
const unsigned int order,
357 const unsigned int *
const offsets,
358 Field *
const values )
const
360 typedef MonomialBasisSize<baseGeometryId> BaseSize;
361 static const BaseSize &size = BaseSize::instance();
362 const unsigned int *
const baseSizes = size.sizes_;
364 static constexpr GeometryType baseGeometry = baseGeometryId;
365 static constexpr GeometryType nextGeometry = GeometryTypes::prismaticExtension(baseGeometry);
367 typedef MonomialBasisSize<nextGeometry.toId()> Size;
368 static const Size &mySize = Size::instance();
370 Field *row0 = values;
371 for(
unsigned int k = 1; k <= order; ++k )
373 Field *
const row1begin = values + offsets[ k-1 ];
374 Field *
const row1End = row1begin + mySize.sizes_[ k ];
375 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
377 Field *row1 = row1begin;
378 Field *it = row1begin + baseSizes[ k ];
379 for(
unsigned int j = 1; j <= k; ++j )
381 Field *
const end = it + baseSizes[ k ];
382 assert( (
unsigned int)(end - values) <= offsets[ k ] );
383 for( ; it != end; ++row1, ++it )
386 for( ; it != row1End; ++row0, ++it )
393 template<GeometryType::Id baseGeometryId>
394 void integrateConical (
const unsigned int order,
395 const unsigned int *
const offsets,
396 Field *
const values)
const
398 typedef MonomialBasisSize<baseGeometryId> BaseSize;
399 static const BaseSize &size = BaseSize::instance();
400 const unsigned int *
const baseSizes = size.sizes_;
402 static constexpr GeometryType baseGeometry = baseGeometryId;
405 Field *
const col0End = values + baseSizes[ 0 ];
406 for(
Field *it = values; it != col0End; ++it )
407 *it *=
Field( 1 ) /
Field(
int(baseGeometry.dim()+1) );
410 Field *row0 = values;
411 for(
unsigned int k = 1; k <= order; ++k )
415 Field *
const row1 = values+offsets[ k-1 ];
416 Field *
const col0End = row1 + baseSizes[ k ];
418 for( ; it != col0End; ++it )
420 for(
unsigned int i = 1; i <= k; ++i )
422 Field *
const end = it + baseSizes[ k-i ];
423 assert( (
unsigned int)(end - values) <= offsets[ k ] );
424 for( ; it != end; ++row0, ++it )
425 *it = (*row0) * (
Field( i ) * factor);
437 template< GeometryType::Id geometryId,
class F >
441 static constexpr GeometryType geometry = geometryId;
460 size_(
Size::instance())
473 return sizes( order_ );
479 return size_( order_ );
482 unsigned int derivSize (
const unsigned int deriv )
const
495 return geometry.id();
499 Field *
const values )
const
501 Base::evaluate( deriv, order_, x,
derivSize( deriv ),
sizes( order_ ), values );
504 template <
unsigned int deriv>
506 Field *
const values )
const
511 template<
unsigned int deriv,
class Vector >
513 Vector &values )
const
515 evaluate<deriv>(x,&(values[0]));
517 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
521 evaluate<deriv>(x,&(values->block()));
523 template<
unsigned int deriv >
530 template<
class Vector >
532 Vector &values )
const
534 evaluate<0>(x,&(values[0]));
537 template<
class DVector,
class RVector >
538 void evaluate (
const DVector &x, RVector &values )
const
540 assert( DVector::dimension ==
dimension);
544 evaluate<0>( bx, values );
549 Base::integrate( order_,
sizes( order_ ), values );
551 template <
class Vector>
558 This& operator=(
const This&);
568 template<
int dim,
class F >
570 :
public MonomialBasis< GeometryTypes::simplex(dim).toId() , F >
576 static constexpr GeometryType
geometry = GeometryTypes::simplex(dim);
589 template<
int dim,
class F >
591 :
public MonomialBasis< GeometryTypes::cube(dim).toId() , F >
597 static constexpr GeometryType
geometry = GeometryTypes::cube(dim);
610 template<
int dim,
class F >
630 virtual const unsigned int *
sizes ( )
const = 0;
648 Field *
const values )
const = 0;
649 template <
unsigned int deriv >
651 Field *
const values )
const
655 template <
unsigned int deriv,
int size >
657 Dune::FieldVector<Field,size> *
const values )
const
659 evaluate( deriv, x, &(values[0][0]) );
661 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
665 evaluate<deriv>(x,&(values->block()));
667 template <
unsigned int deriv,
class Vector>
669 Vector &values )
const
671 evaluate<deriv>( x, &(values[ 0 ]) );
673 template<
class Vector >
675 Vector &values )
const
677 evaluate<0>(x,values);
679 template<
class DVector,
class RVector >
680 void evaluate (
const DVector &x, RVector &values )
const
682 assert( DVector::dimension ==
dimension);
686 evaluate<0>( bx, values );
688 template<
unsigned int deriv,
class DVector,
class RVector >
689 void evaluate (
const DVector &x, RVector &values )
const
691 assert( DVector::dimension ==
dimension);
695 evaluate<deriv>( bx, values );
699 template <
class Vector>
709 template< GeometryType::Id geometryId,
class F >
713 static constexpr GeometryType geometry = geometryId;
727 return basis_.
sizes(order_);
731 Field *
const values )
const
749 template<
int dim,
class F >
758 template <
int dd,
class FF >
764 template< GeometryType::Id geometryId >
777 template<
int dim,
class SF >
779 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
783 template <
int dd,
class FF >
Definition bdfmcube.hh:18
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition field.hh:159
A class representing the unit of a given Field.
Definition field.hh:30
A class representing the zero of a given Field.
Definition field.hh:79
Definition monomialbasis.hh:84
unsigned int * numBaseFunctions_
Definition monomialbasis.hh:100
void computeSizes(unsigned int order)
Definition monomialbasis.hh:126
unsigned int * sizes_
Definition monomialbasis.hh:97
MonomialBasisSize()
Definition monomialbasis.hh:102
~MonomialBasisSize()
Definition monomialbasis.hh:110
unsigned int operator()(const unsigned int order) const
Definition monomialbasis.hh:116
unsigned int maxOrder_
Definition monomialbasis.hh:94
unsigned int maxOrder() const
Definition monomialbasis.hh:121
static This & instance()
Definition monomialbasis.hh:88
Definition monomialbasis.hh:440
unsigned int size() const
Definition monomialbasis.hh:476
static const unsigned int dimension
Definition monomialbasis.hh:446
Dune::FieldVector< Field, dimRange > RangeVector
Definition monomialbasis.hh:453
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition monomialbasis.hh:498
const unsigned int * sizes() const
Definition monomialbasis.hh:471
unsigned int topologyId() const
Definition monomialbasis.hh:493
void integrate(Vector &values) const
Definition monomialbasis.hh:552
Base::Field Field
Definition monomialbasis.hh:449
void integrate(Field *const values) const
Definition monomialbasis.hh:547
static const unsigned int dimRange
Definition monomialbasis.hh:447
Base::DomainVector DomainVector
Definition monomialbasis.hh:451
void evaluate(const DomainVector &x, Field *const values) const
Definition monomialbasis.hh:505
void evaluate(const DomainVector &x, Vector &values) const
Definition monomialbasis.hh:531
void evaluate(const DomainVector &x, FieldVector< Field, Derivatives< Field, dimension, 1, deriv, DerivativeLayoutNS::value >::size > *values) const
Definition monomialbasis.hh:524
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition monomialbasis.hh:518
void evaluate(const DVector &x, RVector &values) const
Definition monomialbasis.hh:538
MonomialBasisSize< geometryId > Size
Definition monomialbasis.hh:455
void evaluate(const DomainVector &x, Vector &values) const
Definition monomialbasis.hh:512
MonomialBasis(unsigned int order)
Definition monomialbasis.hh:457
unsigned int derivSize(const unsigned int deriv) const
Definition monomialbasis.hh:482
const unsigned int * sizes(unsigned int order) const
Definition monomialbasis.hh:465
unsigned int order() const
Definition monomialbasis.hh:488
Definition monomialbasis.hh:177
MonomialBasisSize< GeometryTypes::simplex(dim).toId() > Size
Definition monomialbasis.hh:179
static void copy(const unsigned int deriv, F *&wit, F *&rit, const unsigned int numBaseFunctions, const F &z)
Definition monomialbasis.hh:181
MonomialBasisSize< GeometryTypes::simplex(mydim).toId() > MySize
Definition monomialbasis.hh:178
Definition monomialbasis.hh:234
FieldVector< Field, dimDomain > DomainVector
Definition monomialbasis.hh:242
static constexpr GeometryType geometry
Definition monomialbasis.hh:238
F Field
Definition monomialbasis.hh:236
static const unsigned int dimDomain
Definition monomialbasis.hh:240
Definition monomialbasis.hh:571
static constexpr GeometryType geometry
Definition monomialbasis.hh:576
StandardMonomialBasis(unsigned int order)
Definition monomialbasis.hh:579
static const int dimension
Definition monomialbasis.hh:577
Definition monomialbasis.hh:592
static const int dimension
Definition monomialbasis.hh:598
static constexpr GeometryType geometry
Definition monomialbasis.hh:597
StandardBiMonomialBasis(unsigned int order)
Definition monomialbasis.hh:600
Definition monomialbasis.hh:612
GeometryType geometry_
Definition monomialbasis.hh:706
FieldVector< Field, dimension > DomainVector
Definition monomialbasis.hh:621
unsigned int order_
Definition monomialbasis.hh:705
void evaluate(const DomainVector &x, Vector &values) const
Definition monomialbasis.hh:674
void evaluate(const DomainVector &x, Field *const values) const
Definition monomialbasis.hh:650
F Field
Definition monomialbasis.hh:616
void evaluate(const DomainVector &x, Vector &values) const
Definition monomialbasis.hh:668
void evaluate(const DVector &x, RVector &values) const
Definition monomialbasis.hh:680
unsigned int order() const
Definition monomialbasis.hh:637
static const unsigned int dimRange
Definition monomialbasis.hh:619
F StorageField
Definition monomialbasis.hh:617
static const int dimension
Definition monomialbasis.hh:618
unsigned int size() const
Definition monomialbasis.hh:632
FieldVector< Field, dimRange > RangeVector
Definition monomialbasis.hh:622
virtual ~VirtualMonomialBasis()
Definition monomialbasis.hh:628
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:656
void evaluate(const DVector &x, RVector &values) const
Definition monomialbasis.hh:689
GeometryType type() const
Definition monomialbasis.hh:642
void evaluate(const DomainVector &x, Derivatives< Field, dimension, 1, deriv, layout > *values) const
Definition monomialbasis.hh:662
VirtualMonomialBasis(const GeometryType >, unsigned int order)
Definition monomialbasis.hh:624
virtual const unsigned int * sizes() const =0
void integrate(Vector &values) const
Definition monomialbasis.hh:700
Definition monomialbasis.hh:712
void integrate(Field *const values) const
Definition monomialbasis.hh:736
const unsigned int * sizes() const
Definition monomialbasis.hh:725
Base::DomainVector DomainVector
Definition monomialbasis.hh:719
Base::Field Field
Definition monomialbasis.hh:718
void evaluate(const unsigned int deriv, const DomainVector &x, Field *const values) const
Definition monomialbasis.hh:730
VirtualMonomialBasisImpl(unsigned int order)
Definition monomialbasis.hh:721
Definition monomialbasis.hh:751
static void release(Object *object)
Definition monomialbasis.hh:769
const VirtualMonomialBasis< dimension, F > Object
Definition monomialbasis.hh:756
static Object * create(const Key &order)
Definition monomialbasis.hh:765
F StorageField
Definition monomialbasis.hh:753
static const unsigned int dimension
Definition monomialbasis.hh:752
unsigned int Key
Definition monomialbasis.hh:755
Definition monomialbasis.hh:760
MonomialBasisFactory< dd, FF > Type
Definition monomialbasis.hh:761
Definition monomialbasis.hh:780
static const unsigned int dimension
Definition monomialbasis.hh:781
SF StorageField
Definition monomialbasis.hh:782
Definition monomialbasis.hh:785
MonomialBasisProvider< dd, FF > Type
Definition monomialbasis.hh:786