1 #ifndef DUNE_MONOMIALBASIS_HH
2 #define DUNE_MONOMIALBASIS_HH
6 #include <dune/common/fvector.hh>
7 #include <dune/common/fmatrix.hh>
9 #include <dune/geometry/topologyfactory.hh>
10 #include <dune/geometry/genericgeometry/topologytypes.hh>
67 template<
class Topology >
70 template<
class Topology,
class F >
86 static This _instance;
95 mutable unsigned int maxOrder_;
97 mutable unsigned int *sizes_;
99 mutable unsigned int *numBaseFunctions_;
104 numBaseFunctions_( 0 )
112 delete[] numBaseFunctions_;
115 unsigned int operator() (
const unsigned int order )
const
117 return numBaseFunctions_[ order ];
120 unsigned int maxOrder ()
const
125 void computeSizes (
unsigned int order )
const
127 if (order <= maxOrder_)
133 delete [] numBaseFunctions_;
134 sizes_ =
new unsigned int [ order+1 ];
135 numBaseFunctions_ =
new unsigned int [ order+1 ];
138 numBaseFunctions_[ 0 ] = 1;
139 for(
unsigned int k = 1; k <= order; ++k )
142 numBaseFunctions_[ k ] = 1;
147 template<
class BaseTopology >
155 static This _instance;
159 typedef GenericGeometry::Prism< BaseTopology >
Topology;
164 mutable unsigned int maxOrder_;
166 mutable unsigned int *sizes_;
168 mutable unsigned int *numBaseFunctions_;
173 numBaseFunctions_( 0 )
181 delete[] numBaseFunctions_;
184 unsigned int operator() (
const unsigned int order )
const
186 return numBaseFunctions_[ order ];
189 unsigned int maxOrder()
const
194 void computeSizes (
unsigned int order )
const
196 if (order <= maxOrder_)
202 delete[] numBaseFunctions_;
203 sizes_ =
new unsigned int[ order+1 ];
204 numBaseFunctions_ =
new unsigned int[ order+1 ];
208 baseBasis.computeSizes( order );
209 const unsigned int *
const baseSizes = baseBasis.sizes_;
210 const unsigned int *
const baseNBF = baseBasis.numBaseFunctions_;
213 numBaseFunctions_[ 0 ] = 1;
214 for(
unsigned int k = 1; k <= order; ++k )
216 sizes_[ k ] = baseNBF[ k ] + k*baseSizes[ k ];
217 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
222 template<
class BaseTopology >
230 static This _instance;
234 typedef GenericGeometry::Pyramid< BaseTopology >
Topology;
239 mutable unsigned int maxOrder_;
241 mutable unsigned int *sizes_;
243 mutable unsigned int *numBaseFunctions_;
248 numBaseFunctions_( 0 )
256 delete[] numBaseFunctions_;
259 unsigned int operator() (
const unsigned int order )
const
261 return numBaseFunctions_[ order ];
264 unsigned int maxOrder()
const
269 void computeSizes (
unsigned int order )
const
271 if (order <= maxOrder_)
277 delete[] numBaseFunctions_;
278 sizes_ =
new unsigned int[ order+1 ];
279 numBaseFunctions_ =
new unsigned int[ order+1 ];
284 baseBasis.computeSizes( order );
286 const unsigned int *
const baseNBF = baseBasis.numBaseFunctions_;
288 numBaseFunctions_[ 0 ] = 1;
289 for(
unsigned int k = 1; k <= order; ++k )
291 sizes_[ k ] = baseNBF[ k ];
292 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
303 template<
int mydim,
int dim,
class F >
309 static void copy (
const unsigned int deriv, F *&wit, F *&rit,
310 const unsigned int numBaseFunctions,
const F &z )
313 MySize &mySize = MySize::instance();
314 Size &size = Size::instance();
316 const F *
const rend = rit + size( deriv )*numBaseFunctions;
317 for( ; rit != rend; )
324 for(
unsigned d = 1; d <= deriv; ++d )
327 const F *
const derivEnd = rit + mySize.sizes_[ d ];
329 const F *
const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
330 for( ; rit != drend ; ++rit, ++wit )
332 for (
unsigned int j=1;j<d;++j)
334 const F *
const drend = rit + mySize.sizes_[ d-j ] - mySize.sizes_[ d-j-1 ];
335 for( ; rit != drend ; ++prit, ++rit, ++wit )
336 *wit = F(j) * *prit + z * *rit;
338 *wit = F(d) * *prit + z * *rit;
339 ++prit, ++rit, ++wit;
340 assert(derivEnd == rit);
341 rit += size.sizes_[d] - mySize.sizes_[d];
342 prit += size.sizes_[d-1] - mySize.sizes_[d-1];
343 const F *
const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
344 for ( ; wit != emptyWitEnd; ++wit )
356 template<
class Topology,
class F >
368 static const unsigned int dimDomain = Topology::dimension;
378 void evaluate ( const unsigned int deriv, const unsigned int order,
379 const FieldVector< Field, dimD > &x,
380 const unsigned int block, const unsigned int *const offsets,
381 Field *const values ) const
384 F *
const end = values + block;
385 for(
Field *it = values+1 ; it != end; ++it )
389 void integrate (
const unsigned int order,
390 const unsigned int *
const offsets,
391 Field *
const values )
const
397 template<
class BaseTopology,
class F >
403 typedef GenericGeometry::Prism< BaseTopology >
Topology;
406 static const unsigned int dimDomain = Topology::dimension;
424 void evaluate (
const unsigned int deriv,
const unsigned int order,
425 const FieldVector< Field, dimD > &x,
426 const unsigned int block,
const unsigned int *
const offsets,
427 Field *
const values )
const
430 const BaseSize &size = BaseSize::instance();
432 const Field &z = x[ dimDomain-1 ];
435 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
437 Field *row0 = values;
438 for(
unsigned int k = 1; k <= order; ++k )
440 Field *row1 = values + block*offsets[ k-1 ];
441 Field *wit = row1 + block*size.sizes_[ k ];
442 Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
443 Helper::copy( deriv, wit, row0, size( k-1 ), z );
448 void integrate (
const unsigned int order,
449 const unsigned int *
const offsets,
450 Field *
const values )
const
452 const BaseSize &size = BaseSize::instance();
453 const Size &mySize = Size::instance();
455 baseBasis_.integrate( order, offsets, values );
456 const unsigned int *
const baseSizes = size.sizes_;
458 Field *row0 = values;
459 for(
unsigned int k = 1; k <= order; ++k )
461 Field *
const row1begin = values + offsets[ k-1 ];
462 Field *
const row1End = row1begin + mySize.sizes_[ k ];
463 assert( (
unsigned int)(row1End - values) <= offsets[ k ] );
465 Field *row1 = row1begin;
466 Field *it = row1begin + baseSizes[ k ];
467 for(
unsigned int j = 1; j <= k; ++j )
469 Field *
const end = it + baseSizes[ k ];
470 assert( (
unsigned int)(end - values) <= offsets[ k ] );
471 for( ; it != end; ++row1, ++it )
472 *it = (Field( j ) / Field( j+1 )) * (*row1);
474 for( ; it != row1End; ++row0, ++it )
475 *it = (Field( k ) / Field( k+1 )) * (*row0);
481 template<
class BaseTopology,
class F >
487 typedef GenericGeometry::Pyramid< BaseTopology >
Topology;
490 static const unsigned int dimDomain = Topology::dimension;
509 void evaluateSimplexBase (
const unsigned int deriv,
const unsigned int order,
510 const FieldVector< Field, dimD > &x,
511 const unsigned int block,
const unsigned int *
const offsets,
513 const BaseSize &size )
const
515 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
519 void evaluatePyramidBase (
const unsigned int deriv,
const unsigned int order,
520 const FieldVector< Field, dimD > &x,
521 const unsigned int block,
const unsigned int *
const offsets,
523 const BaseSize &size )
const
525 Field omz = Unity< Field >() - x[ dimDomain-1 ];
527 if( Zero< Field >() < omz )
529 const Field invomz = Unity< Field >() / omz;
530 FieldVector< Field, dimD > y;
531 for(
unsigned int i = 0; i < dimDomain-1; ++i )
532 y[ i ] = x[ i ] * invomz;
535 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
538 for(
unsigned int k = 1; k <= order; ++k )
540 Field *it = values + block*offsets[ k-1 ];
541 Field *
const end = it + block*size.sizes_[ k ];
542 for( ; it != end; ++it )
550 *values = Unity< Field >();
551 for(
unsigned int k = 1; k <= order; ++k )
553 Field *it = values + block*offsets[ k-1 ];
554 Field *
const end = it + block*size.sizes_[ k ];
555 for( ; it != end; ++it )
556 *it = Zero< Field >();
562 void evaluate (
const unsigned int deriv,
const unsigned int order,
563 const FieldVector< Field, dimD > &x,
564 const unsigned int block,
const unsigned int *
const offsets,
565 Field *
const values )
const
567 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
568 const BaseSize &size = BaseSize::instance();
571 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
573 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
575 Field *row0 = values;
576 for(
unsigned int k = 1; k <= order; ++k )
578 Field *row1 = values + block*offsets[ k-1 ];
579 Field *wit = row1 + block*size.sizes_[ k ];
580 Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
585 void integrate (
const unsigned int order,
586 const unsigned int *
const offsets,
587 Field *
const values )
const
589 const BaseSize &size = BaseSize::instance();
592 baseBasis_.integrate( order, offsets, values );
594 const unsigned int *
const baseSizes = size.sizes_;
596 Field *
const col0End = values + baseSizes[ 0 ];
597 for( Field *it = values; it != col0End; ++it )
598 *it *= Field( 1 ) / Field(
int(dimDomain) );
599 Field *row0 = values;
601 for(
unsigned int k = 1; k <= order; ++k )
603 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
605 Field *
const row1 = values+offsets[ k-1 ];
606 Field *
const col0End = row1 + baseSizes[ k ];
608 for( ; it != col0End; ++it )
610 for(
unsigned int i = 1; i <= k; ++i )
612 Field *
const end = it + baseSizes[ k-i ];
613 assert( (
unsigned int)(end - values) <= offsets[ k ] );
614 for( ; it != end; ++row0, ++it )
615 *it = (*row0) * (Field( i ) * factor);
627 template<
class Topology,
class F >
629 :
public MonomialBasisImpl< Topology, F >
631 typedef MonomialBasis< Topology, F > This;
632 typedef MonomialBasisImpl< Topology, F > Base;
649 size_(
Size::instance())
654 const unsigned int *
sizes (
unsigned int order )
const
656 size_.computeSizes( order );
657 return size_.numBaseFunctions_;
662 return sizes( order_ );
665 const unsigned int size ()
const
667 size_.computeSizes( order_ );
668 return size_( order_ );
671 const unsigned int derivSize (
const unsigned int deriv )
const
673 typedef typename GenericGeometry::SimplexTopology< dimension >::type SimplexTopology;
689 Field *
const values )
const
691 Base::evaluate( deriv, order_, x,
derivSize( deriv ),
sizes( order_ ), values );
694 template <
unsigned int deriv>
696 Field *
const values )
const
701 template<
unsigned int deriv,
class Vector >
703 Vector &values )
const
705 evaluate<deriv>(x,&(values[0]));
707 template<
unsigned int deriv, DerivativeLayout layout >
711 evaluate<deriv>(x,&(values->block()));
713 template<
unsigned int deriv >
720 template<
class Vector >
722 Vector &values )
const
724 evaluate<0>(x,&(values[0]));
727 template<
class DVector,
class RVector >
728 void evaluate (
const DVector &x, RVector &values )
const
730 assert( DVector::dimension ==
dimension);
734 evaluate<0>( bx, values );
739 Base::integrate( order_,
sizes( order_ ), values );
741 template <
class Vector>
748 This& operator=(
const This&);
758 template<
int dim,
class F >
760 :
public MonomialBasis< typename GenericGeometry::SimplexTopology< dim >::type, F >
766 typedef typename GenericGeometry::SimplexTopology< dim >::type
Topology;
779 template<
int dim,
class F >
781 :
public MonomialBasis< typename GenericGeometry::CubeTopology< dim >::type, F >
787 typedef typename GenericGeometry::CubeTopology< dim >::type
Topology;
800 template<
int dim,
class F >
820 virtual const unsigned int *
sizes ( )
const = 0;
822 const unsigned int size ( )
const
838 Field *
const values )
const = 0;
839 template <
unsigned int deriv >
841 Field *
const values )
const
845 template <
unsigned int deriv,
int size >
847 Dune::FieldVector<Field,size> *
const values )
const
849 evaluate( deriv, x, &(values[0][0]) );
851 template<
unsigned int deriv, DerivativeLayout layout >
855 evaluate<deriv>(x,&(values->block()));
857 template <
unsigned int deriv,
class Vector>
859 Vector &values )
const
861 evaluate<deriv>( x, &(values[ 0 ]) );
863 template<
class Vector >
865 Vector &values )
const
867 evaluate<0>(x,values);
869 template<
class DVector,
class RVector >
870 void evaluate (
const DVector &x, RVector &values )
const
872 assert( DVector::dimension ==
dimension);
876 evaluate<0>( bx, values );
878 template<
unsigned int deriv,
class DVector,
class RVector >
879 void evaluate (
const DVector &x, RVector &values )
const
881 assert( DVector::dimension ==
dimension);
885 evaluate<deriv>( bx, values );
889 template <
class Vector>
899 template<
class Topology,
class F >
911 :
Base(Topology::id,order), basis_(order)
920 Field *
const values )
const
938 template<
int dim,
class F >
941 template<
int dim,
class F >
950 template<
int dim,
class F >
952 public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
961 template <
int dd,
class FF >
967 template<
class Topology >
979 template<
int dim,
class SF >
981 :
public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
985 template <
int dd,
class FF >