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