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