dune-pdelab 2.7-git
Loading...
Searching...
No Matches
qkdglagrange.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#ifndef DUNE_PDELAB_FINITEELEMENT_QKDGLAGRANGE_HH
5#define DUNE_PDELAB_FINITEELEMENT_QKDGLAGRANGE_HH
6
7#include <numeric>
8
9#include <dune/localfunctions/common/localbasis.hh>
10#include <dune/localfunctions/common/localkey.hh>
11#include <dune/localfunctions/common/localfiniteelementtraits.hh>
12#include <dune/localfunctions/common/localtoglobaladaptors.hh>
13
14namespace Dune
15{
16 namespace QkStuff
17 {
18 // This is the main class
19 // usage: QkSize<2,3>::value
20 // k is the polynomial degree,
21 // n is the space dimension
22 template<int k, int n>
23 struct QkSize
24 {
25 enum{
27 };
28 };
29
30 template<>
31 struct QkSize<0,1>
32 {
33 enum{
34 value=1
35 };
36 };
37
38 template<int k>
39 struct QkSize<k,1>
40 {
41 enum{
42 value=k+1
43 };
44 };
45
46 template<int n>
47 struct QkSize<0,n>
48 {
49 enum{
50 value=1
51 };
52 };
53
54 // ith Lagrange polynomial of degree k in one dimension
55 template<class D, class R, int k>
56 R p (int i, D x)
57 {
58 R result(1.0);
59 for (int j=0; j<=k; j++)
60 if (j!=i) result *= (k*x-j)/(i-j);
61 return result;
62 }
63
64 // derivative of ith Lagrange polynomial of degree k in one dimension
65 template<class D, class R, int k>
66 R dp (int i, D x)
67 {
68 R result(0.0);
69
70 for (int j=0; j<=k; j++)
71 if (j!=i)
72 {
73 R prod( (k*1.0)/(i-j) );
74 for (int l=0; l<=k; l++)
75 if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
76 result += prod;
77 }
78 return result;
79 }
80
82 template<class D, class R, int k>
84 {
85 public:
86 // ith Lagrange polynomial of degree k in one dimension
87 R p (int i, D x) const
88 {
89 R result(1.0);
90 for (int j=0; j<=k; j++)
91 if (j!=i) result *= (k*x-j)/(i-j);
92 return result;
93 }
94
95 // derivative of ith Lagrange polynomial of degree k in one dimension
96 R dp (int i, D x) const
97 {
98 R result(0.0);
99
100 for (int j=0; j<=k; j++)
101 if (j!=i)
102 {
103 R prod( (k*1.0)/(i-j) );
104 for (int l=0; l<=k; l++)
105 if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
106 result += prod;
107 }
108 return result;
109 }
110
111 // get ith Lagrange point
112 R x (int i)
113 {
114 return i/((1.0)*(k+1));
115 }
116 };
117
118 template<int k, int d>
119 Dune::FieldVector<int,d> multiindex (int i)
120 {
121 Dune::FieldVector<int,d> alpha;
122 for (int j=0; j<d; j++)
123 {
124 alpha[j] = i % (k+1);
125 i = i/(k+1);
126 }
127 return alpha;
128 }
129
142 template<class D, class R, int k, int d>
144 {
145 enum{ n = QkSize<k,d>::value };
146 public:
147 typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
148
150 unsigned int size () const
151 {
152 return QkSize<k,d>::value;
153 }
154
156 inline void evaluateFunction (const typename Traits::DomainType& in,
157 std::vector<typename Traits::RangeType>& out) const
158 {
159 out.resize(size());
160 for (size_t i=0; i<size(); i++)
161 {
162 // convert index i to multiindex
163 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
164
165 // initialize product
166 out[i] = 1.0;
167
168 // dimension by dimension
169 for (int j=0; j<d; j++)
170 out[i] *= p<D,R,k>(alpha[j],in[j]);
171 }
172 }
173
175 inline void
176 evaluateJacobian (const typename Traits::DomainType& in, // position
177 std::vector<typename Traits::JacobianType>& out) const // return value
178 {
179 out.resize(size());
180
181 // Loop over all shape functions
182 for (size_t i=0; i<size(); i++)
183 {
184 // convert index i to multiindex
185 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
186
187 // Loop over all coordinate directions
188 for (int j=0; j<d; j++)
189 {
190 // Initialize: the overall expression is a product
191 // if j-th bit of i is set to -1, else 1
192 out[i][0][j] = dp<D,R,k>(alpha[j],in[j]);
193
194 // rest of the product
195 for (int l=0; l<d; l++)
196 if (l!=j)
197 out[i][0][j] *= p<D,R,k>(alpha[l],in[l]);
198 }
199 }
200 }
201
203 void partial(const std::array<unsigned int,Traits::dimDomain>& DUNE_UNUSED(order),
204 const typename Traits::DomainType& DUNE_UNUSED(in),
205 std::vector<typename Traits::RangeType>& DUNE_UNUSED(out)) const
206 {
207 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
208 if (totalOrder == 0) {
209 evaluateFunction(in, out);
210 } else {
211 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
212 }
213 }
214
216 unsigned int order () const
217 {
218 return k;
219 }
220 };
221
228 template<int k, int d>
230 {
231 public:
234 {
235 for (std::size_t i=0; i<QkSize<k,d>::value; i++)
236 li[i] = LocalKey(0,0,i);
237 }
238
240 std::size_t size () const
241 {
242 return QkSize<k,d>::value;
243 }
244
246 const LocalKey& localKey (std::size_t i) const
247 {
248 return li[i];
249 }
250
251 private:
252 std::vector<LocalKey> li;
253 };
254
256 template<int k, int d, class LB>
258 {
259 public:
260
262 template<typename F, typename C>
263 void interpolate (const F& f, std::vector<C>& out) const
264 {
265 typename LB::Traits::DomainType x;
266
267 out.resize(QkSize<k,d>::value);
268
269 for (int i=0; i<QkSize<k,d>::value; i++)
270 {
271 // convert index i to multiindex
272 Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
273
274 // Generate coordinate of the i-th Lagrange point
275 for (int j=0; j<d; j++)
276 x[j] = (1.0*alpha[j])/k;
277
278 out[i] = f(x);
279 }
280 }
281 };
282
284 template<int d, class LB>
286 {
287 public:
289 template<typename F, typename C>
290 void interpolate (const F& f, std::vector<C>& out) const
291 {
292 typename LB::Traits::DomainType x(0.5);
293
294 out.resize(1);
295 out[0] = f(x);
296 }
297 };
298
299 }
300
303 template<class D, class R, int k, int d>
305 {
309
310 public:
311 // static number of basis functions
313
316 typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
317
321 : gt(GeometryTypes::cube(d))
322 {}
323
326 const typename Traits::LocalBasisType& localBasis () const
327 {
328 return basis;
329 }
330
333 const typename Traits::LocalCoefficientsType& localCoefficients () const
334 {
335 return coefficients;
336 }
337
340 const typename Traits::LocalInterpolationType& localInterpolation () const
341 {
342 return interpolation;
343 }
344
347 std::size_t size() const
348 {
349 return basis.size();
350 }
351
354 GeometryType type () const
355 {
356 return gt;
357 }
358
360 {
361 return new QkDGLagrangeLocalFiniteElement(*this);
362 }
363
364 private:
365 LocalBasis basis;
366 LocalCoefficients coefficients;
367 LocalInterpolation interpolation;
368 GeometryType gt;
369 };
370
372
377 template<class Geometry, class RF, int k>
379 public ScalarLocalToGlobalFiniteElementAdaptorFactory<
380 QkDGLagrangeLocalFiniteElement<
381 typename Geometry::ctype, RF, k, Geometry::mydimension
382 >,
383 Geometry
384 >
385 {
387 typename Geometry::ctype, RF, k, Geometry::mydimension
388 > LFE;
389 typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
390
391 static const LFE lfe;
392
393 public:
396 };
397
398 template<class Geometry, class RF, int k>
399 const typename QkDGFiniteElementFactory<Geometry, RF, k>::LFE
400 QkDGFiniteElementFactory<Geometry, RF, k>::lfe;
401
402}
403
404
405#endif // DUNE_PDELAB_FINITEELEMENT_QKDGLAGRANGE_HH
const P & p
Definition constraints.hh:148
For backward compatibility – Do not use this!
Definition adaptivity.hh:28
R dp(int i, D x)
Definition qkdglagrange.hh:66
Dune::FieldVector< int, d > multiindex(int i)
Definition qkdglagrange.hh:119
Definition qkdglagrange.hh:24
@ value
Definition qkdglagrange.hh:26
Lagrange polynomials of degree k at equidistant points as a class.
Definition qkdglagrange.hh:84
R x(int i)
Definition qkdglagrange.hh:112
R dp(int i, D x) const
Definition qkdglagrange.hh:96
R p(int i, D x) const
Definition qkdglagrange.hh:87
Lagrange shape functions of order k on the reference cube.
Definition qkdglagrange.hh:144
unsigned int size() const
number of shape functions
Definition qkdglagrange.hh:150
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition qkdglagrange.hh:176
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 qkdglagrange.hh:203
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition qkdglagrange.hh:147
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition qkdglagrange.hh:156
unsigned int order() const
Polynomial order of the shape functions.
Definition qkdglagrange.hh:216
Layout map for Q1 elements.
Definition qkdglagrange.hh:230
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition qkdglagrange.hh:246
QkDGLocalCoefficients()
Standard constructor.
Definition qkdglagrange.hh:233
std::size_t size() const
number of coefficients
Definition qkdglagrange.hh:240
Definition qkdglagrange.hh:258
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition qkdglagrange.hh:263
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition qkdglagrange.hh:290
Definition qkdglagrange.hh:305
std::size_t size() const
Definition qkdglagrange.hh:347
const Traits::LocalInterpolationType & localInterpolation() const
Definition qkdglagrange.hh:340
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition qkdglagrange.hh:316
const Traits::LocalCoefficientsType & localCoefficients() const
Definition qkdglagrange.hh:333
GeometryType type() const
Definition qkdglagrange.hh:354
QkDGLagrangeLocalFiniteElement()
Definition qkdglagrange.hh:320
QkDGLagrangeLocalFiniteElement * clone() const
Definition qkdglagrange.hh:359
@ n
Definition qkdglagrange.hh:312
const Traits::LocalBasisType & localBasis() const
Definition qkdglagrange.hh:326
Factory for global-valued QkDG elements.
Definition qkdglagrange.hh:385
QkDGFiniteElementFactory()
default constructor
Definition qkdglagrange.hh:395
static const unsigned int value
Definition gridfunctionspace/tags.hh:139