dune-pdelab  2.4.1
dglegendre.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 // DG tensor product basis with Legendre polynomials
5 
6 #ifndef DUNE_PDELAB_FINITEELEMENT_DGLEGENDRE_HH
7 #define DUNE_PDELAB_FINITEELEMENT_DGLEGENDRE_HH
8 
9 #include <vector>
10 
11 #include <dune/common/fvector.hh>
12 #include <dune/common/deprecated.hh>
13 
14 #include <dune/geometry/type.hh>
15 #include <dune/geometry/quadraturerules.hh>
16 
17 #include <dune/localfunctions/common/localbasis.hh>
18 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
19 #include <dune/localfunctions/common/localkey.hh>
20 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
21 
22 namespace Dune
23 {
24 
25  namespace LegendreStuff
26  {
27  // This is the size class
28  // k is the polynomial degree,
29  // n is the space dimension
30  template<int k, int n>
31  struct LegendreSize
32  {
33  enum{
35  };
36  };
37 
38  template<>
39  struct LegendreSize<0,1>
40  {
41  enum{
43  };
44  };
45 
46  template<int k>
47  struct LegendreSize<k,1>
48  {
49  enum{
50  value=k+1
51  };
52  };
53 
54  template<int n>
55  struct LegendreSize<0,n>
56  {
57  enum{
59  };
60  };
61 
62  template<int k, int d>
63  Dune::FieldVector<int,d> multiindex (int i)
64  {
65  Dune::FieldVector<int,d> alpha;
66  for (int j=0; j<d; j++)
67  {
68  alpha[j] = i % (k+1);
69  i = i/(k+1);
70  }
71  return alpha;
72  }
73 
75  template<class D, class R, int k>
77  {
78  public:
80  void p (D x, std::vector<R>& value) const
81  {
82  value.resize(k+1);
83  value[0] = 1;
84  value[1] = 2*x-1;
85  for (int n=2; n<=k; n++)
86  value[n] = ((2*n-1)*(2*x-1)*value[n-1]-(n-1)*value[n-2])/n;
87  }
88 
89  // ith Lagrange polynomial of degree k in one dimension
90  R p (int i, D x) const
91  {
92  std::vector<R> v(k+1);
93  p(x,v);
94  return v[i];
95  }
96 
97  // derivative of all polynomials
98  void dp (D x, std::vector<R>& derivative) const
99  {
100  R value[k+1];
101  derivative.resize(k+1);
102  value[0] = 1;
103  derivative[0] = 0.0;
104  value[1] = 2*x-1;
105  derivative[1] = 2.0;
106  for (int n=2; n<=k; n++)
107  {
108  value[n] = ((2*n-1)*(2*x-1)*value[n-1]-(n-1)*value[n-2])/n;
109  derivative[n] = (2*x-1)*derivative[n-1] + 2*n*value[n-1];
110  }
111  }
112 
113  // value and derivative of all polynomials
114  void pdp (D x, std::vector<R>& value, std::vector<R>& derivative) const
115  {
116  value.resize(k+1);
117  derivative.resize(k+1);
118  value[0] = 1;
119  derivative[0] = 0.0;
120  value[1] = 2*x-1;
121  derivative[1] = 2.0;
122  for (int n=2; n<=k; n++)
123  {
124  value[n] = ((2*n-1)*(2*x-1)*value[n-1]-(n-1)*value[n-2])/n;
125  derivative[n] = (2*x-1)*derivative[n-1] + 2*n*value[n-1];
126  }
127  }
128 
129  // derivative of ith Lagrange polynomial of degree k in one dimension
130  R dp (int i, D x) const
131  {
132  std::vector<R> v(k+1);
133  dp(x,v);
134  return v[i];
135  }
136  };
137 
138  template<class D, class R>
140  {
141  public:
143  void p (D x, std::vector<R>& value) const
144  {
145  value.resize(1);
146  value[0] = 1.0;
147  }
148 
149  // ith Lagrange polynomial of degree k in one dimension
150  R p (int i, D x) const
151  {
152  return 1.0;
153  }
154 
155  // derivative of all polynomials
156  void dp (D x, std::vector<R>& derivative) const
157  {
158  derivative.resize(1);
159  derivative[0] = 0.0;
160  }
161 
162  // derivative of ith Lagrange polynomial of degree k in one dimension
163  R dp (int i, D x) const
164  {
165  return 0.0;
166  }
167 
168  // value and derivative of all polynomials
169  void pdp (D x, std::vector<R>& value, std::vector<R>& derivative) const
170  {
171  value.resize(1);
172  derivative.resize(1);
173  value[0] = 1.0;
174  derivative[0] = 0.0;
175  }
176  };
177 
178  template<class D, class R>
180  {
181  public:
183  void p (D x, std::vector<R>& value) const
184  {
185  value.resize(2);
186  value[0] = 1.0;
187  value[1] = 2*x-1;
188  }
189 
190  // ith Lagrange polynomial of degree k in one dimension
191  R p (int i, D x) const
192  {
193  return (1-i) + i*(2*x-1);
194  }
195 
196  // derivative of all polynomials
197  void dp (D x, std::vector<R>& derivative) const
198  {
199  derivative.resize(2);
200  derivative[0] = 0.0;
201  derivative[1] = 2.0;
202  }
203 
204  // derivative of ith Lagrange polynomial of degree k in one dimension
205  R dp (int i, D x) const
206  {
207  return (1-i)*0 + i*(2);
208  }
209 
210  // value and derivative of all polynomials
211  void pdp (D x, std::vector<R>& value, std::vector<R>& derivative) const
212  {
213  value.resize(2);
214  derivative.resize(2);
215  value[0] = 1.0;
216  value[1] = 2*x-1;
217  derivative[0] = 0.0;
218  derivative[1] = 2.0;
219  }
220  };
221 
234  template<class D, class R, int k, int d>
236  {
237  enum { n = LegendreSize<k,d>::value };
239  mutable std::vector<std::vector<R> > v;
240  mutable std::vector<std::vector<R> > a;
241 
242  public:
243  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
244 
245  DGLegendreLocalBasis () : v(d,std::vector<R>(k+1,0.0)), a(d,std::vector<R>(k+1,0.0))
246  {}
247 
249  unsigned int size () const
250  {
251  return n;
252  }
253 
255  inline void evaluateFunction (const typename Traits::DomainType& x,
256  std::vector<typename Traits::RangeType>& value) const
257  {
258  // resize output vector
259  value.resize(n);
260 
261  // compute values of 1d basis functions in each direction
262  for (size_t j=0; j<d; j++) poly.p(x[j],v[j]);
263 
264  // now compute the product
265  for (size_t i=0; i<n; i++)
266  {
267  // convert index i to multiindex
268  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
269 
270  // initialize product
271  value[i] = 1.0;
272 
273  // dimension by dimension
274  for (int j=0; j<d; j++) value[i] *= v[j][alpha[j]];
275  }
276  }
277 
279  inline void
280  evaluateJacobian (const typename Traits::DomainType& x, // position
281  std::vector<typename Traits::JacobianType>& value) const // return value
282  {
283  // resize output vector
284  value.resize(size());
285 
286  // compute values of 1d basis functions in each direction
287  for (size_t j=0; j<d; j++) poly.pdp(x[j],v[j],a[j]);
288 
289  // Loop over all shape functions
290  for (size_t i=0; i<n; i++)
291  {
292  // convert index i to multiindex
293  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
294 
295  // Loop over all coordinate directions
296  for (int j=0; j<d; j++)
297  {
298  // Initialize: the overall expression is a product
299  value[i][0][j] = a[j][alpha[j]];
300 
301  // rest of the product
302  for (int l=0; l<d; l++)
303  if (l!=j)
304  value[i][0][j] *= v[l][alpha[l]];
305  }
306  }
307  }
308 
310  unsigned int order () const
311  {
312  return k;
313  }
314  };
315 
316 
318  template<class D, class R, int k, int d>
320  {
321  enum { n = LegendreSize<k,d>::value };
323  const LB lb;
324  Dune::GeometryType gt;
325 
326  public:
327 
328  DGLegendreLocalInterpolation () : gt(Dune::GeometryType::cube,d)
329  {}
330 
332  template<typename F, typename C>
333  void interpolate (const F& f, std::vector<C>& out) const
334  {
335  // select quadrature rule
336  typedef typename LB::Traits::RangeType RangeType;
337  const Dune::QuadratureRule<R,d>&
338  rule = Dune::QuadratureRules<R,d>::rule(gt,2*k);
339 
340  // prepare result
341  out.resize(n);
342  std::vector<R> diagonal(n);
343  for (int i=0; i<n; i++) { out[i] = 0.0; diagonal[i] = 0.0; }
344 
345  // loop over quadrature points
346  for (typename Dune::QuadratureRule<R,d>::const_iterator
347  it=rule.begin(); it!=rule.end(); ++it)
348  {
349  // evaluate function at quadrature point
350  typename LB::Traits::DomainType x;
351  RangeType y;
352  for (int i=0; i<d; i++) x[i] = it->position()[i];
353  f.evaluate(x,y);
354 
355  // evaluate the basis
356  std::vector<RangeType> phi(n);
357  lb.evaluateFunction(it->position(),phi);
358 
359  // do integration
360  for (int i=0; i<n; i++) {
361  out[i] += y*phi[i]*it->weight();
362  diagonal[i] += phi[i]*phi[i]*it->weight();
363  }
364  }
365  for (int i=0; i<n; i++) out[i] /= diagonal[i];
366  }
367  };
368 
375  template<int k, int d>
377  {
378  enum { n = LegendreSize<k,d>::value };
379 
380  public:
383  {
384  for (std::size_t i=0; i<n; i++)
385  li[i] = LocalKey(0,0,i);
386  }
387 
389  std::size_t size () const
390  {
391  return n;
392  }
393 
395  const LocalKey& localKey (std::size_t i) const
396  {
397  return li[i];
398  }
399 
400  private:
401  std::vector<LocalKey> li;
402  };
403 
404  } // end of LegendreStuff namespace
405 
408  template<class D, class R, int k, int d>
410  {
414 
415  public:
416  // static number of basis functions
418 
421  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
422 
426  {
427  gt.makeCube(d);
428  }
429 
432  const typename Traits::LocalBasisType& localBasis () const
433  {
434  return basis;
435  }
436 
439  const typename Traits::LocalCoefficientsType& localCoefficients () const
440  {
441  return coefficients;
442  }
443 
446  const typename Traits::LocalInterpolationType& localInterpolation () const
447  {
448  return interpolation;
449  }
450 
453  GeometryType type () const
454  {
455  return gt;
456  }
457 
459  {
460  return new DGLegendreLocalFiniteElement(*this);
461  }
462 
463  private:
464  LocalBasis basis;
465  LocalCoefficients coefficients;
466  LocalInterpolation interpolation;
467  GeometryType gt;
468  };
469 
471 
476  template<class Geometry, class RF, int k>
478  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
479  DGLegendreLocalFiniteElement<
480  typename Geometry::ctype, RF, k, Geometry::mydimension
481  >,
482  Geometry
483  >
484  {
486  typename Geometry::ctype, RF, k, Geometry::mydimension
487  > LFE;
488  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
489 
490  static const LFE lfe;
491 
492  public:
495  };
496 
497  template<class Geometry, class RF, int k>
500 
501 } // End Dune namespace
502 
503 #endif
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: dglegendre.hh:421
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: dglegendre.hh:333
R dp(int i, D x)
Definition: qkdg.hh:72
const Traits::LocalInterpolationType & localInterpolation() const
Definition: dglegendre.hh:446
STL namespace.
The 1d Legendre Polynomials (k=0,1 are specialized below)
Definition: dglegendre.hh:76
unsigned int size() const
number of shape functions
Definition: dglegendre.hh:249
DGLegendreLocalBasis()
Definition: dglegendre.hh:245
R p(int i, D x) const
Definition: dglegendre.hh:191
DGLegendreLocalCoefficients()
Standard constructor.
Definition: dglegendre.hh:382
determine degrees of freedom
Definition: dglegendre.hh:319
Definition: dglegendre.hh:31
std::size_t size() const
number of coefficients
Definition: dglegendre.hh:389
DGLegendreLocalFiniteElement()
Definition: dglegendre.hh:425
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &value) const
Evaluate all shape functions.
Definition: dglegendre.hh:255
Definition: adaptivity.hh:27
Lagrange shape functions of order k on the reference cube.
Definition: dglegendre.hh:235
R dp(int i, D x) const
Definition: dglegendre.hh:130
void p(D x, std::vector< R > &value) const
evaluate all polynomials at point x
Definition: dglegendre.hh:80
R dp(int i, D x) const
Definition: dglegendre.hh:205
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: dglegendre.hh:243
DGLegendreFiniteElementFactory()
default constructor
Definition: dglegendre.hh:494
void dp(D x, std::vector< R > &derivative) const
Definition: dglegendre.hh:98
Dune::FieldVector< int, d > multiindex(int i)
Definition: dglegendre.hh:63
unsigned int order() const
Polynomial order of the shape functions.
Definition: dglegendre.hh:310
R p(int i, D x) const
Definition: dglegendre.hh:90
void pdp(D x, std::vector< R > &value, std::vector< R > &derivative) const
Definition: dglegendre.hh:169
R dp(int i, D x) const
Definition: dglegendre.hh:163
DGLegendreLocalFiniteElement * clone() const
Definition: dglegendre.hh:458
Factory for global-valued DGLegendre elements.
Definition: dglegendre.hh:477
void p(D x, std::vector< R > &value) const
evaluate all polynomials at point x
Definition: dglegendre.hh:143
const P & p
Definition: constraints.hh:147
void dp(D x, std::vector< R > &derivative) const
Definition: dglegendre.hh:156
Definition: dglegendre.hh:409
Definition: dglegendre.hh:34
void pdp(D x, std::vector< R > &value, std::vector< R > &derivative) const
Definition: dglegendre.hh:114
DGLegendreLocalInterpolation()
Definition: dglegendre.hh:328
R p(int i, D x) const
Definition: dglegendre.hh:150
void dp(D x, std::vector< R > &derivative) const
Definition: dglegendre.hh:197
GeometryType type() const
Definition: dglegendre.hh:453
void pdp(D x, std::vector< R > &value, std::vector< R > &derivative) const
Definition: dglegendre.hh:211
void p(D x, std::vector< R > &value) const
evaluate all polynomials at point x
Definition: dglegendre.hh:183
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &value) const
Evaluate Jacobian of all shape functions.
Definition: dglegendre.hh:280
Layout map for Q1 elements.
Definition: dglegendre.hh:376
const LocalKey & localKey(std::size_t i) const
get i&#39;th index
Definition: dglegendre.hh:395
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: dglegendre.hh:439
const Traits::LocalBasisType & localBasis() const
Definition: dglegendre.hh:432