3#ifndef DUNE_MONOMIALBASIS_HH
4#define DUNE_MONOMIALBASIS_HH
12#include <dune/geometry/topologyfactory.hh>
14#include <dune/localfunctions/utility/field.hh>
15#include <dune/localfunctions/utility/multiindex.hh>
16#include <dune/localfunctions/utility/tensor.hh>
69 template< GeometryType::Id geometryId >
70 class MonomialBasisSize;
72 template< GeometryType::Id geometryId,
class F >
80 template< GeometryType::Id geometryId >
81 class MonomialBasisSize
83 typedef MonomialBasisSize< geometryId > This;
86 static This &instance ()
88 static This _instance;
92 unsigned int maxOrder_;
95 mutable unsigned int *sizes_;
98 mutable unsigned int *numBaseFunctions_;
103 numBaseFunctions_( 0 )
108 ~MonomialBasisSize ()
111 delete[] numBaseFunctions_;
114 unsigned int operator() (
const unsigned int order )
const
116 return numBaseFunctions_[ order ];
119 unsigned int maxOrder()
const
124 void computeSizes (
unsigned int order )
126 if (order <= maxOrder_)
132 delete[] numBaseFunctions_;
133 sizes_ =
new unsigned int[ order+1 ];
134 numBaseFunctions_ =
new unsigned int[ order+1 ];
137 constexpr auto dim = geometry.dim();
140 for(
unsigned int k = 1; k <= order; ++k )
143 std::fill(numBaseFunctions_, numBaseFunctions_+order+1, 1);
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 )
151 sizes_[ k ] = numBaseFunctions_[ k ] + k*sizes_[ k ];
152 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
157 for(
unsigned int k = 1; k <= order; ++k )
159 sizes_[ k ] = numBaseFunctions_[ k ];
160 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
173 template<
int mydim,
int dim,
class F >
174 struct MonomialBasisHelper
179 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
180 const unsigned int numBaseFunctions,
const F &z )
183 MySize &mySize = MySize::instance();
184 Size &size = Size::instance();
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);
215 rit += size.sizes_[d] - mySize.sizes_[d];
216 prit += size.sizes_[d-1] - mySize.sizes_[d-1];
217 const F *
const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
218 for ( ; wit != emptyWitEnd; ++wit )
230 template< GeometryType::Id geometryId,
class F>
231 class MonomialBasisImpl
238 static const unsigned int dimDomain = geometry.dim();
240 typedef FieldVector< Field, dimDomain > DomainVector;
243 friend class MonomialBasis< geometryId, Field >;
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
255 *values = Unity< F >();
256 F *
const end = values + block;
257 for( Field *it = values+1 ; it != end; ++it )
262 if constexpr ( geometry ==
gt)
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)
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 >();
319 if constexpr ( geometry ==
gt)
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)
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;
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 )
382 *it = (Field( j ) / Field( j+1 )) * (*row1);
384 for( ; it != row1End; ++row0, ++it )
385 *it = (Field( k ) / Field( k+1 )) * (*row0);
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 )
411 const Field factor = (Field( 1 ) / Field( k + baseGeometry.dim()+1));
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 >
437 :
public MonomialBasisImpl< geometryId, F >
440 typedef MonomialBasis< geometryId, F > This;
441 typedef MonomialBasisImpl< geometryId, F > Base;
444 static const unsigned int dimension = Base::dimDomain;
445 static const unsigned int dimRange = 1;
447 typedef typename Base::Field Field;
449 typedef typename Base::DomainVector DomainVector;
453 typedef MonomialBasisSize<geometryId> Size;
455 MonomialBasis (
unsigned int order)
458 size_(Size::instance())
463 const unsigned int *sizes (
unsigned int order )
const
465 size_.computeSizes( order );
466 return size_.numBaseFunctions_;
469 const unsigned int *sizes ()
const
471 return sizes( order_ );
474 unsigned int size ()
const
476 size_.computeSizes( order_ );
477 return size_( order_ );
480 unsigned int derivSize (
const unsigned int deriv )
const
486 unsigned int order ()
const
491 unsigned int topologyId ( )
const
493 return geometry.id();
496 void evaluate (
const unsigned int deriv,
const DomainVector &x,
497 Field *
const values )
const
499 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
502 template <
unsigned int deriv>
503 void evaluate (
const DomainVector &x,
504 Field *
const values )
const
506 evaluate( deriv, x, values );
509 template<
unsigned int deriv,
class Vector >
510 void evaluate (
const DomainVector &x,
511 Vector &values )
const
513 evaluate<deriv>(x,&(values[0]));
515 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
516 void evaluate (
const DomainVector &x,
517 Derivatives<Field,dimension,1,deriv,layout> *values )
const
519 evaluate<deriv>(x,&(values->block()));
521 template<
unsigned int deriv >
522 void evaluate (
const DomainVector &x,
523 FieldVector<Field,Derivatives<Field,dimension,1,deriv,DerivativeLayoutNS::value>::size> *values )
const
525 evaluate(0,x,&(values[0][0]));
528 template<
class Vector >
529 void evaluate (
const DomainVector &x,
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);
540 for(
int d = 0; d < dimension; ++d )
542 evaluate<0>( bx, values );
545 void integrate ( Field *
const values )
const
547 Base::integrate( order_, sizes( order_ ), values );
549 template <
class Vector>
550 void integrate ( Vector &values )
const
552 integrate( &(values[ 0 ]) );
555 MonomialBasis(
const This&);
556 This& operator=(
const This&);
566 template<
int dim,
class F >
567 class StandardMonomialBasis
568 :
public MonomialBasis< GeometryTypes::simplex(dim).toId() , F >
570 typedef StandardMonomialBasis< dim, F > This;
575 static const int dimension = dim;
577 StandardMonomialBasis (
unsigned int order )
587 template<
int dim,
class F >
588 class StandardBiMonomialBasis
589 :
public MonomialBasis< GeometryTypes::cube(dim).toId() , F >
591 typedef StandardBiMonomialBasis< dim, F > This;
596 static const int dimension = dim;
598 StandardBiMonomialBasis (
unsigned int order )
608 template<
int dim,
class F >
609 class VirtualMonomialBasis
611 typedef VirtualMonomialBasis< dim, F > This;
615 typedef F StorageField;
616 static const int dimension = dim;
617 static const unsigned int dimRange = 1;
619 typedef FieldVector<Field,dimension> DomainVector;
620 typedef FieldVector<Field,dimRange> RangeVector;
622 [[deprecated(
"Use VirtualMonomialBasis(GeometryType gt, unsigned int order) instead.")]]
623 explicit VirtualMonomialBasis(
unsigned int topologyId,
625 : order_(order), geometry_(
GeometryType(topologyId,dim)) {}
627 explicit VirtualMonomialBasis(
const GeometryType&
gt,
629 : order_(order), geometry_(
gt) {}
631 virtual ~VirtualMonomialBasis() {}
633 virtual const unsigned int *sizes ( )
const = 0;
635 unsigned int size ( )
const
637 return sizes( )[ order_ ];
640 unsigned int order ()
const
645 [[deprecated(
"Use type().id() instead.")]]
646 unsigned int topologyId ( )
const
656 virtual void evaluate (
const unsigned int deriv,
const DomainVector &x,
657 Field *
const values )
const = 0;
658 template <
unsigned int deriv >
659 void evaluate (
const DomainVector &x,
660 Field *
const values )
const
662 evaluate( deriv, x, values );
664 template <
unsigned int deriv,
int size >
665 void evaluate (
const DomainVector &x,
668 evaluate( deriv, x, &(values[0][0]) );
670 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
671 void evaluate (
const DomainVector &x,
672 Derivatives<Field,dimension,1,deriv,layout> *values )
const
674 evaluate<deriv>(x,&(values->block()));
676 template <
unsigned int deriv,
class Vector>
677 void evaluate (
const DomainVector &x,
678 Vector &values )
const
680 evaluate<deriv>( x, &(values[ 0 ]) );
682 template<
class Vector >
683 void evaluate (
const DomainVector &x,
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);
693 for(
int d = 0; d < dimension; ++d )
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);
702 for(
int d = 0; d < dimension; ++d )
704 evaluate<deriv>( bx, values );
707 virtual void integrate ( Field *
const values )
const = 0;
708 template <
class Vector>
709 void integrate ( Vector &values )
const
711 integrate( &(values[ 0 ]) );
718 template< GeometryType::Id geometryId,
class F >
719 class VirtualMonomialBasisImpl
720 :
public VirtualMonomialBasis< static_cast<GeometryType>(geometryId).dim(), F >
723 typedef VirtualMonomialBasis< geometry.dim(), F > Base;
724 typedef VirtualMonomialBasisImpl< geometryId, F > This;
727 typedef typename Base::Field Field;
728 typedef typename Base::DomainVector DomainVector;
730 VirtualMonomialBasisImpl(
unsigned int order)
731 : Base(geometry,order), basis_(order)
734 const unsigned int *sizes ( )
const
736 return basis_.sizes(order_);
739 void evaluate (
const unsigned int deriv,
const DomainVector &x,
740 Field *
const values )
const
742 basis_.evaluate(deriv,x,values);
745 void integrate ( Field *
const values )
const
747 basis_.integrate(values);
751 MonomialBasis<geometryId,Field> basis_;
758 template<
int dim,
class F >
759 struct MonomialBasisFactory
761 static const unsigned int dimension = dim;
762 typedef F StorageField;
764 typedef unsigned int Key;
765 typedef const VirtualMonomialBasis< dimension, F > Object;
767 template <
int dd,
class FF >
768 struct EvaluationBasisFactory
770 typedef MonomialBasisFactory<dd,FF> Type;
773 template< GeometryType::Id geometryId >
774 static Object* create (
const Key &order )
776 return new VirtualMonomialBasisImpl< geometryId, StorageField >( order );
778 static void release( Object *
object ) {
delete object; }
786 template<
int dim,
class SF >
787 struct MonomialBasisProvider
788 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
790 static const unsigned int dimension = dim;
791 typedef SF StorageField;
792 template <
int dd,
class FF >
793 struct EvaluationBasisFactory
795 typedef MonomialBasisProvider<dd,FF> Type;
vector space out of a tensor product of fields.
Definition: fvector.hh:95
constexpr Id toId() const
Create an Id representation of this GeometryType.
Definition: type.hh:219
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:130
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:470
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:461
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:504
constexpr GeometryType prismaticExtension(const GeometryType >)
Return GeometryType of a prismatic construction with gt as base
Definition: type.hh:491
constexpr GeometryType conicalExtension(const GeometryType >)
Return GeometryType of a conical construction with gt as base
Definition: type.hh:485
Dune namespace.
Definition: alignedallocator.hh:11
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
A unique label for each type of element that can occur in a grid.