3#ifndef DUNE_MONOMIALBASIS_HH
4#define DUNE_MONOMIALBASIS_HH
11#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<
class Topology >
70 class MonomialBasisSize;
72 template<
class Topology,
class F >
80 template<
class TopologyType >
81 class MonomialBasisSize
83 typedef MonomialBasisSize< TopologyType > 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 ];
136 constexpr auto dim = TopologyType::dimension;
139 for(
unsigned int k = 1; k <= order; ++k )
142 std::fill(numBaseFunctions_, numBaseFunctions_+order+1, 1);
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 )
150 sizes_[ k ] = numBaseFunctions_[ k ] + k*sizes_[ k ];
151 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
156 for(
unsigned int k = 1; k <= order; ++k )
158 sizes_[ k ] = numBaseFunctions_[ k ];
159 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
172 template<
int mydim,
int dim,
class F >
173 struct MonomialBasisHelper
175 typedef MonomialBasisSize< typename Impl::SimplexTopology< mydim >::type > MySize;
176 typedef MonomialBasisSize< typename Impl::SimplexTopology< dim >::type > Size;
178 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
179 const unsigned int numBaseFunctions,
const F &z )
182 MySize &mySize = MySize::instance();
183 Size &size = Size::instance();
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);
214 rit += size.sizes_[d] - mySize.sizes_[d];
215 prit += size.sizes_[d-1] - mySize.sizes_[d-1];
216 const F *
const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
217 for ( ; wit != emptyWitEnd; ++wit )
229 template<
class Topology,
class F >
230 class MonomialBasisImpl;
233 class MonomialBasisImpl< Impl::Point, F >
235 typedef MonomialBasisImpl< Impl::Point, F > This;
238 typedef Impl::Point Topology;
241 static const unsigned int dimDomain = Topology::dimension;
243 typedef FieldVector< Field, dimDomain > DomainVector;
246 friend class MonomialBasis< Topology, Field >;
247 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
248 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, 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
256 *values = Unity< F >();
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
266 values[ 0 ] = Unity< Field >();
270 template<
class BaseTopology,
class F >
271 class MonomialBasisImpl< Impl::Prism< BaseTopology >, F >
273 typedef MonomialBasisImpl< Impl::Prism< BaseTopology >, F > This;
276 typedef Impl::Prism< BaseTopology > Topology;
279 static const unsigned int dimDomain = Topology::dimension;
281 typedef FieldVector< Field, dimDomain > DomainVector;
284 friend class MonomialBasis< Topology, Field >;
285 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
286 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
288 typedef MonomialBasisSize< BaseTopology > BaseSize;
289 typedef MonomialBasisSize< Topology > Size;
291 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
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
302 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
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 >
356 class MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F >
358 typedef MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F > This;
361 typedef Impl::Pyramid< BaseTopology > Topology;
364 static const unsigned int dimDomain = Topology::dimension;
366 typedef FieldVector< Field, dimDomain > DomainVector;
369 friend class MonomialBasis< Topology, Field >;
370 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
371 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
373 typedef MonomialBasisSize< BaseTopology > BaseSize;
374 typedef MonomialBasisSize< Topology > Size;
376 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
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);
444 if( Impl::IsSimplex< Topology >::value )
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 >
505 :
public MonomialBasisImpl< Topology, F >
507 typedef MonomialBasis< Topology, F > This;
508 typedef MonomialBasisImpl< Topology, F > Base;
511 static const unsigned int dimension = Base::dimDomain;
512 static const unsigned int dimRange = 1;
514 typedef typename Base::Field Field;
516 typedef typename Base::DomainVector DomainVector;
520 typedef MonomialBasisSize<Topology> Size;
522 MonomialBasis (
unsigned int order)
525 size_(Size::instance())
530 const unsigned int *sizes (
unsigned int order )
const
532 size_.computeSizes( order );
533 return size_.numBaseFunctions_;
536 const unsigned int *sizes ()
const
538 return sizes( order_ );
541 unsigned int size ()
const
543 size_.computeSizes( order_ );
544 return size_( order_ );
547 unsigned int derivSize (
const unsigned int deriv )
const
549 typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology;
550 MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
551 return MonomialBasisSize< SimplexTopology >::instance() ( deriv );
554 unsigned int order ()
const
559 unsigned int topologyId ( )
const
564 void evaluate (
const unsigned int deriv,
const DomainVector &x,
565 Field *
const values )
const
567 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
570 template <
unsigned int deriv>
571 void evaluate (
const DomainVector &x,
572 Field *
const values )
const
574 evaluate( deriv, x, values );
577 template<
unsigned int deriv,
class Vector >
578 void evaluate (
const DomainVector &x,
579 Vector &values )
const
581 evaluate<deriv>(x,&(values[0]));
583 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
584 void evaluate (
const DomainVector &x,
585 Derivatives<Field,dimension,1,deriv,layout> *values )
const
587 evaluate<deriv>(x,&(values->block()));
589 template<
unsigned int deriv >
590 void evaluate (
const DomainVector &x,
591 FieldVector<Field,Derivatives<Field,dimension,1,deriv,DerivativeLayoutNS::value>::size> *values )
const
593 evaluate(0,x,&(values[0][0]));
596 template<
class Vector >
597 void evaluate (
const DomainVector &x,
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);
608 for(
int d = 0; d < dimension; ++d )
610 evaluate<0>( bx, values );
613 void integrate ( Field *
const values )
const
615 Base::integrate( order_, sizes( order_ ), values );
617 template <
class Vector>
618 void integrate ( Vector &values )
const
620 integrate( &(values[ 0 ]) );
623 MonomialBasis(
const This&);
624 This& operator=(
const This&);
634 template<
int dim,
class F >
635 class StandardMonomialBasis
636 :
public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F >
638 typedef StandardMonomialBasis< dim, F > This;
639 typedef MonomialBasis< typename Impl::SimplexTopology< dim >::type, F > Base;
642 typedef typename Impl::SimplexTopology< dim >::type Topology;
643 static const int dimension = dim;
645 StandardMonomialBasis (
unsigned int order )
655 template<
int dim,
class F >
656 class StandardBiMonomialBasis
657 :
public MonomialBasis< typename Impl::CubeTopology< dim >::type, F >
659 typedef StandardBiMonomialBasis< dim, F > This;
660 typedef MonomialBasis< typename Impl::CubeTopology< dim >::type, F > Base;
663 typedef typename Impl::CubeTopology< dim >::type Topology;
664 static const int dimension = dim;
666 StandardBiMonomialBasis (
unsigned int order )
676 template<
int dim,
class F >
677 class VirtualMonomialBasis
679 typedef VirtualMonomialBasis< dim, F > This;
683 typedef F StorageField;
684 static const int dimension = dim;
685 static const unsigned int dimRange = 1;
687 typedef FieldVector<Field,dimension> DomainVector;
688 typedef FieldVector<Field,dimRange> RangeVector;
690 explicit VirtualMonomialBasis(
unsigned int topologyId,
692 : order_(order), topologyId_(topologyId) {}
694 virtual ~VirtualMonomialBasis() {}
696 virtual const unsigned int *sizes ( )
const = 0;
698 unsigned int size ( )
const
700 return sizes( )[ order_ ];
703 unsigned int order ()
const
708 unsigned int topologyId ( )
const
713 virtual void evaluate (
const unsigned int deriv,
const DomainVector &x,
714 Field *
const values )
const = 0;
715 template <
unsigned int deriv >
716 void evaluate (
const DomainVector &x,
717 Field *
const values )
const
719 evaluate( deriv, x, values );
721 template <
unsigned int deriv,
int size >
722 void evaluate (
const DomainVector &x,
725 evaluate( deriv, x, &(values[0][0]) );
727 template<
unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
728 void evaluate (
const DomainVector &x,
729 Derivatives<Field,dimension,1,deriv,layout> *values )
const
731 evaluate<deriv>(x,&(values->block()));
733 template <
unsigned int deriv,
class Vector>
734 void evaluate (
const DomainVector &x,
735 Vector &values )
const
737 evaluate<deriv>( x, &(values[ 0 ]) );
739 template<
class Vector >
740 void evaluate (
const DomainVector &x,
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);
750 for(
int d = 0; d < dimension; ++d )
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);
759 for(
int d = 0; d < dimension; ++d )
761 evaluate<deriv>( bx, values );
764 virtual void integrate ( Field *
const values )
const = 0;
765 template <
class Vector>
766 void integrate ( Vector &values )
const
768 integrate( &(values[ 0 ]) );
772 unsigned int topologyId_;
775 template<
class Topology,
class F >
776 class VirtualMonomialBasisImpl
777 :
public VirtualMonomialBasis< Topology::dimension, F >
779 typedef VirtualMonomialBasis< Topology::dimension, F > Base;
780 typedef VirtualMonomialBasisImpl< Topology, F > This;
783 typedef typename Base::Field Field;
784 typedef typename Base::DomainVector DomainVector;
786 VirtualMonomialBasisImpl(
unsigned int order)
787 : Base(Topology::id,order), basis_(order)
790 const unsigned int *sizes ( )
const
792 return basis_.sizes(order_);
795 void evaluate (
const unsigned int deriv,
const DomainVector &x,
796 Field *
const values )
const
798 basis_.evaluate(deriv,x,values);
801 void integrate ( Field *
const values )
const
803 basis_.integrate(values);
807 MonomialBasis<Topology,Field> basis_;
814 template<
int dim,
class F >
815 struct MonomialBasisFactory
817 static const unsigned int dimension = dim;
818 typedef F StorageField;
820 typedef unsigned int Key;
821 typedef const VirtualMonomialBasis< dimension, F > Object;
823 template <
int dd,
class FF >
824 struct EvaluationBasisFactory
826 typedef MonomialBasisFactory<dd,FF> Type;
829 template<
class Topology >
830 static Object* create (
const Key &order )
832 return new VirtualMonomialBasisImpl< Topology, StorageField >( order );
834 static void release( Object *
object ) {
delete object; }
842 template<
int dim,
class SF >
843 struct MonomialBasisProvider
844 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
846 static const unsigned int dimension = dim;
847 typedef SF StorageField;
848 template <
int dd,
class FF >
849 struct EvaluationBasisFactory
851 typedef MonomialBasisProvider<dd,FF> Type;
vector space out of a tensor product of fields.
Definition: fvector.hh:96
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.
Dune namespace.
Definition: alignedallocator.hh:14
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.