Dune Core Modules (unstable)

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