DUNE PDELab (git)

l2orthonormal.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3
4#ifndef DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
5#define DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
6
7#include<array>
8#include<iostream>
9#include<algorithm>
10#include<memory>
11#include<numeric>
12
18
19#include<dune/geometry/referenceelements.hh>
22
23#include<dune/localfunctions/common/localbasis.hh>
24#include<dune/localfunctions/common/localkey.hh>
25#include<dune/localfunctions/common/localfiniteelementtraits.hh>
26#include<dune/localfunctions/common/localinterpolation.hh>
27
38namespace Dune {
39 namespace PB {
40
41#if HAVE_GMP
42 typedef GMPField<512> DefaultComputationalFieldType;
43#else
44 typedef double DefaultComputationalFieldType;
45#endif
46
47 //=====================================================
48 // Type to represent a multiindex in d dimensions
49 //=====================================================
50
51 template<int d>
52 class MultiIndex : public std::array<int,d>
53 {
54 public:
55
56 MultiIndex () : std::array<int,d>()
57 {
58 }
59
60 };
61
62 // the number of polynomials of at most degree k in space dimension d
63 constexpr int pk_size (int k, int d)
64 {
65 if (k==0) return 1;
66 if (d==1) return k+1;
67 int n=0;
68 for (int j=0; j<=k; j++)
69 n += pk_size(k-j,d-1);
70 return n;
71 }
72
73 // the number of polynomials of exactly degree k in space dimension d (as run-time function)
74 inline int pk_size_exact (int k, int d)
75 {
76 if (k==0)
77 return 1;
78 else
79 return pk_size(k,d)-pk_size(k-1,d);
80 }
81
82 // k is the polynomial degree and d is the space dimension
83 // then PkSize<k,d> is the number of polynomials of at most total degree k.
84 template<int k, int d>
85 struct PkSize
86 {
87 enum{
88 value=pk_size(k,d)
89 };
90 };
91
92 // enumerate all multiindices of degree k and find the i'th
93 template<int d>
94 void pk_enumerate_multiindex (MultiIndex<d>& alpha, int& count, int norm, int dim, int k, int i)
95 {
96 // check if dim is valid
97 if (dim>=d) return;
98
99 // put all k to current dimension and check
100 alpha[dim]=k; count++; // alpha has correct norm
101 // std::cout << "dadi alpha=" << alpha << " count=" << count << " norm=" << norm+k << " dim=" << dim << " k=" << k << " i=" << i << std::endl;
102 if (count==i) return; // found the index
103
104 // search recursively
105 for (int m=k-1; m>=0; m--)
106 {
107 alpha[dim]=m;
108 //std::cout << "dada alpha=" << alpha << " count=" << count << " norm=" << norm+m << " dim=" << dim << " k=" << k << " i=" << i << std::endl;
109 pk_enumerate_multiindex(alpha,count,norm+m,dim+1,k-m,i);
110 if (count==i) return;
111 }
112
113 // reset if index is not yet found
114 alpha[dim]=0;
115 }
116
117 // map integer 0<=i<pk_size(k,d) to multiindex
118 template<int d>
119 void pk_multiindex (int i, MultiIndex<d>& alpha)
120 {
121 for (int j=0; j<d; j++) alpha[j] = 0; // set alpha to 0
122 if (i==0) return; // handle case i==0 and k==0
123 for (int k=1; k<25; k++)
124 if (i>=pk_size(k-1,d) && i<pk_size(k,d))
125 {
126 int count = -1;
127 pk_enumerate_multiindex(alpha,count,0,0,k,i-pk_size(k-1,d));
128 return;
129 }
130 DUNE_THROW(Dune::NotImplemented,"Polynomial degree greater than 24 in pk_multiindex");
131 }
132
133 // the number of polynomials of at most degree k in space dimension d (as run-time function)
134 constexpr int qk_size (int k, int d)
135 {
136 int n = 1;
137 for (int i=0; i<d; ++i)
138 n *= (k+1);
139 return n;
140 }
141
142 // map integer 0<=i<qk_size(k,d) to multiindex
143 template<int d>
144 void qk_multiindex (int i, int k, MultiIndex<d>& alpha)
145 {
146 for (int j = 0; j<d; ++j) {
147 alpha[j] = i % (k+1);
148 i /= (k+1);
149 }
150 }
151
152 //=====================================================
153 // Traits classes to group Pk and Qk specifics
154 //=====================================================
155 enum BasisType {
156 Pk, Qk
157 };
158
159 template <BasisType basisType>
160 struct BasisTraits;
161
162 template <>
163 struct BasisTraits<BasisType::Pk>
164 {
165 template <int k, int d>
166 struct Size
167 {
168 enum{
169 value = pk_size(k,d)
170 };
171 };
172
173 template <int k, int d>
174 struct Order
175 {
176 enum{
177 value = k
178 };
179 };
180
181 static int size(int k, int d)
182 {
183 return pk_size(k,d);
184 }
185
186 template <int d>
187 static void multiindex(int i, int k, MultiIndex<d>& alpha)
188 {
189 pk_multiindex(i,alpha);
190 }
191 };
192
193 template <>
194 struct BasisTraits<BasisType::Qk>
195 {
196 template <int k, int d>
197 struct Size
198 {
199 enum{
200 value = qk_size(k,d)
201 };
202 };
203
204 template <int k, int d>
205 struct Order
206 {
207 enum{
208 // value = d*k
209 value = k
210 };
211 };
212
213 static int size(int k, int d)
214 {
215 return qk_size(k,d);
216 }
217
218 template <int d>
219 static void multiindex(int i, int k, MultiIndex<d>& alpha)
220 {
221 return qk_multiindex(i,k,alpha);
222 }
223 };
224
225 //=====================================================
226 // Integration kernels for monomials
227 //=====================================================
228
230 inline long binomial (long n, long k)
231 {
232 // pick the shorter version of
233 // n*(n-1)*...*(k+1)/((n-k)*(n-k-1)*...*1)
234 // and
235 // n*(n-1)*...*(n-k+1)/(k*(k-1)*...*1)
236 if (2*k>=n)
237 {
238 long nominator=1;
239 for (long i=k+1; i<=n; i++) nominator *= i;
240 long denominator=1;
241 for (long i=2; i<=n-k; i++) denominator *= i;
242 return nominator/denominator;
243 }
244 else
245 {
246 long nominator=1;
247 for (long i=n-k+1; i<=n; i++) nominator *= i;
248 long denominator=1;
249 for (long i=2; i<=k; i++) denominator *= i;
250 return nominator/denominator;
251 }
252 }
253
260 template<typename ComputationFieldType, Dune::GeometryType::BasicType bt, int d>
262 {
263 public:
265 ComputationFieldType integrate (const MultiIndex<d>& a) const
266 {
267 DUNE_THROW(Dune::NotImplemented,"non-specialized version of MonomalIntegrator called. Please implement.");
268 }
269 };
270
273 template<typename ComputationFieldType, int d>
274 class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::cube,d>
275 {
276 public:
278 ComputationFieldType integrate (const MultiIndex<d>& a) const
279 {
280 ComputationFieldType result(1.0);
281 for (int i=0; i<d; i++)
282 {
283 ComputationFieldType exponent(a[i]+1);
284 result = result/exponent;
285 }
286 return result;
287 }
288 };
289
292 template<typename ComputationFieldType>
293 class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,1>
294 {
295 public:
297 ComputationFieldType integrate (const MultiIndex<1>& a) const
298 {
299 ComputationFieldType one(1.0);
300 ComputationFieldType exponent0(a[0]+1);
301 return one/exponent0;
302 }
303 };
304
307 template<typename ComputationFieldType>
308 class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,2>
309 {
310 public:
312 ComputationFieldType integrate (const MultiIndex<2>& a) const
313 {
314 ComputationFieldType sum(0.0);
315 for (int k=0; k<=a[1]+1; k++)
316 {
317 int sign=1;
318 if (k%2==1) sign=-1;
319 ComputationFieldType nom(sign*binomial(a[1]+1,k));
320 ComputationFieldType denom(a[0]+k+1);
321 sum = sum + (nom/denom);
322 }
323 ComputationFieldType denom(a[1]+1);
324 return sum/denom;
325 }
326 };
327
330 template<typename ComputationFieldType>
331 class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,3>
332 {
333 public:
335 ComputationFieldType integrate (const MultiIndex<3>& a) const
336 {
337 ComputationFieldType sumk(0.0);
338 for (int k=0; k<=a[2]+1; k++)
339 {
340 int sign=1;
341 if (k%2==1) sign=-1;
342 ComputationFieldType nom(sign*binomial(a[2]+1,k));
343 ComputationFieldType denom(a[1]+k+1);
344 sumk = sumk + (nom/denom);
345 }
346 ComputationFieldType sumj(0.0);
347 for (int j=0; j<=a[1]+a[2]+2; j++)
348 {
349 int sign=1;
350 if (j%2==1) sign=-1;
351 ComputationFieldType nom(sign*binomial(a[1]+a[2]+2,j));
352 ComputationFieldType denom(a[0]+j+1);
353 sumj = sumj + (nom/denom);
354 }
355 ComputationFieldType denom(a[2]+1);
356 return sumk*sumj/denom;
357 }
358 };
359
360 //=====================================================
361 // construct orthonormal basis
362 //=====================================================
363
365 template<typename F, int d>
367 {
368 template<typename X, typename A>
369 static F compute (const X& x, const A& a)
370 {
371 F prod(1.0);
372 for (int i=0; i<a[d]; i++)
373 prod = prod*x[d];
374 return prod*MonomialEvaluate<F,d-1>::compute(x,a);
375 }
376 };
377
378 template<typename F>
379 struct MonomialEvaluate<F,0>
380 {
381 template<typename X, typename A>
382 static F compute (const X& x, const A& a)
383 {
384 F prod(1.0);
385 for (int i=0; i<a[0]; i++)
386 prod = prod*x[0];
387 return prod;
388 }
389 };
390
420 template<typename FieldType, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
422 {
423 typedef BasisTraits<basisType> Traits;
424 public:
425 enum{ n = Traits::template Size<k,d>::value };
428
429 // construct orthonormal basis
431 : coeffs(new LowprecMat)
432 {
433 for (int i=0; i<d; ++i)
434 gradcoeffs[i].reset(new LowprecMat());
435 // compute index to multiindex map
436 for (int i=0; i<n; i++)
437 {
438 alpha[i].reset(new MultiIndex<d>());
439 Traits::multiindex(i,k,*alpha[i]);
440 //std::cout << "i=" << i << " alpha_i=" << alpha[i] << std::endl;
441 }
442
443 orthonormalize();
444 }
445
446 // construct orthonormal basis from an other basis
447 template<class LFE>
448 OrthonormalPolynomialBasis (const LFE & lfe)
449 : coeffs(new LowprecMat)
450 {
451 for (int i=0; i<d; ++i)
452 gradcoeffs[i].reset(new LowprecMat());
453 // compute index to multiindex map
454 for (int i=0; i<n; i++)
455 {
456 alpha[i].reset(new MultiIndex<d>());
457 Traits::multiindex(i,k,*alpha[i]);
458 //std::cout << "i=" << i << " alpha_i=" << alpha[i] << std::endl;
459 }
460
461 orthonormalize();
462 }
463
464 // return dimension of P_l
465 int size (int l)
466 {
467 return Traits::size(l,d);
468 }
469
470 // evaluate all basis polynomials at given point
471 template<typename Point, typename Result>
472 void evaluateFunction (const Point& x, Result& r) const
473 {
474 std::fill(r.begin(),r.end(),0.0);
475 for (int j=0; j<n; ++j)
476 {
477 const FieldType monomial_value = MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
478 for (int i=j; i<n; ++i)
479 r[i] += (*coeffs)[i][j] * monomial_value;
480 }
481 }
482
483 // evaluate all basis polynomials at given point
484 template<typename Point, typename Result>
485 void evaluateJacobian (const Point& x, Result& r) const
486 {
487 std::fill(r.begin(),r.end(),0.0);
488
489 for (int j=0; j<n; ++j)
490 {
491 const FieldType monomial_value = MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
492 for (int i=j; i<n; ++i)
493 for (int s=0; s<d; ++s)
494 r[i][0][s] += (*gradcoeffs[s])[i][j]*monomial_value;
495 }
496 }
497
498 // evaluate all basis polynomials at given point up to order l <= k
499 template<typename Point, typename Result>
500 void evaluateFunction (int l, const Point& x, Result& r) const
501 {
502 if (l>k)
503 DUNE_THROW(Dune::RangeError,"l>k in OrthonormalPolynomialBasis::evaluateFunction");
504
505 for (int i=0; i<Traits::size(l,d); i++)
506 {
507 FieldType sum(0.0);
508 for (int j=0; j<=i; j++)
509 sum = sum + (*coeffs)[i][j]*MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
510 r[i] = sum;
511 }
512 }
513
514 // evaluate all basis polynomials at given point
515 template<typename Point, typename Result>
516 void evaluateJacobian (int l, const Point& x, Result& r) const
517 {
518 if (l>k)
519 DUNE_THROW(Dune::RangeError,"l>k in OrthonormalPolynomialBasis::evaluateFunction");
520
521 for (int i=0; i<Traits::size(l,d); i++)
522 {
523 FieldType sum[d];
524 for (int s=0; s<d; s++)
525 {
526 sum[s] = 0.0;
527 for (int j=0; j<=i; j++)
528 sum[s] += (*gradcoeffs[s])[i][j]*MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
529 }
530 for (int s=0; s<d; s++) r[i][0][s] = sum[s];
531 }
532 }
533
534 private:
535 // store multiindices and coefficients on heap
536 std::array<std::shared_ptr<MultiIndex<d> >,n> alpha; // store index to multiindex map
537 std::shared_ptr<LowprecMat> coeffs; // coefficients with respect to monomials
538 std::array<std::shared_ptr<LowprecMat>,d > gradcoeffs; // coefficients of gradient
539
540 // compute orthonormalized shapefunctions from a given set of coefficients
541 void orthonormalize()
542 {
543 // run Gram-Schmidt orthogonalization procedure in high precission
544 gram_schmidt();
545
546 // std::cout << "orthogonal basis monomial representation" << std::endl;
547 // for (int i=0; i<n; i++)
548 // {
549 // std::cout << "phi_" << i << ":" ;
550 // for (int j=0; j<=i; j++)
551 // std::cout << " (" << alpha[j] << "," << coeffs[i][j] << ")";
552 // std::cout << std::endl;
553 // }
554
555 // compute coefficients of gradient
556 for (int s=0; s<d; s++)
557 for (int i=0; i<n; i++)
558 for (int j=0; j<=i; j++)
559 (*gradcoeffs[s])[i][j] = 0;
560 for (int i=0; i<n; i++)
561 for (int j=0; j<=i; j++)
562 for (int s=0; s<d; s++)
563 if ((*alpha[j])[s]>0)
564 {
565 MultiIndex<d> beta = *alpha[j]; // get exponents
566 FieldType factor = beta[s];
567 beta[s] -= 1;
568 int l = invert_index(beta);
569 (*gradcoeffs[s])[i][l] += (*coeffs)[i][j]*factor;
570 }
571
572 // for (int s=0; s<d; s++)
573 // {
574 // std::cout << "derivative in direction " << s << std::endl;
575 // for (int i=0; i<n; i++)
576 // {
577 // std::cout << "phi_" << i << ":" ;
578 // for (int j=0; j<=i; j++)
579 // std::cout << " (" << alpha[j] << "," << gradcoeffs[s][i][j] << ")";
580 // std::cout << std::endl;
581 // }
582 // }
583 }
584
585 // get index from a given multiindex
586 int invert_index (MultiIndex<d>& a)
587 {
588 for (int i=0; i<n; i++)
589 {
590 bool found(true);
591 for (int j=0; j<d; j++)
592 if (a[j]!=(*alpha[i])[j]) found=false;
593 if (found) return i;
594 }
595 DUNE_THROW(Dune::RangeError,"index not found in invertindex");
596 }
597
598 void gram_schmidt ()
599 {
600 // allocate a high precission matrix on the heap
601 HighprecMat *p = new HighprecMat();
602 HighprecMat& c = *p;
603
604 // fill identity matrix
605 for (int i=0; i<n; i++)
606 for (int j=0; j<n; j++)
607 if (i==j)
608 c[i][j] = ComputationFieldType(1.0);
609 else
610 c[i][j] = ComputationFieldType(0.0);
611
612 // the Gram-Schmidt loop
614 for (int i=0; i<n; i++)
615 {
616 // store orthogonalization coefficients for scaling
617 ComputationFieldType bi[n];
618
619 // make p_i orthogonal to previous polynomials p_j
620 for (int j=0; j<i; j++)
621 {
622 // p_j is represented with monomials
623 bi[j] = ComputationFieldType(0.0);
624 for (int l=0; l<=j; l++)
625 {
626 MultiIndex<d> a;
627 for (int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
628 bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
629 }
630 for (int l=0; l<=j; l++)
631 c[i][l] = c[i][l] - bi[j]*c[j][l];
632 }
633
634 // scale ith polynomial
635 ComputationFieldType s2(0.0);
636 MultiIndex<d> a;
637 for (int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
638 s2 = s2 + integrator.integrate(a);
639 for (int j=0; j<i; j++)
640 s2 = s2 - bi[j]*bi[j];
641 ComputationFieldType s(1.0);
642 using std::sqrt;
643 s = s/sqrt(s2);
644 for (int l=0; l<=i; l++)
645 c[i][l] = s*c[i][l];
646 }
647
648 // store coefficients in low precission type
649 for (int i=0; i<n; i++)
650 for (int j=0; j<n; j++)
651 (*coeffs)[i][j] = c[i][j];
652
653 delete p;
654
655 //std::cout << coeffs << std::endl;
656 }
657 };
658
659 } // PB namespace
660
661 // define the local finite element here
662
663 template<class D, class R, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
664 class OPBLocalBasis
665 {
666 typedef PB::BasisTraits<basisType> BasisTraits;
668 PolynomialBasis opb;
670
671 public:
673 enum{ n = BasisTraits::template Size<k,d>::value };
674
676
677 OPBLocalBasis (int order_) : opb(), gt(bt,d) {}
678
679 template<class LFE>
680 OPBLocalBasis (int order_, const LFE & lfe) : opb(lfe), gt(bt,d) {}
681
683
684 unsigned int size () const { return n; }
685
687 inline void evaluateFunction (const typename Traits::DomainType& in,
688 std::vector<typename Traits::RangeType>& out) const {
689 out.resize(n);
690 opb.evaluateFunction(in,out);
691 }
692
694 inline void
695 evaluateJacobian (const typename Traits::DomainType& in,
696 std::vector<typename Traits::JacobianType>& out) const {
697 out.resize(n);
698 opb.evaluateJacobian(in,out);
699 }
700
702 void partial(const std::array<unsigned int, Traits::dimDomain>& order,
703 const typename Traits::DomainType& in,
704 std::vector<typename Traits::RangeType>& out) const {
705 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
706 if (totalOrder == 0) {
707 evaluateFunction(in, out);
708 } else {
709 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
710 }
711 }
712
714 unsigned int order () const {
715 return BasisTraits::template Order<k,d>::value;
716 }
717
718 Dune::GeometryType type () const { return gt; }
719 };
720
721 template<int k, int d, PB::BasisType basisType = PB::BasisType::Pk>
722 class OPBLocalCoefficients
723 {
724 enum{ n = PB::BasisTraits<basisType>::template Size<k,d>::value };
725 public:
726 OPBLocalCoefficients (int order_) : li(n) {
727 for (int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
728 }
729
731 std::size_t size () const { return n; }
732
734 const Dune::LocalKey& localKey (int i) const {
735 return li[i];
736 }
737
738 private:
739 std::vector<Dune::LocalKey> li;
740 };
741
742 template<class LB>
743 class OPBLocalInterpolation
744 {
745 const LB& lb;
746
747 public:
748 OPBLocalInterpolation (const LB& lb_, int order_) : lb(lb_) {}
749
751 template<typename F, typename C>
752 void interpolate (const F& ff, std::vector<C>& out) const
753 {
754 // select quadrature rule
755 typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
756
757 typedef typename LB::Traits::RangeType RangeType;
758 const int d = LB::Traits::dimDomain;
760 rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
761
762 auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
763
764 // prepare result
765 out.resize(LB::n);
766 for (int i=0; i<LB::n; i++) out[i] = 0.0;
767
768 // loop over quadrature points
770 it=rule.begin(); it!=rule.end(); ++it)
771 {
772 // evaluate function at quadrature point
773 typename LB::Traits::DomainType x;
774 RangeType y;
775 for (int i=0; i<d; i++) x[i] = it->position()[i];
776 y = f(x);
777
778 // evaluate the basis
779 std::vector<RangeType> phi(LB::n);
780 lb.evaluateFunction(it->position(),phi);
781
782 // do integration
783 for (int i=0; i<LB::n; i++)
784 out[i] += y*phi[i]*it->weight();
785 }
786 }
787 };
788
789 template<class D, class R, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=Dune::PB::DefaultComputationalFieldType, PB::BasisType basisType = PB::BasisType::Pk>
790 class OPBLocalFiniteElement
791 {
793 OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> basis;
794 OPBLocalCoefficients<k,d,basisType> coefficients;
795 OPBLocalInterpolation<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> > interpolation;
796 public:
798 OPBLocalCoefficients<k,d,basisType>,
799 OPBLocalInterpolation<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType> > > Traits;
800
802
803 OPBLocalFiniteElement ()
804 : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
805 {}
806
807 template<class LFE>
808 explicit OPBLocalFiniteElement (const LFE & lfe)
809 : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
810 {}
811
813
814 OPBLocalFiniteElement (const OPBLocalFiniteElement & other)
815 : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
816 {}
817
818 const typename Traits::LocalBasisType& localBasis () const
819 {
820 return basis;
821 }
822
823 const typename Traits::LocalCoefficientsType& localCoefficients () const
824 {
825 return coefficients;
826 }
827
828 const typename Traits::LocalInterpolationType& localInterpolation () const
829 {
830 return interpolation;
831 }
832
835 std::size_t size() const
836 {
837 return basis.size();
838 }
839
840 Dune::GeometryType type () const { return gt; }
841
842 OPBLocalFiniteElement* clone () const {
843 return new OPBLocalFiniteElement(*this);
844 }
845 };
846
847}
848
851#endif // DUNE_PDELAB_FINITEELEMENT_L2ORTHONORMAL_HH
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:122
@ simplex
Simplicial element in any nonnegative dimension.
Definition: type.hh:121
Describe position of one degree of freedom.
Definition: localkey.hh:24
Default exception for dummy implementations.
Definition: exceptions.hh:263
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:278
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:297
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:312
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:335
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:262
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:265
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:422
Definition: polynomialbasis.hh:65
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: polynomialbasis.hh:119
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: polynomialbasis.hh:127
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
Default exception class for range errors.
Definition: exceptions.hh:254
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Implements a vector constructed from a given type representing a field and a compile-time given size.
Wrapper for the GNU multiprecision (GMP) library.
#define DUNE_NO_DEPRECATED_END
Ignore deprecation warnings (end)
Definition: deprecated.hh:38
#define DUNE_NO_DEPRECATED_BEGIN
Ignore deprecation warnings (start)
Definition: deprecated.hh:32
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
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
static constexpr T binomial(const T &n, const T &k) noexcept
calculate the binomial coefficient n over k as a constexpr
Definition: math.hh:131
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
STL namespace.
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
compute
Definition: l2orthonormal.hh:367
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 (Nov 23, 23:29, 2024)