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