dune-pdelab  2.0.0
l2orthonormal.hh
Go to the documentation of this file.
1 #ifndef DUNE_PM_ORTHONORMAL_HH
2 #define DUNE_PM_ORTHONORMAL_HH
3 
4 #include<iostream>
5 #include<algorithm>
6 #include<dune/common/fvector.hh>
7 #include<dune/common/fmatrix.hh>
8 #include<dune/common/exceptions.hh>
9 #include<dune/common/fvector.hh>
10 #include<dune/common/gmpfield.hh>
11 #include<dune/common/shared_ptr.hh>
12 #include<dune/common/array.hh>
13 
14 #include<dune/geometry/referenceelements.hh>
15 #include<dune/geometry/quadraturerules.hh>
16 #include<dune/geometry/type.hh>
17 
18 #include<dune/localfunctions/common/localbasis.hh>
19 #include<dune/localfunctions/common/localkey.hh>
20 #include<dune/localfunctions/common/localfiniteelementtraits.hh>
21 
32 namespace Dune {
33  namespace PB {
34 
35  //=====================================================
36  // TMPs for computing number of basis functions in P_k
37  //=====================================================
38 
39  template<int k, int d> struct PkSize;
40 
41  template<int j, int k, int d>
42  struct PkSizeHelper
43  {
44  enum{
45  value = PkSizeHelper<j-1,k,d>::value + PkSize<k-j,d-1>::value
46  };
47  };
48 
49  template<int k, int d>
50  struct PkSizeHelper<0,k,d>
51  {
52  enum{
54  };
55  };
56 
57  // This is the main class
58  // k is the polynomial degree and d is the space dimension
59  // then PkSize<k,d> is the number of polynomials of at most total degree k.
60  template<int k, int d>
61  struct PkSize
62  {
63  enum{
65  };
66  };
67 
68  template<>
69  struct PkSize<0,1>
70  {
71  enum{
73  };
74  };
75 
76  template<int k>
77  struct PkSize<k,1>
78  {
79  enum{
80  value=k+1
81  };
82  };
83 
84  template<int d>
85  struct PkSize<0,d>
86  {
87  enum{
89  };
90  };
91 
92  // number of polynomials of exactly degree k
93  template<int k, int d>
94  struct PkExactSize
95  {
96  enum{
98  };
99  };
100 
101  template<int d>
102  struct PkExactSize<0,d>
103  {
104  enum{
106  };
107  };
108 
109  //=====================================================
110  // TMPs for computing number of basis functions in Q_k
111  //=====================================================
112 
113  // This is the main class
114  // usage: QkSize<2,3>::value
115  // k is the polynomial degree, d is the space dimension
116  template<int k, int d>
117  struct QkSize
118  {
119  enum{
121  };
122  };
123 
124  template<>
125  struct QkSize<0,1>
126  {
127  enum{
129  };
130  };
131 
132  template<int k>
133  struct QkSize<k,1>
134  {
135  enum{
136  value=k+1
137  };
138  };
139 
140  template<int d>
141  struct QkSize<0,d>
142  {
143  enum{
145  };
146  };
147 
148  //=====================================================
149  // Type to represent a multiindex in d dimensions
150  //=====================================================
151 
152  template<int d>
153  class MultiIndex : public Dune::array<int,d>
154  {
155  public:
156 
157  MultiIndex () : Dune::array<int,d>()
158  {
159  }
160 
161  };
162 
163  // the number of polynomials of at most degree k in space dimension d (as run-time function)
164  inline int pk_size (int k, int d)
165  {
166  if (k==0) return 1;
167  if (d==1) return k+1;
168  int n=0;
169  for (int j=0; j<=k; j++)
170  n += pk_size(k-j,d-1);
171  return n;
172  }
173 
174  // the number of polynomials of exactly degree k in space dimension d (as run-time function)
175  inline int pk_size_exact (int k, int d)
176  {
177  if (k==0)
178  return 1;
179  else
180  return pk_size(k,d)-pk_size(k-1,d);
181  }
182 
183  // enumerate all multiindices of degree k and find the i'th
184  template<int d>
185  void pk_enumerate_multiindex (MultiIndex<d>& alpha, int& count, int norm, int dim, int k, int i)
186  {
187  // check if dim is valid
188  if (dim>=d) return;
189 
190  // put all k to current dimension and check
191  alpha[dim]=k; count++; // alpha has correct norm
192  // std::cout << "dadi alpha=" << alpha << " count=" << count << " norm=" << norm+k << " dim=" << dim << " k=" << k << " i=" << i << std::endl;
193  if (count==i) return; // found the index
194 
195  // search recursively
196  for (int m=k-1; m>=0; m--)
197  {
198  alpha[dim]=m;
199  //std::cout << "dada alpha=" << alpha << " count=" << count << " norm=" << norm+m << " dim=" << dim << " k=" << k << " i=" << i << std::endl;
200  pk_enumerate_multiindex(alpha,count,norm+m,dim+1,k-m,i);
201  if (count==i) return;
202  }
203 
204  // reset if index is not yet found
205  alpha[dim]=0;
206  }
207 
208  // map integer 0<=i<pk_size(k,d) to multiindex
209  template<int d>
210  void pk_multiindex (int i, MultiIndex<d>& alpha)
211  {
212  for (int j=0; j<d; j++) alpha[j] = 0; // set alpha to 0
213  if (i==0) return; // handle case i==0 and k==0
214  for (int k=1; k<25; k++)
215  if (i>=pk_size(k-1,d) && i<pk_size(k,d))
216  {
217  int count = -1;
218  pk_enumerate_multiindex(alpha,count,0,0,k,i-pk_size(k-1,d));
219  return;
220  }
221  DUNE_THROW(Dune::NotImplemented,"Polynomial degree greater than 24 in pk_multiindex");
222  }
223 
224  // the number of polynomials of at most degree k in space dimension d (as run-time function)
225  inline int qk_size (int k, int d)
226  {
227  int n = 1;
228  for (int i=0; i<d; ++i)
229  n *= (k+1);
230  return n;
231  }
232 
233  // map integer 0<=i<qk_size(k,d) to multiindex
234  template<int d>
235  void qk_multiindex (int i, int k, MultiIndex<d>& alpha)
236  {
237  for (int j = 0; j<d; ++j) {
238  alpha[j] = i % (k+1);
239  i /= (k+1);
240  }
241  }
242 
243  //=====================================================
244  // Traits classes to group Pk and Qk specifics
245  //=====================================================
246  enum BasisType {
248  };
249 
250  template <BasisType basisType>
251  struct BasisTraits;
252 
253  template <>
255  {
256  template <int k, int d>
257  struct Size
258  {
259  enum{
261  };
262  };
263 
264  template <int k, int d>
265  struct Order
266  {
267  enum{
268  value = k
269  };
270  };
271 
272  static int size(int k, int d)
273  {
274  return pk_size(k,d);
275  }
276 
277  template <int d>
278  static void multiindex(int i, int k, MultiIndex<d>& alpha)
279  {
280  pk_multiindex(i,alpha);
281  }
282  };
283 
284  template <>
286  {
287  template <int k, int d>
288  struct Size
289  {
290  enum{
292  };
293  };
294 
295  template <int k, int d>
296  struct Order
297  {
298  enum{
299  // value = d*k
300  value = k
301  };
302  };
303 
304  static int size(int k, int d)
305  {
306  return qk_size(k,d);
307  }
308 
309  template <int d>
310  static void multiindex(int i, int k, MultiIndex<d>& alpha)
311  {
312  return qk_multiindex(i,k,alpha);
313  }
314  };
315 
316  //=====================================================
317  // Integration kernels for monomials
318  //=====================================================
319 
321  inline long binomial (long n, long k)
322  {
323  // pick the shorter version of
324  // n*(n-1)*...*(k+1)/((n-k)*(n-k-1)*...*1)
325  // and
326  // n*(n-1)*...*(n-k+1)/(k*(k-1)*...*1)
327  if (2*k>=n)
328  {
329  long nominator=1;
330  for (long i=k+1; i<=n; i++) nominator *= i;
331  long denominator=1;
332  for (long i=2; i<=n-k; i++) denominator *= i;
333  return nominator/denominator;
334  }
335  else
336  {
337  long nominator=1;
338  for (long i=n-k+1; i<=n; i++) nominator *= i;
339  long denominator=1;
340  for (long i=2; i<=k; i++) denominator *= i;
341  return nominator/denominator;
342  }
343  }
344 
351  template<typename ComputationFieldType, Dune::GeometryType::BasicType bt, int d>
353  {
354  public:
356  ComputationFieldType integrate (const MultiIndex<d>& a) const
357  {
358  DUNE_THROW(Dune::NotImplemented,"non-specialized version of MonomalIntegrator called. Please implement.");
359  }
360  };
361 
364  template<typename ComputationFieldType, int d>
365  class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::cube,d>
366  {
367  public:
369  ComputationFieldType integrate (const MultiIndex<d>& a) const
370  {
371  ComputationFieldType result(1.0);
372  for (int i=0; i<d; i++)
373  {
374  ComputationFieldType exponent(a[i]+1);
375  result = result/exponent;
376  }
377  return result;
378  }
379  };
380 
383  template<typename ComputationFieldType>
384  class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,1>
385  {
386  public:
388  ComputationFieldType integrate (const MultiIndex<1>& a) const
389  {
390  ComputationFieldType one(1.0);
391  ComputationFieldType exponent0(a[0]+1);
392  return one/exponent0;
393  }
394  };
395 
398  template<typename ComputationFieldType>
399  class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,2>
400  {
401  public:
403  ComputationFieldType integrate (const MultiIndex<2>& a) const
404  {
405  ComputationFieldType sum(0.0);
406  for (int k=0; k<=a[1]+1; k++)
407  {
408  int sign=1;
409  if (k%2==1) sign=-1;
410  ComputationFieldType nom(sign*binomial(a[1]+1,k));
411  ComputationFieldType denom(a[0]+k+1);
412  sum = sum + (nom/denom);
413  }
414  ComputationFieldType denom(a[1]+1);
415  return sum/denom;
416  }
417  };
418 
421  template<typename ComputationFieldType>
422  class MonomialIntegrator<ComputationFieldType,Dune::GeometryType::simplex,3>
423  {
424  public:
426  ComputationFieldType integrate (const MultiIndex<3>& a) const
427  {
428  ComputationFieldType sumk(0.0);
429  for (int k=0; k<=a[2]+1; k++)
430  {
431  int sign=1;
432  if (k%2==1) sign=-1;
433  ComputationFieldType nom(sign*binomial(a[2]+1,k));
434  ComputationFieldType denom(a[1]+k+1);
435  sumk = sumk + (nom/denom);
436  }
437  ComputationFieldType sumj(0.0);
438  for (int j=0; j<=a[1]+a[2]+2; j++)
439  {
440  int sign=1;
441  if (j%2==1) sign=-1;
442  ComputationFieldType nom(sign*binomial(a[1]+a[2]+2,j));
443  ComputationFieldType denom(a[0]+j+1);
444  sumj = sumj + (nom/denom);
445  }
446  ComputationFieldType denom(a[2]+1);
447  return sumk*sumj/denom;
448  }
449  };
450 
451  //=====================================================
452  // construct orthonormal basis
453  //=====================================================
454 
456  template<typename F, int d>
458  {
459  template<typename X, typename A>
460  static F compute (const X& x, const A& a)
461  {
462  F prod(1.0);
463  for (int i=0; i<a[d]; i++)
464  prod = prod*x[d];
465  return prod*MonomialEvaluate<F,d-1>::compute(x,a);
466  }
467  };
468 
469  template<typename F>
470  struct MonomialEvaluate<F,0>
471  {
472  template<typename X, typename A>
473  static F compute (const X& x, const A& a)
474  {
475  F prod(1.0);
476  for (int i=0; i<a[0]; i++)
477  prod = prod*x[0];
478  return prod;
479  }
480  };
481 
511  template<typename FieldType, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=FieldType, BasisType basisType = BasisType::Pk>
513  {
515  public:
516  enum{ n = Traits::template Size<k,d>::value };
517  typedef Dune::FieldMatrix<FieldType,n,n> LowprecMat;
518  typedef Dune::FieldMatrix<ComputationFieldType,n,n> HighprecMat;
519 
520  // construct orthonormal basis
522  : coeffs(new LowprecMat)
523  {
524  for (int i=0; i<d; ++i)
525  gradcoeffs[i].reset(new LowprecMat());
526  // compute index to multiindex map
527  for (int i=0; i<n; i++)
528  {
529  alpha[i].reset(new MultiIndex<d>());
530  Traits::multiindex(i,k,*alpha[i]);
531  //std::cout << "i=" << i << " alpha_i=" << alpha[i] << std::endl;
532  }
533 
534  orthonormalize();
535  }
536 
537  // construct orthonormal basis from an other basis
538  template<class LFE>
539  OrthonormalPolynomialBasis (const LFE & lfe)
540  : coeffs(new LowprecMat)
541  {
542  for (int i=0; i<d; ++i)
543  gradcoeffs[i].reset(new LowprecMat());
544  // compute index to multiindex map
545  for (int i=0; i<n; i++)
546  {
547  alpha[i].reset(new MultiIndex<d>());
548  Traits::multiindex(i,k,*alpha[i]);
549  //std::cout << "i=" << i << " alpha_i=" << alpha[i] << std::endl;
550  }
551 
552  orthonormalize();
553  }
554 
555  // return dimension of P_l
556  int size (int l)
557  {
558  return Traits::size(l,d);
559  }
560 
561  // evaluate all basis polynomials at given point
562  template<typename Point, typename Result>
563  void evaluateFunction (const Point& x, Result& r) const
564  {
565  std::fill(r.begin(),r.end(),0.0);
566  for (int j=0; j<n; ++j)
567  {
568  const FieldType monomial_value = MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
569  for (int i=j; i<n; ++i)
570  r[i] += (*coeffs)[i][j] * monomial_value;
571  }
572  }
573 
574  // evaluate all basis polynomials at given point
575  template<typename Point, typename Result>
576  void evaluateJacobian (const Point& x, Result& r) const
577  {
578  std::fill(r.begin(),r.end(),0.0);
579 
580  for (int j=0; j<n; ++j)
581  {
582  const FieldType monomial_value = MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
583  for (int i=j; i<n; ++i)
584  for (int s=0; s<d; ++s)
585  r[i][0][s] += (*gradcoeffs[s])[i][j]*monomial_value;
586  }
587  }
588 
589  // evaluate all basis polynomials at given point up to order l <= k
590  template<typename Point, typename Result>
591  void evaluateFunction (int l, const Point& x, Result& r) const
592  {
593  if (l>k)
594  DUNE_THROW(Dune::RangeError,"l>k in OrthonormalPolynomialBasis::evaluateFunction");
595 
596  for (int i=0; i<Traits::size(l,d); i++)
597  {
598  FieldType sum(0.0);
599  for (int j=0; j<=i; j++)
600  sum = sum + (*coeffs)[i][j]*MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
601  r[i] = sum;
602  }
603  }
604 
605  // evaluate all basis polynomials at given point
606  template<typename Point, typename Result>
607  void evaluateJacobian (int l, const Point& x, Result& r) const
608  {
609  if (l>k)
610  DUNE_THROW(Dune::RangeError,"l>k in OrthonormalPolynomialBasis::evaluateFunction");
611 
612  for (int i=0; i<Traits::size(l,d); i++)
613  {
614  FieldType sum[d];
615  for (int s=0; s<d; s++)
616  {
617  sum[s] = 0.0;
618  for (int j=0; j<=i; j++)
619  sum[s] += (*gradcoeffs[s])[i][j]*MonomialEvaluate<FieldType,d-1>::compute(x,*alpha[j]);
620  }
621  for (int s=0; s<d; s++) r[i][0][s] = sum[s];
622  }
623  }
624 
625  private:
626  // store multiindices and coefficients on heap
627  Dune::array<Dune::shared_ptr<MultiIndex<d> >,n> alpha; // store index to multiindex map
628  Dune::shared_ptr<LowprecMat> coeffs; // coefficients with respect to monomials
629  Dune::array<Dune::shared_ptr<LowprecMat>,d > gradcoeffs; // coefficients of gradient
630 
631  // compute orthonormalized shapefunctions from a given set of coefficients
632  void orthonormalize()
633  {
634  // run Gram-Schmidt orthogonalization procedure in high precission
635  gram_schmidt();
636 
637  // std::cout << "orthogonal basis monomial representation" << std::endl;
638  // for (int i=0; i<n; i++)
639  // {
640  // std::cout << "phi_" << i << ":" ;
641  // for (int j=0; j<=i; j++)
642  // std::cout << " (" << alpha[j] << "," << coeffs[i][j] << ")";
643  // std::cout << std::endl;
644  // }
645 
646  // compute coefficients of gradient
647  for (int s=0; s<d; s++)
648  for (int i=0; i<n; i++)
649  for (int j=0; j<=i; j++)
650  (*gradcoeffs[s])[i][j] = 0;
651  for (int i=0; i<n; i++)
652  for (int j=0; j<=i; j++)
653  for (int s=0; s<d; s++)
654  if ((*alpha[j])[s]>0)
655  {
656  MultiIndex<d> beta = *alpha[j]; // get exponents
657  FieldType factor = beta[s];
658  beta[s] -= 1;
659  int l = invert_index(beta);
660  (*gradcoeffs[s])[i][l] += (*coeffs)[i][j]*factor;
661  }
662 
663  // for (int s=0; s<d; s++)
664  // {
665  // std::cout << "derivative in direction " << s << std::endl;
666  // for (int i=0; i<n; i++)
667  // {
668  // std::cout << "phi_" << i << ":" ;
669  // for (int j=0; j<=i; j++)
670  // std::cout << " (" << alpha[j] << "," << gradcoeffs[s][i][j] << ")";
671  // std::cout << std::endl;
672  // }
673  // }
674  }
675 
676  // get index from a given multiindex
677  int invert_index (MultiIndex<d>& a)
678  {
679  for (int i=0; i<n; i++)
680  {
681  bool found(true);
682  for (int j=0; j<d; j++)
683  if (a[j]!=(*alpha[i])[j]) found=false;
684  if (found) return i;
685  }
686  DUNE_THROW(Dune::RangeError,"index not found in invertindex");
687  }
688 
689  void gram_schmidt ()
690  {
691  // allocate a high precission matrix on the heap
692  HighprecMat *p = new HighprecMat();
693  HighprecMat& c = *p;
694 
695  // fill identity matrix
696  for (int i=0; i<n; i++)
697  for (int j=0; j<n; j++)
698  if (i==j)
699  c[i][j] = ComputationFieldType(1.0);
700  else
701  c[i][j] = ComputationFieldType(0.0);
702 
703  // the Gram-Schmidt loop
704  MonomialIntegrator<ComputationFieldType,bt,d> integrator;
705  for (int i=0; i<n; i++)
706  {
707  // store orthogonalization coefficients for scaling
708  ComputationFieldType bi[n];
709 
710  // make p_i orthogonal to previous polynomials p_j
711  for (int j=0; j<i; j++)
712  {
713  // p_j is represented with monomials
714  bi[j] = ComputationFieldType(0.0);
715  for (int l=0; l<=j; l++)
716  {
717  MultiIndex<d> a;
718  for (int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[l])[m];
719  bi[j] = bi[j] + c[j][l]*integrator.integrate(a);
720  }
721  for (int l=0; l<=j; l++)
722  c[i][l] = c[i][l] - bi[j]*c[j][l];
723  }
724 
725  // scale ith polynomial
726  ComputationFieldType s2(0.0);
727  MultiIndex<d> a;
728  for (int m=0; m<d; m++) a[m] = (*alpha[i])[m] + (*alpha[i])[m];
729  s2 = s2 + integrator.integrate(a);
730  for (int j=0; j<i; j++)
731  s2 = s2 - bi[j]*bi[j];
732  ComputationFieldType s(1.0);
733  s = s/std::sqrt(s2);
734  for (int l=0; l<=i; l++)
735  c[i][l] = s*c[i][l];
736  }
737 
738  // store coefficients in low precission type
739  for (int i=0; i<n; i++)
740  for (int j=0; j<n; j++)
741  (*coeffs)[i][j] = c[i][j];
742 
743  delete p;
744 
745  //std::cout << coeffs << std::endl;
746  }
747  };
748 
749  } // PB namespace
750 
751  // define the local finite element here
752 
753  template<class D, class R, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=R, PB::BasisType basisType = PB::BasisType::Pk>
755  {
758  PolynomialBasis opb;
759  Dune::GeometryType gt;
760 
761  public:
762  typedef Dune::LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d>, 0> Traits;
763  enum{ n = BasisTraits::template Size<k,d>::value };
764 
765  OPBLocalBasis (int order_) : opb(), gt(bt,d) {}
766 
767  template<class LFE>
768  OPBLocalBasis (int order_, const LFE & lfe) : opb(lfe), gt(bt,d) {}
769 
770  unsigned int size () const { return n; }
771 
773  inline void evaluateFunction (const typename Traits::DomainType& in,
774  std::vector<typename Traits::RangeType>& out) const {
775  out.resize(n);
776  opb.evaluateFunction(in,out);
777  }
778 
780  inline void
781  evaluateJacobian (const typename Traits::DomainType& in,
782  std::vector<typename Traits::JacobianType>& out) const {
783  out.resize(n);
784  opb.evaluateJacobian(in,out);
785  }
786 
788  unsigned int order () const {
789  return BasisTraits::template Order<k,d>::value;
790  }
791 
792  Dune::GeometryType type () const { return gt; }
793  };
794 
795  template<int k, int d, PB::BasisType basisType = PB::BasisType::Pk>
797  {
799  public:
800  OPBLocalCoefficients (int order_) : li(n) {
801  for (int i=0; i<n; i++) li[i] = Dune::LocalKey(0,0,i);
802  }
803 
805  std::size_t size () const { return n; }
806 
808  const Dune::LocalKey& localKey (int i) const {
809  return li[i];
810  }
811 
812  private:
813  std::vector<Dune::LocalKey> li;
814  };
815 
816  template<class LB>
818  {
819  const LB& lb;
820 
821  public:
822  OPBLocalInterpolation (const LB& lb_, int order_) : lb(lb_) {}
823 
825  template<typename F, typename C>
826  void interpolate (const F& f, std::vector<C>& out) const
827  {
828  // select quadrature rule
829  typedef typename FieldTraits<typename LB::Traits::RangeFieldType>::real_type RealFieldType;
830 
831  typedef typename LB::Traits::RangeType RangeType;
832  const int d = LB::Traits::dimDomain;
833  const Dune::QuadratureRule<RealFieldType,d>&
834  rule = Dune::QuadratureRules<RealFieldType,d>::rule(lb.type(),2*lb.order());
835 
836  // prepare result
837  out.resize(LB::n);
838  for (int i=0; i<LB::n; i++) out[i] = 0.0;
839 
840  // loop over quadrature points
841  for (typename Dune::QuadratureRule<RealFieldType,d>::const_iterator
842  it=rule.begin(); it!=rule.end(); ++it)
843  {
844  // evaluate function at quadrature point
845  typename LB::Traits::DomainType x;
846  RangeType y;
847  for (int i=0; i<d; i++) x[i] = it->position()[i];
848  f.evaluate(x,y);
849 
850  // evaluate the basis
851  std::vector<RangeType> phi(LB::n);
852  lb.evaluateFunction(it->position(),phi);
853 
854  // do integration
855  for (int i=0; i<LB::n; i++)
856  out[i] += y*phi[i]*it->weight();
857  }
858  }
859  };
860 
861  template<class D, class R, int k, int d, Dune::GeometryType::BasicType bt, typename ComputationFieldType=R, PB::BasisType basisType = PB::BasisType::Pk>
863  {
864  Dune::GeometryType gt;
868  public:
869  typedef Dune::LocalFiniteElementTraits<OPBLocalBasis<D,R,k,d,bt,ComputationFieldType,basisType>,
872 
874  : gt(bt,d), basis(k), coefficients(k), interpolation(basis,k)
875  {}
876 
877  template<class LFE>
878  explicit OPBLocalFiniteElement (const LFE & lfe)
879  : gt(bt,d), basis(k, lfe), coefficients(k), interpolation(basis,k)
880  {}
881 
883  : gt(other.gt), basis(other.basis), coefficients(other.coefficients), interpolation(basis,k)
884  {}
885 
886  const typename Traits::LocalBasisType& localBasis () const
887  {
888  return basis;
889  }
890 
891  const typename Traits::LocalCoefficientsType& localCoefficients () const
892  {
893  return coefficients;
894  }
895 
896  const typename Traits::LocalInterpolationType& localInterpolation () const
897  {
898  return interpolation;
899  }
900 
901  Dune::GeometryType type () const { return gt; }
902 
904  return new OPBLocalFiniteElement(*this);
905  }
906  };
907 
908 }
909 
912 #endif
Definition: l2orthonormal.hh:64
Dune::LocalFiniteElementTraits< OPBLocalBasis< D, R, k, d, bt, ComputationFieldType, basisType >, OPBLocalCoefficients< k, d, basisType >, OPBLocalInterpolation< OPBLocalBasis< D, R, k, d, bt, ComputationFieldType, basisType > > > Traits
Definition: l2orthonormal.hh:871
ComputationFieldType integrate(const MultiIndex< 3 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:426
OrthonormalPolynomialBasis(const LFE &lfe)
Definition: l2orthonormal.hh:539
ComputationFieldType integrate(const MultiIndex< 1 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:388
Definition: l2orthonormal.hh:94
std::size_t size() const
number of coefficients
Definition: l2orthonormal.hh:805
Definition: l2orthonormal.hh:251
OPBLocalBasis(int order_)
Definition: l2orthonormal.hh:765
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:460
Definition: l2orthonormal.hh:120
void evaluateJacobian(const Point &x, Result &r) const
Definition: l2orthonormal.hh:576
Definition: l2orthonormal.hh:247
Dune::GeometryType type() const
Definition: l2orthonormal.hh:901
const Traits::LocalInterpolationType & localInterpolation() const
Definition: l2orthonormal.hh:896
const Traits::LocalBasisType & localBasis() const
Definition: l2orthonormal.hh:886
Definition: l2orthonormal.hh:39
OPBLocalFiniteElement()
Definition: l2orthonormal.hh:873
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:512
OPBLocalFiniteElement(const OPBLocalFiniteElement &other)
Definition: l2orthonormal.hh:882
void evaluateJacobian(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:607
Definition: l2orthonormal.hh:45
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:310
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: l2orthonormal.hh:891
Definition: l2orthonormal.hh:862
void evaluateFunction(const Point &x, Result &r) const
Definition: l2orthonormal.hh:563
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: l2orthonormal.hh:773
unsigned int size() const
Definition: l2orthonormal.hh:770
int pk_size_exact(int k, int d)
Definition: l2orthonormal.hh:175
static const int dim
Definition: adaptivity.hh:82
OPBLocalInterpolation(const LB &lb_, int order_)
Definition: l2orthonormal.hh:822
OPBLocalFiniteElement * clone() const
Definition: l2orthonormal.hh:903
static F compute(const X &x, const A &a)
Definition: l2orthonormal.hh:473
OPBLocalBasis(int order_, const LFE &lfe)
Definition: l2orthonormal.hh:768
unsigned int order() const
Polynomial order of the shape functions.
Definition: l2orthonormal.hh:788
static int size(int k, int d)
Definition: l2orthonormal.hh:304
Dune::FieldMatrix< FieldType, n, n > LowprecMat
Definition: l2orthonormal.hh:517
Definition: l2orthonormal.hh:153
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:369
int size(int l)
Definition: l2orthonormal.hh:556
Dune::GeometryType type() const
Definition: l2orthonormal.hh:792
static int size(int k, int d)
Definition: l2orthonormal.hh:272
Definition: l2orthonormal.hh:817
Definition: l2orthonormal.hh:42
OPBLocalCoefficients(int order_)
Definition: l2orthonormal.hh:800
void pk_multiindex(int i, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:210
OrthonormalPolynomialBasis()
Definition: l2orthonormal.hh:521
Dune::FieldMatrix< ComputationFieldType, n, n > HighprecMat
Definition: l2orthonormal.hh:518
Definition: l2orthonormal.hh:796
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdg.hh:125
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
Definition: l2orthonormal.hh:247
const F & f
Definition: common/constraints.hh:145
ComputationFieldType integrate(const MultiIndex< 2 > &a) const
integrate one monomial
Definition: l2orthonormal.hh:403
void pk_enumerate_multiindex(MultiIndex< d > &alpha, int &count, int norm, int dim, int k, int i)
Definition: l2orthonormal.hh:185
MultiIndex()
Definition: l2orthonormal.hh:157
static void multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:278
int qk_size(int k, int d)
Definition: l2orthonormal.hh:225
OPBLocalFiniteElement(const LFE &lfe)
Definition: l2orthonormal.hh:878
Dune::LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d >, 0 > Traits
Definition: l2orthonormal.hh:762
BasisType
Definition: l2orthonormal.hh:246
long binomial(long n, long k)
compute binomial coefficient "n over k"
Definition: l2orthonormal.hh:321
Definition: l2orthonormal.hh:754
Definition: l2orthonormal.hh:97
R p(int i, D x)
Definition: qkdg.hh:62
ComputationFieldType integrate(const MultiIndex< d > &a) const
integrate one monomial
Definition: l2orthonormal.hh:356
void qk_multiindex(int i, int k, MultiIndex< d > &alpha)
Definition: l2orthonormal.hh:235
int pk_size(int k, int d)
Definition: l2orthonormal.hh:164
const Dune::LocalKey & localKey(int i) const
map index i to local key
Definition: l2orthonormal.hh:808
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: l2orthonormal.hh:781
Integrate monomials over the reference element.
Definition: l2orthonormal.hh:352
void evaluateFunction(int l, const Point &x, Result &r) const
Definition: l2orthonormal.hh:591
Definition: l2orthonormal.hh:117
compute
Definition: l2orthonormal.hh:457
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: l2orthonormal.hh:826
const std::string s
Definition: function.hh:1103