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