Dune Core Modules (2.6.0)

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<>
81 class MonomialBasisSize< Impl::Point >
82 {
83 typedef MonomialBasisSize< Impl::Point > This;
84
85 public:
86 static This &instance ()
87 {
88 static This _instance;
89 return _instance;
90 }
91
92 typedef Impl::Point Topology;
93
94 friend class MonomialBasisSize< Impl::Prism< Topology > >;
95 friend class MonomialBasisSize< Impl::Pyramid< Topology > >;
96
97 mutable unsigned int maxOrder_;
98 // sizes_[ k ]: number of basis functions of exactly order k
99 mutable unsigned int *sizes_;
100 // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
101 mutable unsigned int *numBaseFunctions_;
102
103 MonomialBasisSize ()
104 : maxOrder_( 0 ),
105 sizes_( 0 ),
106 numBaseFunctions_( 0 )
107 {
108 computeSizes( 2 );
109 }
110
111 ~MonomialBasisSize ()
112 {
113 delete[] sizes_;
114 delete[] numBaseFunctions_;
115 }
116
117 unsigned int operator() ( const unsigned int order ) const
118 {
119 return numBaseFunctions_[ order ];
120 }
121
122 unsigned int maxOrder () const
123 {
124 return maxOrder_;
125 }
126
127 void computeSizes ( unsigned int order ) const
128 {
129 if (order <= maxOrder_)
130 return;
131
132 maxOrder_ = order;
133
134 delete [] sizes_;
135 delete [] numBaseFunctions_;
136 sizes_ = new unsigned int [ order+1 ];
137 numBaseFunctions_ = new unsigned int [ order+1 ];
138
139 sizes_[ 0 ] = 1;
140 numBaseFunctions_[ 0 ] = 1;
141 for( unsigned int k = 1; k <= order; ++k )
142 {
143 sizes_[ k ] = 0;
144 numBaseFunctions_[ k ] = 1;
145 }
146 }
147 };
148
149 template< class BaseTopology >
150 class MonomialBasisSize< Impl::Prism< BaseTopology > >
151 {
152 typedef MonomialBasisSize< Impl::Prism< BaseTopology > > This;
153
154 public:
155 static This &instance ()
156 {
157 static This _instance;
158 return _instance;
159 }
160
161 typedef Impl::Prism< BaseTopology > Topology;
162
163 friend class MonomialBasisSize< Impl::Prism< Topology > >;
164 friend class MonomialBasisSize< Impl::Pyramid< Topology > >;
165
166 mutable unsigned int maxOrder_;
167 // sizes_[ k ]: number of basis functions of exactly order k
168 mutable unsigned int *sizes_;
169 // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
170 mutable unsigned int *numBaseFunctions_;
171
172 MonomialBasisSize ()
173 : maxOrder_( 0 ),
174 sizes_( 0 ),
175 numBaseFunctions_( 0 )
176 {
177 computeSizes( 2 );
178 }
179
180 ~MonomialBasisSize ()
181 {
182 delete[] sizes_;
183 delete[] numBaseFunctions_;
184 }
185
186 unsigned int operator() ( const unsigned int order ) const
187 {
188 return numBaseFunctions_[ order ];
189 }
190
191 unsigned int maxOrder() const
192 {
193 return maxOrder_;
194 }
195
196 void computeSizes ( unsigned int order ) const
197 {
198 if (order <= maxOrder_)
199 return;
200
201 maxOrder_ = order;
202
203 delete[] sizes_;
204 delete[] numBaseFunctions_;
205 sizes_ = new unsigned int[ order+1 ];
206 numBaseFunctions_ = new unsigned int[ order+1 ];
207
208 MonomialBasisSize<BaseTopology> &baseBasis =
209 MonomialBasisSize<BaseTopology>::instance();
210 baseBasis.computeSizes( order );
211 const unsigned int *const baseSizes = baseBasis.sizes_;
212 const unsigned int *const baseNBF = baseBasis.numBaseFunctions_;
213
214 sizes_[ 0 ] = 1;
215 numBaseFunctions_[ 0 ] = 1;
216 for( unsigned int k = 1; k <= order; ++k )
217 {
218 sizes_[ k ] = baseNBF[ k ] + k*baseSizes[ k ];
219 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
220 }
221 }
222 };
223
224 template< class BaseTopology >
225 class MonomialBasisSize< Impl::Pyramid< BaseTopology > >
226 {
227 typedef MonomialBasisSize< Impl::Pyramid< BaseTopology > > This;
228
229 public:
230 static This &instance ()
231 {
232 static This _instance;
233 return _instance;
234 }
235
236 typedef Impl::Pyramid< BaseTopology > Topology;
237
238 friend class MonomialBasisSize< Impl::Prism< Topology > >;
239 friend class MonomialBasisSize< Impl::Pyramid< Topology > >;
240
241 mutable unsigned int maxOrder_;
242 // sizes_[ k ]: number of basis functions of exactly order k
243 mutable unsigned int *sizes_;
244 // numBaseFunctions_[ k ] = sizes_[ 0 ] + ... + sizes_[ k ]
245 mutable unsigned int *numBaseFunctions_;
246
247 MonomialBasisSize ()
248 : maxOrder_( 0 ),
249 sizes_( 0 ),
250 numBaseFunctions_( 0 )
251 {
252 computeSizes( 2 );
253 }
254
255 ~MonomialBasisSize ()
256 {
257 delete[] sizes_;
258 delete[] numBaseFunctions_;
259 }
260
261 unsigned int operator() ( const unsigned int order ) const
262 {
263 return numBaseFunctions_[ order ];
264 }
265
266 unsigned int maxOrder() const
267 {
268 return maxOrder_;
269 }
270
271 void computeSizes ( unsigned int order ) const
272 {
273 if (order <= maxOrder_)
274 return;
275
276 maxOrder_ = order;
277
278 delete[] sizes_;
279 delete[] numBaseFunctions_;
280 sizes_ = new unsigned int[ order+1 ];
281 numBaseFunctions_ = new unsigned int[ order+1 ];
282
283 MonomialBasisSize<BaseTopology> &baseBasis =
284 MonomialBasisSize<BaseTopology>::instance();
285
286 baseBasis.computeSizes( order );
287
288 const unsigned int *const baseNBF = baseBasis.numBaseFunctions_;
289 sizes_[ 0 ] = 1;
290 numBaseFunctions_[ 0 ] = 1;
291 for( unsigned int k = 1; k <= order; ++k )
292 {
293 sizes_[ k ] = baseNBF[ k ];
294 numBaseFunctions_[ k ] = numBaseFunctions_[ k-1 ] + sizes_[ k ];
295 }
296 }
297 };
298
299
300
301 // MonomialBasisHelper
302 // -------------------
303
304
305 template< int mydim, int dim, class F >
306 struct MonomialBasisHelper
307 {
308 typedef MonomialBasisSize< typename Impl::SimplexTopology< mydim >::type > MySize;
309 typedef MonomialBasisSize< typename Impl::SimplexTopology< dim >::type > Size;
310
311 static void copy ( const unsigned int deriv, F *&wit, F *&rit,
312 const unsigned int numBaseFunctions, const F &z )
313 {
314 // n(d,k) = size<k>[d];
315 MySize &mySize = MySize::instance();
316 Size &size = Size::instance();
317
318 const F *const rend = rit + size( deriv )*numBaseFunctions;
319 for( ; rit != rend; )
320 {
321 F *prit = rit;
322
323 *wit = z * *rit;
324 ++rit, ++wit;
325
326 for( unsigned d = 1; d <= deriv; ++d )
327 {
328 #ifndef NDEBUG
329 const F *const derivEnd = rit + mySize.sizes_[ d ];
330 #endif
331
332 {
333 const F *const drend = rit + mySize.sizes_[ d ] - mySize.sizes_[ d-1 ];
334 for( ; rit != drend ; ++rit, ++wit )
335 *wit = z * *rit;
336 }
337
338 for (unsigned int j=1; j<d; ++j)
339 {
340 const F *const drend = rit + mySize.sizes_[ d-j ] - mySize.sizes_[ d-j-1 ];
341 for( ; rit != drend ; ++prit, ++rit, ++wit )
342 *wit = F(j) * *prit + z * *rit;
343 }
344 *wit = F(d) * *prit + z * *rit;
345 ++prit, ++rit, ++wit;
346 assert(derivEnd == rit);
347 rit += size.sizes_[d] - mySize.sizes_[d];
348 prit += size.sizes_[d-1] - mySize.sizes_[d-1];
349 const F *const emptyWitEnd = wit + size.sizes_[d] - mySize.sizes_[d];
350 for ( ; wit != emptyWitEnd; ++wit )
351 *wit = Zero<F>();
352 }
353 }
354 }
355 };
356
357
358
359 // MonomialBasisImpl
360 // -----------------
361
362 template< class Topology, class F >
363 class MonomialBasisImpl;
364
365 template< class F >
366 class MonomialBasisImpl< Impl::Point, F >
367 {
368 typedef MonomialBasisImpl< Impl::Point, F > This;
369
370 public:
371 typedef Impl::Point Topology;
372 typedef F Field;
373
374 static const unsigned int dimDomain = Topology::dimension;
375
376 typedef FieldVector< Field, dimDomain > DomainVector;
377
378 private:
379 friend class MonomialBasis< Topology, Field >;
380 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
381 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
382
383 template< int dimD >
384 void evaluate ( const unsigned int deriv, const unsigned int order,
385 const FieldVector< Field, dimD > &x,
386 const unsigned int block, const unsigned int *const offsets,
387 Field *const values ) const
388 {
389 *values = Unity< F >();
390 F *const end = values + block;
391 for( Field *it = values+1 ; it != end; ++it )
392 *it = Zero< F >();
393 }
394
395 void integrate ( const unsigned int order,
396 const unsigned int *const offsets,
397 Field *const values ) const
398 {
399 values[ 0 ] = Unity< Field >();
400 }
401 };
402
403 template< class BaseTopology, class F >
404 class MonomialBasisImpl< Impl::Prism< BaseTopology >, F >
405 {
406 typedef MonomialBasisImpl< Impl::Prism< BaseTopology >, F > This;
407
408 public:
409 typedef Impl::Prism< BaseTopology > Topology;
410 typedef F Field;
411
412 static const unsigned int dimDomain = Topology::dimension;
413
414 typedef FieldVector< Field, dimDomain > DomainVector;
415
416 private:
417 friend class MonomialBasis< Topology, Field >;
418 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
419 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
420
421 typedef MonomialBasisSize< BaseTopology > BaseSize;
422 typedef MonomialBasisSize< Topology > Size;
423
424 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
425
426 MonomialBasisImpl ()
427 {}
428
429 template< int dimD >
430 void evaluate ( const unsigned int deriv, const unsigned int order,
431 const FieldVector< Field, dimD > &x,
432 const unsigned int block, const unsigned int *const offsets,
433 Field *const values ) const
434 {
435 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
436 const BaseSize &size = BaseSize::instance();
437
438 const Field &z = x[ dimDomain-1 ];
439
440 // fill first column
441 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
442
443 Field *row0 = values;
444 for( unsigned int k = 1; k <= order; ++k )
445 {
446 Field *row1 = values + block*offsets[ k-1 ];
447 Field *wit = row1 + block*size.sizes_[ k ];
448 Helper::copy( deriv, wit, row1, k*size.sizes_[ k ], z );
449 Helper::copy( deriv, wit, row0, size( k-1 ), z );
450 row0 = row1;
451 }
452 }
453
454 void integrate ( const unsigned int order,
455 const unsigned int *const offsets,
456 Field *const values ) const
457 {
458 const BaseSize &size = BaseSize::instance();
459 const Size &mySize = Size::instance();
460 // fill first column
461 baseBasis_.integrate( order, offsets, values );
462 const unsigned int *const baseSizes = size.sizes_;
463
464 Field *row0 = values;
465 for( unsigned int k = 1; k <= order; ++k )
466 {
467 Field *const row1begin = values + offsets[ k-1 ];
468 Field *const row1End = row1begin + mySize.sizes_[ k ];
469 assert( (unsigned int)(row1End - values) <= offsets[ k ] );
470
471 Field *row1 = row1begin;
472 Field *it = row1begin + baseSizes[ k ];
473 for( unsigned int j = 1; j <= k; ++j )
474 {
475 Field *const end = it + baseSizes[ k ];
476 assert( (unsigned int)(end - values) <= offsets[ k ] );
477 for( ; it != end; ++row1, ++it )
478 *it = (Field( j ) / Field( j+1 )) * (*row1);
479 }
480 for( ; it != row1End; ++row0, ++it )
481 *it = (Field( k ) / Field( k+1 )) * (*row0);
482 row0 = row1;
483 }
484 }
485 };
486
487 template< class BaseTopology, class F >
488 class MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F >
489 {
490 typedef MonomialBasisImpl< Impl::Pyramid< BaseTopology >, F > This;
491
492 public:
493 typedef Impl::Pyramid< BaseTopology > Topology;
494 typedef F Field;
495
496 static const unsigned int dimDomain = Topology::dimension;
497
498 typedef FieldVector< Field, dimDomain > DomainVector;
499
500 private:
501 friend class MonomialBasis< Topology, Field >;
502 friend class MonomialBasisImpl< Impl::Prism< Topology >, Field >;
503 friend class MonomialBasisImpl< Impl::Pyramid< Topology >, Field >;
504
505 typedef MonomialBasisSize< BaseTopology > BaseSize;
506 typedef MonomialBasisSize< Topology > Size;
507
508 MonomialBasisImpl< BaseTopology, Field > baseBasis_;
509
510 MonomialBasisImpl ()
511 {}
512
513 template< int dimD >
514 void evaluateSimplexBase ( const unsigned int deriv, const unsigned int order,
515 const FieldVector< Field, dimD > &x,
516 const unsigned int block, const unsigned int *const offsets,
517 Field *const values,
518 const BaseSize &size ) const
519 {
520 baseBasis_.evaluate( deriv, order, x, block, offsets, values );
521 }
522
523 template< int dimD >
524 void evaluatePyramidBase ( const unsigned int deriv, const unsigned int order,
525 const FieldVector< Field, dimD > &x,
526 const unsigned int block, const unsigned int *const offsets,
527 Field *const values,
528 const BaseSize &size ) const
529 {
530 Field omz = Unity< Field >() - x[ dimDomain-1 ];
531
532 if( Zero< Field >() < omz )
533 {
534 const Field invomz = Unity< Field >() / omz;
535 FieldVector< Field, dimD > y;
536 for( unsigned int i = 0; i < dimDomain-1; ++i )
537 y[ i ] = x[ i ] * invomz;
538
539 // fill first column
540 baseBasis_.evaluate( deriv, order, y, block, offsets, values );
541
542 Field omzk = omz;
543 for( unsigned int k = 1; k <= order; ++k )
544 {
545 Field *it = values + block*offsets[ k-1 ];
546 Field *const end = it + block*size.sizes_[ k ];
547 for( ; it != end; ++it )
548 *it *= omzk;
549 omzk *= omz;
550 }
551 }
552 else
553 {
554 assert( deriv==0 );
555 *values = Unity< Field >();
556 for( unsigned int k = 1; k <= order; ++k )
557 {
558 Field *it = values + block*offsets[ k-1 ];
559 Field *const end = it + block*size.sizes_[ k ];
560 for( ; it != end; ++it )
561 *it = Zero< Field >();
562 }
563 }
564 }
565
566 template< int dimD >
567 void evaluate ( const unsigned int deriv, const unsigned int order,
568 const FieldVector< Field, dimD > &x,
569 const unsigned int block, const unsigned int *const offsets,
570 Field *const values ) const
571 {
572 typedef MonomialBasisHelper< dimDomain, dimD, Field > Helper;
573 const BaseSize &size = BaseSize::instance();
574
575 if( Impl::IsSimplex< Topology >::value )
576 evaluateSimplexBase( deriv, order, x, block, offsets, values, size );
577 else
578 evaluatePyramidBase( deriv, order, x, block, offsets, values, size );
579
580 Field *row0 = values;
581 for( unsigned int k = 1; k <= order; ++k )
582 {
583 Field *row1 = values + block*offsets[ k-1 ];
584 Field *wit = row1 + block*size.sizes_[ k ];
585 Helper::copy( deriv, wit, row0, size( k-1 ), x[ dimDomain-1 ] );
586 row0 = row1;
587 }
588 }
589
590 void integrate ( const unsigned int order,
591 const unsigned int *const offsets,
592 Field *const values ) const
593 {
594 const BaseSize &size = BaseSize::instance();
595
596 // fill first column
597 baseBasis_.integrate( order, offsets, values );
598
599 const unsigned int *const baseSizes = size.sizes_;
600
601 {
602 Field *const col0End = values + baseSizes[ 0 ];
603 for( Field *it = values; it != col0End; ++it )
604 *it *= Field( 1 ) / Field( int(dimDomain) );
605 }
606
607 Field *row0 = values;
608 for( unsigned int k = 1; k <= order; ++k )
609 {
610 const Field factor = (Field( 1 ) / Field( k + dimDomain ));
611
612 Field *const row1 = values+offsets[ k-1 ];
613 Field *const col0End = row1 + baseSizes[ k ];
614 Field *it = row1;
615 for( ; it != col0End; ++it )
616 *it *= factor;
617 for( unsigned int i = 1; i <= k; ++i )
618 {
619 Field *const end = it + baseSizes[ k-i ];
620 assert( (unsigned int)(end - values) <= offsets[ k ] );
621 for( ; it != end; ++row0, ++it )
622 *it = (*row0) * (Field( i ) * factor);
623 }
624 row0 = row1;
625 }
626 }
627 };
628
629
630
631 // MonomialBasis
632 // -------------
633
634 template< class Topology, class F >
635 class MonomialBasis
636 : public MonomialBasisImpl< Topology, F >
637 {
638 typedef MonomialBasis< Topology, F > This;
639 typedef MonomialBasisImpl< Topology, F > Base;
640
641 public:
642 static const unsigned int dimension = Base::dimDomain;
643 static const unsigned int dimRange = 1;
644
645 typedef typename Base::Field Field;
646
647 typedef typename Base::DomainVector DomainVector;
648
649 typedef Dune::FieldVector<Field,dimRange> RangeVector;
650
651 typedef MonomialBasisSize<Topology> Size;
652
653 MonomialBasis (unsigned int order)
654 : Base(),
655 order_(order),
656 size_(Size::instance())
657 {
658 assert(order<=1024); // avoid wrapping of unsigned int (0-1) order=1024 is quite hight...)
659 }
660
661 const unsigned int *sizes ( unsigned int order ) const
662 {
663 size_.computeSizes( order );
664 return size_.numBaseFunctions_;
665 }
666
667 const unsigned int *sizes () const
668 {
669 return sizes( order_ );
670 }
671
672 unsigned int size () const
673 {
674 size_.computeSizes( order_ );
675 return size_( order_ );
676 }
677
678 unsigned int derivSize ( const unsigned int deriv ) const
679 {
680 typedef typename Impl::SimplexTopology< dimension >::type SimplexTopology;
681 MonomialBasisSize< SimplexTopology >::instance().computeSizes( deriv );
682 return MonomialBasisSize< SimplexTopology >::instance() ( deriv );
683 }
684
685 unsigned int order () const
686 {
687 return order_ ;
688 }
689
690 unsigned int topologyId ( ) const
691 {
692 return Topology::id;
693 }
694
695 void evaluate ( const unsigned int deriv, const DomainVector &x,
696 Field *const values ) const
697 {
698 Base::evaluate( deriv, order_, x, derivSize( deriv ), sizes( order_ ), values );
699 }
700
701 template <unsigned int deriv>
702 void evaluate ( const DomainVector &x,
703 Field *const values ) const
704 {
705 evaluate( deriv, x, values );
706 }
707
708 template<unsigned int deriv, class Vector >
709 void evaluate ( const DomainVector &x,
710 Vector &values ) const
711 {
712 evaluate<deriv>(x,&(values[0]));
713 }
714 template<unsigned int deriv, DerivativeLayout layout >
715 void evaluate ( const DomainVector &x,
716 Derivatives<Field,dimension,1,deriv,layout> *values ) const
717 {
718 evaluate<deriv>(x,&(values->block()));
719 }
720 template< unsigned int deriv >
721 void evaluate ( const DomainVector &x,
722 FieldVector<Field,Derivatives<Field,dimension,1,deriv,value>::size> *values ) const
723 {
724 evaluate(0,x,&(values[0][0]));
725 }
726
727 template<class Vector >
728 void evaluate ( const DomainVector &x,
729 Vector &values ) const
730 {
731 evaluate<0>(x,&(values[0]));
732 }
733
734 template< class DVector, class RVector >
735 void evaluate ( const DVector &x, RVector &values ) const
736 {
737 assert( DVector::dimension == dimension);
738 DomainVector bx;
739 for( int d = 0; d < dimension; ++d )
740 field_cast( x[ d ], bx[ d ] );
741 evaluate<0>( bx, values );
742 }
743
744 void integrate ( Field *const values ) const
745 {
746 Base::integrate( order_, sizes( order_ ), values );
747 }
748 template <class Vector>
749 void integrate ( Vector &values ) const
750 {
751 integrate( &(values[ 0 ]) );
752 }
753 private:
754 MonomialBasis(const This&);
755 This& operator=(const This&);
756 unsigned int order_;
757 Size &size_;
758 };
759
760
761
762 // StdMonomialBasis
763 // ----------------
764
765 template< int dim,class F >
766 class StandardMonomialBasis
767 : public MonomialBasis< typename Impl::SimplexTopology< dim >::type, F >
768 {
769 typedef StandardMonomialBasis< dim, F > This;
770 typedef MonomialBasis< typename Impl::SimplexTopology< dim >::type, F > Base;
771
772 public:
773 typedef typename Impl::SimplexTopology< dim >::type Topology;
774 static const int dimension = dim;
775
776 StandardMonomialBasis ( unsigned int order )
777 : Base( order )
778 {}
779 };
780
781
782
783 // StandardBiMonomialBasis
784 // -----------------------
785
786 template< int dim, class F >
787 class StandardBiMonomialBasis
788 : public MonomialBasis< typename Impl::CubeTopology< dim >::type, F >
789 {
790 typedef StandardBiMonomialBasis< dim, F > This;
791 typedef MonomialBasis< typename Impl::CubeTopology< dim >::type, F > Base;
792
793 public:
794 typedef typename Impl::CubeTopology< dim >::type Topology;
795 static const int dimension = dim;
796
797 StandardBiMonomialBasis ( unsigned int order )
798 : Base( order )
799 {}
800 };
801
802 // -----------------------------------------------------------
803 // -----------------------------------------------------------
804 // VirtualMonomialBasis
805 // -------------------
806
807 template< int dim, class F >
808 class VirtualMonomialBasis
809 {
810 typedef VirtualMonomialBasis< dim, F > This;
811
812 public:
813 typedef F Field;
814 typedef F StorageField;
815 static const int dimension = dim;
816 static const unsigned int dimRange = 1;
817
818 typedef FieldVector<Field,dimension> DomainVector;
819 typedef FieldVector<Field,dimRange> RangeVector;
820
821 explicit VirtualMonomialBasis(unsigned int topologyId,
822 unsigned int order)
823 : order_(order), topologyId_(topologyId) {}
824
825 virtual ~VirtualMonomialBasis() {}
826
827 virtual const unsigned int *sizes ( ) const = 0;
828
829 unsigned int size ( ) const
830 {
831 return sizes( )[ order_ ];
832 }
833
834 unsigned int order () const
835 {
836 return order_;
837 }
838
839 unsigned int topologyId ( ) const
840 {
841 return topologyId_;
842 }
843
844 virtual void evaluate ( const unsigned int deriv, const DomainVector &x,
845 Field *const values ) const = 0;
846 template < unsigned int deriv >
847 void evaluate ( const DomainVector &x,
848 Field *const values ) const
849 {
850 evaluate( deriv, x, values );
851 }
852 template < unsigned int deriv, int size >
853 void evaluate ( const DomainVector &x,
854 Dune::FieldVector<Field,size> *const values ) const
855 {
856 evaluate( deriv, x, &(values[0][0]) );
857 }
858 template<unsigned int deriv, DerivativeLayout layout >
859 void evaluate ( const DomainVector &x,
860 Derivatives<Field,dimension,1,deriv,layout> *values ) const
861 {
862 evaluate<deriv>(x,&(values->block()));
863 }
864 template <unsigned int deriv, class Vector>
865 void evaluate ( const DomainVector &x,
866 Vector &values ) const
867 {
868 evaluate<deriv>( x, &(values[ 0 ]) );
869 }
870 template< class Vector >
871 void evaluate ( const DomainVector &x,
872 Vector &values ) const
873 {
874 evaluate<0>(x,values);
875 }
876 template< class DVector, class RVector >
877 void evaluate ( const DVector &x, RVector &values ) const
878 {
879 assert( DVector::dimension == dimension);
880 DomainVector bx;
881 for( int d = 0; d < dimension; ++d )
882 field_cast( x[ d ], bx[ d ] );
883 evaluate<0>( bx, values );
884 }
885 template< unsigned int deriv, class DVector, class RVector >
886 void evaluate ( const DVector &x, RVector &values ) const
887 {
888 assert( DVector::dimension == dimension);
889 DomainVector bx;
890 for( int d = 0; d < dimension; ++d )
891 field_cast( x[ d ], bx[ d ] );
892 evaluate<deriv>( bx, values );
893 }
894
895 virtual void integrate ( Field *const values ) const = 0;
896 template <class Vector>
897 void integrate ( Vector &values ) const
898 {
899 integrate( &(values[ 0 ]) );
900 }
901 protected:
902 unsigned int order_;
903 unsigned int topologyId_;
904 };
905
906 template< class Topology, class F >
907 class VirtualMonomialBasisImpl
908 : public VirtualMonomialBasis< Topology::dimension, F >
909 {
910 typedef VirtualMonomialBasis< Topology::dimension, F > Base;
911 typedef VirtualMonomialBasisImpl< Topology, F > This;
912
913 public:
914 typedef typename Base::Field Field;
915 typedef typename Base::DomainVector DomainVector;
916
917 VirtualMonomialBasisImpl(unsigned int order)
918 : Base(Topology::id,order), basis_(order)
919 {}
920
921 const unsigned int *sizes ( ) const
922 {
923 return basis_.sizes(order_);
924 }
925
926 void evaluate ( const unsigned int deriv, const DomainVector &x,
927 Field *const values ) const
928 {
929 basis_.evaluate(deriv,x,values);
930 }
931
932 void integrate ( Field *const values ) const
933 {
934 basis_.integrate(values);
935 }
936
937 private:
938 MonomialBasis<Topology,Field> basis_;
939 using Base::order_;
940 };
941
942 // MonomialBasisFactory
943 // --------------------
944
945 template< int dim, class F >
946 struct MonomialBasisFactory;
947
948 template< int dim, class F >
949 struct MonomialBasisFactoryTraits
950 {
951 static const unsigned int dimension = dim;
952 typedef unsigned int Key;
953 typedef const VirtualMonomialBasis< dimension, F > Object;
954 typedef MonomialBasisFactory<dim,F> Factory;
955 };
956
957 template< int dim, class F >
958 struct MonomialBasisFactory :
959 public TopologyFactory< MonomialBasisFactoryTraits<dim,F> >
960 {
961 static const unsigned int dimension = dim;
962 typedef F StorageField;
963 typedef MonomialBasisFactoryTraits<dim,F> Traits;
964
965 typedef typename Traits::Key Key;
966 typedef typename Traits::Object Object;
967
968 template < int dd, class FF >
969 struct EvaluationBasisFactory
970 {
971 typedef MonomialBasisFactory<dd,FF> Type;
972 };
973
974 template< class Topology >
975 static Object* createObject ( const Key &order )
976 {
977 return new VirtualMonomialBasisImpl< Topology, StorageField >( order );
978 }
979 };
980
981
982
983 // MonomialBasisProvider
984 // ---------------------
985
986 template< int dim, class SF >
987 struct MonomialBasisProvider
988 : public TopologySingletonFactory< MonomialBasisFactory< dim, SF > >
989 {
990 static const unsigned int dimension = dim;
991 typedef SF StorageField;
992 template < int dd, class FF >
993 struct EvaluationBasisFactory
994 {
995 typedef MonomialBasisProvider<dd,FF> Type;
996 };
997 };
998
999}
1000
1001#endif
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
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:10
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)