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