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