DUNE PDELab (2.7)

monomialbasis.hh
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_MONOMIALBASIS_HH
4#define DUNE_MONOMIALBASIS_HH
5
6#include <vector>
7
8#include <dune/common/fvector.hh>
9#include <dune/common/fmatrix.hh>
10
11#include <dune/geometry/topologyfactory.hh>
12#include <dune/geometry/type.hh>
13
14#include <dune/localfunctions/utility/field.hh>
15#include <dune/localfunctions/utility/multiindex.hh>
16#include <dune/localfunctions/utility/tensor.hh>
17
18namespace Dune
19{
20 /************************************************
21 * Classes for evaluating ''Monomials'' on any order
22 * for all reference element type.
23 * For a simplex topology these are the normal
24 * monomials for cube topologies the bimonomials.
25 * The construction follows the construction of the
26 * generic geometries using tensor products for
27 * prism generation and duffy transform for pyramid
28 * construction.
29 * A derivative argument can be applied, in which case
30 * all derivatives up to the desired order are
31 * evaluated. Note that for higher order derivatives
32 * only the ''lower'' part of the symmetric tensor
33 * is evaluated, e.g., passing derivative equal to 2
34 * to the class will provide the vector
35 * (d/dxdx p, d/dxydx p, d/dydy p,
36 * d/dx p, d/dy p, p)
37 * Important:
38 * So far the computation of the derivatives has not
39 * been fully implemented for general pyramid
40 * construction, i.e., in the case where a pyramid is
41 * build over a non simplex base geometry.
42 *
43 * Central classes:
44 * 1) template< class Topology, class F >
45 * class MonomialBasisImpl;
46 * Implementation of the monomial evaluation for
47 * a given topology and field type.
48 * The method evaluate fills a F* vector
49 * 2) template< class Topology, class F >
50 * class MonomialBasis
51 * The base class for the static monomial evaluation
52 * providing addiional evaluate methods including
53 * one taking std::vector<F>.
54 * 3) template< int dim, class F >
55 * class VirtualMonomialBasis
56 * Virtualization of the MonomialBasis.
57 * 4) template< int dim, class F >
58 * struct MonomialBasisFactory;
59 * A factory class for the VirtualMonomialBasis
60 * 5) template< int dim, class F >
61 * struct MonomialBasisProvider
62 * A singleton container for the virtual monomial
63 * basis
64 ************************************************/
65
66 // Internal Forward Declarations
67 // -----------------------------
68
69 template< class Topology >
70 class MonomialBasisSize;
71
72 template< class Topology, class F >
73 class MonomialBasis;
74
75
76
77 // MonomialBasisSize
78 // -----------------
79
80 template< class TopologyType >
81 class MonomialBasisSize
82 {
83 typedef MonomialBasisSize< TopologyType > This;
84
85 public:
86 static This &instance ()
87 {
88 static This _instance;
89 return _instance;
90 }
91
92 unsigned int maxOrder_;
93
94 // sizes_[ k ]: number of basis functions of exactly order k
95 mutable unsigned int *sizes_;
96
97 // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
98 mutable unsigned int *numBaseFunctions_;
99
100 MonomialBasisSize ()
101 : maxOrder_( 0 ),
102 sizes_( 0 ),
103 numBaseFunctions_( 0 )
104 {
105 computeSizes( 2 );
106 }
107
108 ~MonomialBasisSize ()
109 {
110 delete[] sizes_;
111 delete[] numBaseFunctions_;
112 }
113
114 unsigned int operator() ( const unsigned int order ) const
115 {
116 return numBaseFunctions_[ order ];
117 }
118
119 unsigned int maxOrder() const
120 {
121 return maxOrder_;
122 }
123
124 void computeSizes ( unsigned int order )
125 {
126 if (order <= maxOrder_)
127 return;
128
129 maxOrder_ = order;
130
131 delete[] sizes_;
132 delete[] numBaseFunctions_;
133 sizes_ = new unsigned int[ order+1 ];
134 numBaseFunctions_ = new unsigned int[ order+1 ];
135
136 constexpr auto dim = TopologyType::dimension;
137
138 sizes_[ 0 ] = 1;
139 for( unsigned int k = 1; k <= order; ++k )
140 sizes_[ k ] = 0;
141
142 std::fill(numBaseFunctions_, numBaseFunctions_+order+1, 1);
143
144 for( int codim=dim-1; codim>=0; codim--)
145 {
146 if (Impl::isPrism(TopologyType::id,dim,codim))
147 {
148 for( unsigned int k = 1; k <= order; ++k )
149 {
150 sizes_[ k ] = numBaseFunctions_[ k ] + k*sizes_[ k ];
151 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
152 }
153 }
154 else
155 {
156 for( unsigned int k = 1; k <= order; ++k )
157 {
158 sizes_[ k ] = numBaseFunctions_[ k ];
159 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
160 }
161 }
162 }
163 }
164 };
165
166
167
168 // MonomialBasisHelper
169 // -------------------
170
171
172 template< int mydim, int dim, class F >
173 struct MonomialBasisHelper
174 {
175 typedef MonomialBasisSize< typename Impl::SimplexTopology< mydim >::type > MySize;
176 typedef MonomialBasisSize< typename Impl::SimplexTopology< dim >::type > Size;
177
178 static void copy ( const unsigned int deriv, F *&wit, F *&rit,
179 const unsigned int numBaseFunctions, const F &z )
180 {
181 // n(d,k) = size<k>[d];
182 MySize &mySize = MySize::instance();
183 Size &size = Size::instance();
184
185 const F *const rend = rit + size( deriv )*numBaseFunctions;
186 for( ; rit != rend; )
187 {
188 F *prit = rit;
189
190 *wit = z * *rit;
191 ++rit, ++wit;
192
193 for( unsigned d = 1; d <= deriv; ++d )
194 {
195 #ifndef NDEBUG
196 const F *const derivEnd = rit + mySize.sizes_[ d ];
197 #endif
198
199 {
200 const F *const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
201 for( ; rit != drend ; ++rit, ++wit )
202 *wit = z * *rit;
203 }
204
205 for (unsigned int j=1; j<d; ++j)
206 {
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;
210 }
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 )
218 *wit = Zero<F>();
219 }
220 }
221 }
222 };
223
224
225
226 // MonomialBasisImpl
227 // -----------------
228
229 template< class Topology, class F >
230 class MonomialBasisImpl;
231
232 template< class F >
233 class MonomialBasisImpl< Impl::Point, F >
234 {
235 typedef MonomialBasisImpl< Impl::Point, F > This;
236
237 public:
238 typedef Impl::Point Topology;
239 typedef F Field;
240
241 static const unsigned int dimDomain = Topology::dimension;
242
243 typedef FieldVector< Field, dimDomain > DomainVector;
244
245 private:
246 friend class MonomialBasis< Topology, Field >;
247 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
248 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
249
250 template< int dimD >
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
255 {
256 *values = Unity< F >();
257 F *const end = values + block;
258 for( Field *it = values+1 ; it != end; ++it )
259 *it = Zero< F >();
260 }
261
262 void integrate ( const unsigned int order,
263 const unsigned int *const offsets,
264 Field *const values ) const
265 {
266 values[ 0 ] = Unity< Field >();
267 }
268 };
269
270 template< class BaseTopology, class F >
271 class MonomialBasisImpl< Impl::Prism< BaseTopology >, F >
272 {
273 typedef MonomialBasisImpl< Impl::Prism< BaseTopology >, F > This;
274
275 public:
276 typedef Impl::Prism< BaseTopology > Topology;
277 typedef F Field;
278
279 static const unsigned int dimDomain = Topology::dimension;
280
281 typedef FieldVector< Field, dimDomain > DomainVector;
282
283 private:
284 friend class MonomialBasis< Topology, Field >;
285 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
286 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
287
288 typedef MonomialBasisSize< BaseTopology > BaseSize;
289 typedef MonomialBasisSize< Topology > Size;
290
291 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
292
293 MonomialBasisImpl ()
294 {}
295
296 template< int dimD >
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
301 {
302 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
303 const BaseSize &size = BaseSize::instance();
304 const_cast<BaseSize&>(size).computeSizes(order);
305
306 const Field &z = x[ dimDomain-1 ];
307
308 // fill first column
309 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
310
311 Field *row0 = values;
312 for( unsigned int k = 1; k <= order; ++k )
313 {
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 );
318 row0 = row1;
319 }
320 }
321
322 void integrate ( const unsigned int order,
323 const unsigned int *const offsets,
324 Field *const values ) const
325 {
326 const BaseSize &size = BaseSize::instance();
327 const Size &mySize = Size::instance();
328 // fill first column
329 baseBasis_.integrate( order, offsets, values );
330 const unsigned int *const baseSizes = size.sizes_;
331
332 Field *row0 = values;
333 for( unsigned int k = 1; k <= order; ++k )
334 {
335 Field *const row1begin = values + offsets[ k-1 ];
336 Field *const row1End = row1begin + mySize.sizes_[ k ];
337 assert( (unsigned int)(row1End - values) <= offsets[ k ] );
338
339 Field *row1 = row1begin;
340 Field *it = row1begin + baseSizes[ k ];
341 for( unsigned int j = 1; j <= k; ++j )
342 {
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);
347 }
348 for( ; it != row1End; ++row0, ++it )
349 *it = (Field( k ) / Field( k+1 )) * (*row0);
350 row0 = row1;
351 }
352 }
353 };
354
355 template< class BaseTopology, class F >
356 class MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F >
357 {
358 typedef MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F > This;
359
360 public:
361 typedef Impl::Pyramid< BaseTopology > Topology;
362 typedef F Field;
363
364 static const unsigned int dimDomain = Topology::dimension;
365
366 typedef FieldVector< Field, dimDomain > DomainVector;
367
368 private:
369 friend class MonomialBasis< Topology, Field >;
370 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
371 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
372
373 typedef MonomialBasisSize< BaseTopology > BaseSize;
374 typedef MonomialBasisSize< Topology > Size;
375
376 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
377
378 MonomialBasisImpl ()
379 {}
380
381 template< int dimD >
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,
385 Field *const values,
386 const BaseSize &size ) const
387 {
388 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
389 }
390
391 template< int dimD >
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,
395 Field *const values,
396 const BaseSize &size ) const
397 {
398 Field omz = Unity< Field >() - x[ dimDomain-1 ];
399
400 if( Zero< Field >() < omz )
401 {
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;
406
407 // fill first column
408 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
409
410 Field omzk = omz;
411 for( unsigned int k = 1; k <= order; ++k )
412 {
413 Field *it = values + block*offsets[ k-1 ];
414 Field *const end = it + block*size.sizes_[ k ];
415 for( ; it != end; ++it )
416 *it *= omzk;
417 omzk *= omz;
418 }
419 }
420 else
421 {
422 assert( deriv==0 );
423 *values = Unity< Field >();
424 for( unsigned int k = 1; k <= order; ++k )
425 {
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 >();
430 }
431 }
432 }
433
434 template< int dimD >
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
439 {
440 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
441 const BaseSize &size = BaseSize::instance();
442 const_cast<BaseSize&>(size).computeSizes(order);
443
444 if( Impl::IsSimplex< Topology >::value )
445 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
446 else
447 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
448
449 Field *row0 = values;
450 for( unsigned int k = 1; k <= order; ++k )
451 {
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 ] );
455 row0 = row1;
456 }
457 }
458
459 void integrate ( const unsigned int order,
460 const unsigned int *const offsets,
461 Field *const values ) const
462 {
463 const BaseSize &size = BaseSize::instance();
464
465 // fill first column
466 baseBasis_.integrate( order, offsets, values );
467
468 const unsigned int *const baseSizes = size.sizes_;
469
470 {
471 Field *const col0End = values + baseSizes[ 0 ];
472 for( Field *it = values; it != col0End; ++it )
473 *it *= Field( 1 ) / Field( int(dimDomain) );
474 }
475
476 Field *row0 = values;
477 for( unsigned int k = 1; k <= order; ++k )
478 {
479 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
480
481 Field *const row1 = values+offsets[ k-1 ];
482 Field *const col0End = row1 + baseSizes[ k ];
483 Field *it = row1;
484 for( ; it != col0End; ++it )
485 *it *= factor;
486 for( unsigned int i = 1; i <= k; ++i )
487 {
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);
492 }
493 row0 = row1;
494 }
495 }
496 };
497
498
499
500 // MonomialBasis
501 // -------------
502
503 template< class Topology, class F >
504 class MonomialBasis
505 : public MonomialBasisImpl< Topology, F >
506 {
507 typedef MonomialBasis< Topology, F > This;
508 typedef MonomialBasisImpl< Topology, F > Base;
509
510 public:
511 static const unsigned int dimension = Base::dimDomain;
512 static const unsigned int dimRange = 1;
513
514 typedef typename Base::Field Field;
515
516 typedef typename Base::DomainVector DomainVector;
517
518 typedef Dune::FieldVector<Field,dimRange> RangeVector;
519
520 typedef MonomialBasisSize<Topology> Size;
521
522 MonomialBasis (unsigned int order)
523 : Base(),
524 order_(order),
525 size_(Size::instance())
526 {
527 assert(order<=1024); // avoid wrapping of unsigned int (0-1) order=1024 is quite hight...)
528 }
529
530 const unsigned int *sizes ( unsigned int order ) const
531 {
532 size_.computeSizes( order );
533 return size_.numBaseFunctions_;
534 }
535
536 const unsigned int *sizes () const
537 {
538 return sizes( order_ );
539 }
540
541 unsigned int size () const
542 {
543 size_.computeSizes( order_ );
544 return size_( order_ );
545 }
546
547 unsigned int derivSize ( const unsigned int deriv ) const
548 {
549 typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology;
550 MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
551 return MonomialBasisSize< SimplexTopology >::instance() ( deriv );
552 }
553
554 unsigned int order () const
555 {
556 return order_ ;
557 }
558
559 unsigned int topologyId ( ) const
560 {
561 return Topology::id;
562 }
563
564 void evaluate ( const unsigned int deriv, const DomainVector &x,
565 Field *const values ) const
566 {
567 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
568 }
569
570 template <unsigned int deriv>
571 void evaluate ( const DomainVector &x,
572 Field *const values ) const
573 {
574 evaluate( deriv, x, values );
575 }
576
577 template<unsigned int deriv, class Vector >
578 void evaluate ( const DomainVector &x,
579 Vector &values ) const
580 {
581 evaluate<deriv>(x,&(values[0]));
582 }
583 template<unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
584 void evaluate ( const DomainVector &x,
585 Derivatives<Field,dimension,1,deriv,layout> *values ) const
586 {
587 evaluate<deriv>(x,&(values->block()));
588 }
589 template< unsigned int deriv >
590 void evaluate ( const DomainVector &x,
591 FieldVector<Field,Derivatives<Field,dimension,1,deriv,DerivativeLayoutNS::value>::size> *values ) const
592 {
593 evaluate(0,x,&(values[0][0]));
594 }
595
596 template<class Vector >
597 void evaluate ( const DomainVector &x,
598 Vector &values ) const
599 {
600 evaluate<0>(x,&(values[0]));
601 }
602
603 template< class DVector, class RVector >
604 void evaluate ( const DVector &x, RVector &values ) const
605 {
606 assert( DVector::dimension == dimension);
607 DomainVector bx;
608 for( int d = 0; d < dimension; ++d )
609 field_cast( x[ d ], bx[ d ] );
610 evaluate<0>( bx, values );
611 }
612
613 void integrate ( Field *const values ) const
614 {
615 Base::integrate( order_, sizes( order_ ), values );
616 }
617 template <class Vector>
618 void integrate ( Vector &values ) const
619 {
620 integrate( &(values[ 0 ]) );
621 }
622 private:
623 MonomialBasis(const This&);
624 This& operator=(const This&);
625 unsigned int order_;
626 Size &size_;
627 };
628
629
630
631 // StdMonomialBasis
632 // ----------------
633
634 template< int dim,class F >
635 class StandardMonomialBasis
636 : public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F >
637 {
638 typedef StandardMonomialBasis< dim, F > This;
639 typedef MonomialBasis< typename Impl::SimplexTopology< dim >::type, F > Base;
640
641 public:
642 typedef typename Impl::SimplexTopology< dim >::type Topology;
643 static const int dimension = dim;
644
645 StandardMonomialBasis ( unsigned int order )
646 : Base( order )
647 {}
648 };
649
650
651
652 // StandardBiMonomialBasis
653 // -----------------------
654
655 template< int dim, class F >
656 class StandardBiMonomialBasis
657 : public MonomialBasis< typename Impl::CubeTopology< dim >::type, F >
658 {
659 typedef StandardBiMonomialBasis< dim, F > This;
660 typedef MonomialBasis< typename Impl::CubeTopology< dim >::type, F > Base;
661
662 public:
663 typedef typename Impl::CubeTopology< dim >::type Topology;
664 static const int dimension = dim;
665
666 StandardBiMonomialBasis ( unsigned int order )
667 : Base( order )
668 {}
669 };
670
671 // -----------------------------------------------------------
672 // -----------------------------------------------------------
673 // VirtualMonomialBasis
674 // -------------------
675
676 template< int dim, class F >
677 class VirtualMonomialBasis
678 {
679 typedef VirtualMonomialBasis< dim, F > This;
680
681 public:
682 typedef F Field;
683 typedef F StorageField;
684 static const int dimension = dim;
685 static const unsigned int dimRange = 1;
686
687 typedef FieldVector<Field,dimension> DomainVector;
688 typedef FieldVector<Field,dimRange> RangeVector;
689
690 explicit VirtualMonomialBasis(unsigned int topologyId,
691 unsigned int order)
692 : order_(order), topologyId_(topologyId) {}
693
694 virtual ~VirtualMonomialBasis() {}
695
696 virtual const unsigned int *sizes ( ) const = 0;
697
698 unsigned int size ( ) const
699 {
700 return sizes( )[ order_ ];
701 }
702
703 unsigned int order () const
704 {
705 return order_;
706 }
707
708 unsigned int topologyId ( ) const
709 {
710 return topologyId_;
711 }
712
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
718 {
719 evaluate( deriv, x, values );
720 }
721 template < unsigned int deriv, int size >
722 void evaluate ( const DomainVector &x,
723 Dune::FieldVector<Field,size> *const values ) const
724 {
725 evaluate( deriv, x, &(values[0][0]) );
726 }
727 template<unsigned int deriv, DerivativeLayoutNS::DerivativeLayout layout >
728 void evaluate ( const DomainVector &x,
729 Derivatives<Field,dimension,1,deriv,layout> *values ) const
730 {
731 evaluate<deriv>(x,&(values->block()));
732 }
733 template <unsigned int deriv, class Vector>
734 void evaluate ( const DomainVector &x,
735 Vector &values ) const
736 {
737 evaluate<deriv>( x, &(values[ 0 ]) );
738 }
739 template< class Vector >
740 void evaluate ( const DomainVector &x,
741 Vector &values ) const
742 {
743 evaluate<0>(x,values);
744 }
745 template< class DVector, class RVector >
746 void evaluate ( const DVector &x, RVector &values ) const
747 {
748 assert( DVector::dimension == dimension);
749 DomainVector bx;
750 for( int d = 0; d < dimension; ++d )
751 field_cast( x[ d ], bx[ d ] );
752 evaluate<0>( bx, values );
753 }
754 template< unsigned int deriv, class DVector, class RVector >
755 void evaluate ( const DVector &x, RVector &values ) const
756 {
757 assert( DVector::dimension == dimension);
758 DomainVector bx;
759 for( int d = 0; d < dimension; ++d )
760 field_cast( x[ d ], bx[ d ] );
761 evaluate<deriv>( bx, values );
762 }
763
764 virtual void integrate ( Field *const values ) const = 0;
765 template <class Vector>
766 void integrate ( Vector &values ) const
767 {
768 integrate( &(values[ 0 ]) );
769 }
770 protected:
771 unsigned int order_;
772 unsigned int topologyId_;
773 };
774
775 template< class Topology, class F >
776 class VirtualMonomialBasisImpl
777 : public VirtualMonomialBasis< Topology::dimension, F >
778 {
779 typedef VirtualMonomialBasis< Topology::dimension, F > Base;
780 typedef VirtualMonomialBasisImpl< Topology, F > This;
781
782 public:
783 typedef typename Base::Field Field;
784 typedef typename Base::DomainVector DomainVector;
785
786 VirtualMonomialBasisImpl(unsigned int order)
787 : Base(Topology::id,order), basis_(order)
788 {}
789
790 const unsigned int *sizes ( ) const
791 {
792 return basis_.sizes(order_);
793 }
794
795 void evaluate ( const unsigned int deriv, const DomainVector &x,
796 Field *const values ) const
797 {
798 basis_.evaluate(deriv,x,values);
799 }
800
801 void integrate ( Field *const values ) const
802 {
803 basis_.integrate(values);
804 }
805
806 private:
807 MonomialBasis<Topology,Field> basis_;
808 using Base::order_;
809 };
810
811 // MonomialBasisFactory
812 // --------------------
813
814 template< int dim, class F >
815 struct MonomialBasisFactory
816 {
817 static const unsigned int dimension = dim;
818 typedef F StorageField;
819
820 typedef unsigned int Key;
821 typedef const VirtualMonomialBasis< dimension, F > Object;
822
823 template < int dd, class FF >
824 struct EvaluationBasisFactory
825 {
826 typedef MonomialBasisFactory<dd,FF> Type;
827 };
828
829 template< class Topology >
830 static Object* create ( const Key &order )
831 {
832 return new VirtualMonomialBasisImpl< Topology, StorageField >( order );
833 }
834 static void release( Object *object ) { delete object; }
835 };
836
837
838
839 // MonomialBasisProvider
840 // ---------------------
841
842 template< int dim, class SF >
843 struct MonomialBasisProvider
844 : public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
845 {
846 static const unsigned int dimension = dim;
847 typedef SF StorageField;
848 template < int dd, class FF >
849 struct EvaluationBasisFactory
850 {
851 typedef MonomialBasisProvider<dd,FF> Type;
852 };
853 };
854
855}
856
857#endif
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:96
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Dune::field_cast
void field_cast(const F1 &f1, F2 &f2)
a helper class to cast from one field to another
Definition: field.hh:157
type.hh
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)