dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
23namespace 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>
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
const P & p
Definition constraints.hh:148
STL namespace.
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
DGLegendreLocalCoefficients()
Standard constructor.
Definition qkdglegendre.hh:395
std::size_t size() const
number of coefficients
Definition qkdglegendre.hh:402
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition qkdglegendre.hh:408
Definition qkdglegendre.hh:423
const Traits::LocalCoefficientsType & localCoefficients() const
Definition qkdglegendre.hh:451
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition qkdglegendre.hh:434
const Traits::LocalInterpolationType & localInterpolation() const
Definition qkdglegendre.hh:458
QkDGLegendreLocalFiniteElement * clone() const
Definition qkdglegendre.hh:477
GeometryType type() const
Definition qkdglegendre.hh:472
std::size_t size() const
Definition qkdglegendre.hh:465
const Traits::LocalBasisType & localBasis() const
Definition qkdglegendre.hh:444
@ 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