DUNE PDELab (git)

qkdglegendre.hh
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
14
15#include <dune/geometry/type.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
25namespace 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{
37 value=(k+1)*LegendreSize<k,n-1>::value
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 {
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>
142 class LegendrePolynomials1d<D,R,0>
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>
182 class LegendrePolynomials1d<D,R,1>
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:
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;
340
341 public:
342
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);
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
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
434 enum { n = LegendreStuff::LegendreSize<k,d>::value };
435
439
443 : gt(GeometryTypes::cube(d))
444 {}
445
448 const typename Traits::LocalBasisType& localBasis () const
449 {
450 return basis;
451 }
452
456 {
457 return coefficients;
458 }
459
463 {
464 return interpolation;
465 }
466
469 std::size_t size() const
470 {
471 return basis.size();
472 }
473
477 {
478 return gt;
479 }
480
481 QkDGLegendreLocalFiniteElement* clone () const
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>
502 QkDGLegendreLocalFiniteElement<
503 typename Geometry::ctype, RF, k, Geometry::mydimension
504 >,
505 Geometry
506 >
507 {
510 > LFE;
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
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
GridImp::ctype ctype
define type used for coordinates in grid module
Definition: geometry.hh:100
static constexpr int mydimension
geometry dimension
Definition: geometry.hh:94
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
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglegendre.hh:325
Layout map for Q1 elements.
Definition: qkdglegendre.hh:394
DGLegendreLocalCoefficients()
Standard constructor.
Definition: qkdglegendre.hh:399
std::size_t size() const
number of coefficients
Definition: qkdglegendre.hh:406
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: qkdglegendre.hh:412
determine degrees of freedom
Definition: qkdglegendre.hh:335
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglegendre.hh:348
The 1d Legendre Polynomials (k=0,1 are specialized below)
Definition: qkdglegendre.hh:80
void p(D x, std::vector< R > &value) const
evaluate all polynomials at point x
Definition: qkdglegendre.hh:83
Describe position of one degree of freedom.
Definition: localkey.hh:24
Default exception for dummy implementations.
Definition: exceptions.hh:263
Factory for global-valued DGLegendre elements.
Definition: qkdglegendre.hh:507
QkDGLegendreFiniteElementFactory()
default constructor
Definition: qkdglegendre.hh:517
Definition: qkdglegendre.hh:427
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglegendre.hh:455
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglegendre.hh:438
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglegendre.hh:462
GeometryType type() const
Definition: qkdglegendre.hh:476
std::size_t size() const
Definition: qkdglegendre.hh:469
const Traits::LocalBasisType & localBasis() const
Definition: qkdglegendre.hh:448
QkDGLegendreLocalFiniteElement()
Definition: qkdglegendre.hh:442
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
Factory for ScalarLocalToGlobalFiniteElementAdaptor objects.
Definition: localtoglobaladaptors.hh:244
Definition of the DUNE_NO_DEPRECATED_* macros.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
TrigonometricFunction< K, -cosFactor, sinFactor > derivative(const TrigonometricFunction< K, sinFactor, cosFactor > &f)
Obtain derivative of TrigonometricFunction function.
Definition: trigonometricfunction.hh:43
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
Dune namespace.
Definition: alignedallocator.hh:13
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
R RangeType
range type
Definition: localbasis.hh:52
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)