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