5#ifndef DUNE_MONOMIALBASIS_HH
6#define DUNE_MONOMIALBASIS_HH
14#include <dune/geometry/topologyfactory.hh>
16#include <dune/localfunctions/utility/field.hh>
17#include <dune/localfunctions/utility/multiindex.hh>
18#include <dune/localfunctions/utility/tensor.hh>
71 template< GeometryType::Id geometryId >
72 class MonomialBasisSize;
74 template< GeometryType::Id geometryId,
class F >
82 template< GeometryType::Id geometryId >
83 class MonomialBasisSize
85 typedef MonomialBasisSize< geometryId > This;
88 static This &instance ()
90 static This _instance;
94 unsigned int maxOrder_;
97 mutable unsigned int *sizes_;
100 mutable unsigned int *numBaseFunctions_;
105 numBaseFunctions_( 0 )
110 ~MonomialBasisSize ()
113 delete[] numBaseFunctions_;
116 unsigned int operator() (
const unsigned int order )
const
118 return numBaseFunctions_[ order ];
121 unsigned int maxOrder()
const
126 void computeSizes (
unsigned int order )
128 if (order <= maxOrder_)
134 delete[] numBaseFunctions_;
135 sizes_ =
new unsigned int[ order+1 ];
136 numBaseFunctions_ =
new unsigned int[ order+1 ];
139 constexpr auto dim = geometry.dim();
142 for(
unsigned int k = 1; k <= order; ++k )
145 std::fill(numBaseFunctions_, numBaseFunctions_+order+1, 1);
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 )
153 sizes_[ k ] = numBaseFunctions_[ k ] + k*sizes_[ k ];
154 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
159 for(
unsigned int k = 1; k <= order; ++k )
161 sizes_[ k ] = numBaseFunctions_[ k ];
162 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
175 template<
int mydim,
int dim,
class F >
176 struct MonomialBasisHelper
181 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
182 const unsigned int numBaseFunctions,
const F &z )
185 MySize &mySize = MySize::instance();
186 Size &
size = Size::instance();
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);
217 rit +=
size.sizes_[d] - mySize.sizes_[d];
218 prit +=
size.sizes_[d-1] - mySize.sizes_[d-1];
219 const F *
const emptyWitEnd = wit +
size.sizes_[d] - mySize.sizes_[d];
220 for ( ; wit != emptyWitEnd; ++wit )
232 template< GeometryType::Id geometryId,
class F>
233 class MonomialBasisImpl
240 static const unsigned int dimDomain = geometry.dim();
242 typedef FieldVector< Field, dimDomain > DomainVector;
245 friend class MonomialBasis< geometryId, Field >;
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 *values = Unity< F >();
258 F *
const end = values + block;
259 for( Field *it = values+1 ; it != end; ++it )
264 if constexpr ( geometry ==
gt)
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)
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 >();
321 if constexpr ( geometry ==
gt)
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)
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;
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 )
384 *it = (Field( j ) / Field( j+1 )) * (*row1);
386 for( ; it != row1End; ++row0, ++it )
387 *it = (Field( k ) / Field( k+1 )) * (*row0);
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 )
413 const Field factor = (Field( 1 ) / Field( k + baseGeometry.dim()+1));
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 >
439 :
public MonomialBasisImpl< geometryId, F >
442 typedef MonomialBasis< geometryId, F > This;
443 typedef MonomialBasisImpl< geometryId, F > Base;
446 static const unsigned int dimension = Base::dimDomain;
447 static const unsigned int dimRange = 1;
449 typedef typename Base::Field Field;
451 typedef typename Base::DomainVector DomainVector;
455 typedef MonomialBasisSize<geometryId> Size;
457 MonomialBasis (
unsigned int order)
460 size_(Size::instance())
465 const unsigned int *sizes (
unsigned int order )
const
467 size_.computeSizes( order );
468 return size_.numBaseFunctions_;
471 const unsigned int *sizes ()
const
473 return sizes( order_ );
476 unsigned int size ()
const
478 size_.computeSizes( order_ );
479 return size_( order_ );
482 unsigned int derivSize (
const unsigned int deriv )
const
488 unsigned int order ()
const
493 unsigned int topologyId ( )
const
495 return geometry.id();
498 void evaluate (
const unsigned int deriv,
const DomainVector &x,
499 Field *
const values )
const
501 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
504 template <
unsigned int deriv>
505 void evaluate (
const DomainVector &x,
506 Field *
const values )
const
508 evaluate( deriv, x, values );
511 template<
unsigned int deriv,
class Vector >
512 void evaluate (
const DomainVector &x,
513 Vector &values )
const
515 evaluate<deriv>(x,&(values[0]));
517 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
518 void evaluate (
const DomainVector &x,
519 Derivatives<Field,dimension,1,deriv,layout> *values )
const
521 evaluate<deriv>(x,&(values->block()));
523 template<
unsigned int deriv >
524 void evaluate (
const DomainVector &x,
525 FieldVector<Field,Derivatives<Field,dimension,1,deriv,DerivativeLayoutNS::value>::size> *values )
const
527 evaluate(0,x,&(values[0][0]));
530 template<
class Vector >
531 void evaluate (
const DomainVector &x,
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);
542 for(
int d = 0; d < dimension; ++d )
544 evaluate<0>( bx, values );
547 void integrate ( Field *
const values )
const
549 Base::integrate( order_, sizes( order_ ), values );
551 template <
class Vector>
552 void integrate ( Vector &values )
const
554 integrate( &(values[ 0 ]) );
557 MonomialBasis(
const This&);
558 This& operator=(
const This&);
568 template<
int dim,
class F >
569 class StandardMonomialBasis
570 :
public MonomialBasis< GeometryTypes::simplex(dim).toId() , F >
572 typedef StandardMonomialBasis< dim, F > This;
577 static const int dimension = dim;
579 StandardMonomialBasis (
unsigned int order )
589 template<
int dim,
class F >
590 class StandardBiMonomialBasis
591 :
public MonomialBasis< GeometryTypes::cube(dim).toId() , F >
593 typedef StandardBiMonomialBasis< dim, F > This;
598 static const int dimension = dim;
600 StandardBiMonomialBasis (
unsigned int order )
610 template<
int dim,
class F >
611 class VirtualMonomialBasis
613 typedef VirtualMonomialBasis< dim, F > This;
617 typedef F StorageField;
618 static const int dimension = dim;
619 static const unsigned int dimRange = 1;
621 typedef FieldVector<Field,dimension> DomainVector;
622 typedef FieldVector<Field,dimRange> RangeVector;
624 explicit VirtualMonomialBasis(
const GeometryType&
gt,
626 : order_(order), geometry_(
gt) {}
628 virtual ~VirtualMonomialBasis() {}
630 virtual const unsigned int *sizes ( )
const = 0;
632 unsigned int size ( )
const
634 return sizes( )[ order_ ];
637 unsigned int order ()
const
647 virtual void evaluate (
const unsigned int deriv,
const DomainVector &x,
648 Field *
const values )
const = 0;
649 template <
unsigned int deriv >
650 void evaluate (
const DomainVector &x,
651 Field *
const values )
const
653 evaluate( deriv, x, values );
655 template <
unsigned int deriv,
int size >
656 void evaluate (
const DomainVector &x,
659 evaluate( deriv, x, &(values[0][0]) );
661 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
662 void evaluate (
const DomainVector &x,
663 Derivatives<Field,dimension,1,deriv,layout> *values )
const
665 evaluate<deriv>(x,&(values->block()));
667 template <
unsigned int deriv,
class Vector>
668 void evaluate (
const DomainVector &x,
669 Vector &values )
const
671 evaluate<deriv>( x, &(values[ 0 ]) );
673 template<
class Vector >
674 void evaluate (
const DomainVector &x,
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);
684 for(
int d = 0; d < dimension; ++d )
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);
693 for(
int d = 0; d < dimension; ++d )
695 evaluate<deriv>( bx, values );
698 virtual void integrate ( Field *
const values )
const = 0;
699 template <
class Vector>
700 void integrate ( Vector &values )
const
702 integrate( &(values[ 0 ]) );
709 template< GeometryType::Id geometryId,
class F >
710 class VirtualMonomialBasisImpl
711 :
public VirtualMonomialBasis< static_cast<GeometryType>(geometryId).dim(), F >
714 typedef VirtualMonomialBasis< geometry.dim(), F > Base;
715 typedef VirtualMonomialBasisImpl< geometryId, F > This;
718 typedef typename Base::Field Field;
719 typedef typename Base::DomainVector DomainVector;
721 VirtualMonomialBasisImpl(
unsigned int order)
722 : Base(geometry,order), basis_(order)
725 const unsigned int *sizes ( )
const
727 return basis_.sizes(order_);
730 void evaluate (
const unsigned int deriv,
const DomainVector &x,
731 Field *
const values )
const
733 basis_.evaluate(deriv,x,values);
736 void integrate ( Field *
const values )
const
738 basis_.integrate(values);
742 MonomialBasis<geometryId,Field> basis_;
749 template<
int dim,
class F >
750 struct MonomialBasisFactory
752 static const unsigned int dimension = dim;
753 typedef F StorageField;
755 typedef unsigned int Key;
756 typedef const VirtualMonomialBasis< dimension, F > Object;
758 template <
int dd,
class FF >
759 struct EvaluationBasisFactory
761 typedef MonomialBasisFactory<dd,FF> Type;
764 template< GeometryType::Id geometryId >
765 static Object* create (
const Key &order )
767 return new VirtualMonomialBasisImpl< geometryId, StorageField >( order );
769 static void release( Object *
object ) {
delete object; }
777 template<
int dim,
class SF >
778 struct MonomialBasisProvider
779 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
781 static const unsigned int dimension = dim;
782 typedef SF StorageField;
783 template <
int dd,
class FF >
784 struct EvaluationBasisFactory
786 typedef MonomialBasisProvider<dd,FF> Type;
vector space out of a tensor product of fields.
Definition: fvector.hh:92
constexpr Id toId() const
Create an Id representation of this GeometryType.
Definition: type.hh:210
GeometryType
Type representing VTK's entity geometry types.
Definition: common.hh:132
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:158
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr GeometryType simplex(unsigned int dim)
Returns a GeometryType representing a simplex of dimension dim.
Definition: type.hh:453
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
constexpr GeometryType prismaticExtension(const GeometryType >)
Return GeometryType of a prismatic construction with gt as base
Definition: type.hh:483
constexpr GeometryType conicalExtension(const GeometryType >)
Return GeometryType of a conical construction with gt as base
Definition: type.hh:477
Dune namespace.
Definition: alignedallocator.hh:13
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:159
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
A unique label for each type of element that can occur in a grid.