dune-pdelab  2.5-dev
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 <dune/localfunctions/common/localbasis.hh>
8 #include <dune/localfunctions/common/localkey.hh>
9 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
10 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
11 
12 namespace Dune
13 {
14  namespace QkStuff
15  {
16  // This is the main class
17  // usage: QkSize<2,3>::value
18  // k is the polynomial degree,
19  // n is the space dimension
20  template<int k, int n>
21  struct QkSize
22  {
23  enum{
25  };
26  };
27 
28  template<>
29  struct QkSize<0,1>
30  {
31  enum{
33  };
34  };
35 
36  template<int k>
37  struct QkSize<k,1>
38  {
39  enum{
40  value=k+1
41  };
42  };
43 
44  template<int n>
45  struct QkSize<0,n>
46  {
47  enum{
49  };
50  };
51 
52  // ith Lagrange polynomial of degree k in one dimension
53  template<class D, class R, int k>
54  R p (int i, D x)
55  {
56  R result(1.0);
57  for (int j=0; j<=k; j++)
58  if (j!=i) result *= (k*x-j)/(i-j);
59  return result;
60  }
61 
62  // derivative of ith Lagrange polynomial of degree k in one dimension
63  template<class D, class R, int k>
64  R dp (int i, D x)
65  {
66  R result(0.0);
67 
68  for (int j=0; j<=k; j++)
69  if (j!=i)
70  {
71  R prod( (k*1.0)/(i-j) );
72  for (int l=0; l<=k; l++)
73  if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
74  result += prod;
75  }
76  return result;
77  }
78 
80  template<class D, class R, int k>
82  {
83  public:
84  // ith Lagrange polynomial of degree k in one dimension
85  R p (int i, D x) const
86  {
87  R result(1.0);
88  for (int j=0; j<=k; j++)
89  if (j!=i) result *= (k*x-j)/(i-j);
90  return result;
91  }
92 
93  // derivative of ith Lagrange polynomial of degree k in one dimension
94  R dp (int i, D x) const
95  {
96  R result(0.0);
97 
98  for (int j=0; j<=k; j++)
99  if (j!=i)
100  {
101  R prod( (k*1.0)/(i-j) );
102  for (int l=0; l<=k; l++)
103  if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
104  result += prod;
105  }
106  return result;
107  }
108 
109  // get ith Lagrange point
110  R x (int i)
111  {
112  return i/((1.0)*(k+1));
113  }
114  };
115 
116  template<int k, int d>
117  Dune::FieldVector<int,d> multiindex (int i)
118  {
119  Dune::FieldVector<int,d> alpha;
120  for (int j=0; j<d; j++)
121  {
122  alpha[j] = i % (k+1);
123  i = i/(k+1);
124  }
125  return alpha;
126  }
127 
140  template<class D, class R, int k, int d>
142  {
143  enum{ n = QkSize<k,d>::value };
144  public:
145  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
146 
148  unsigned int size () const
149  {
150  return QkSize<k,d>::value;
151  }
152 
154  inline void evaluateFunction (const typename Traits::DomainType& in,
155  std::vector<typename Traits::RangeType>& out) const
156  {
157  out.resize(size());
158  for (size_t i=0; i<size(); i++)
159  {
160  // convert index i to multiindex
161  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
162 
163  // initialize product
164  out[i] = 1.0;
165 
166  // dimension by dimension
167  for (int j=0; j<d; j++)
168  out[i] *= p<D,R,k>(alpha[j],in[j]);
169  }
170  }
171 
173  inline void
174  evaluateJacobian (const typename Traits::DomainType& in, // position
175  std::vector<typename Traits::JacobianType>& out) const // return value
176  {
177  out.resize(size());
178 
179  // Loop over all shape functions
180  for (size_t i=0; i<size(); i++)
181  {
182  // convert index i to multiindex
183  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
184 
185  // Loop over all coordinate directions
186  for (int j=0; j<d; j++)
187  {
188  // Initialize: the overall expression is a product
189  // if j-th bit of i is set to -1, else 1
190  out[i][0][j] = dp<D,R,k>(alpha[j],in[j]);
191 
192  // rest of the product
193  for (int l=0; l<d; l++)
194  if (l!=j)
195  out[i][0][j] *= p<D,R,k>(alpha[l],in[l]);
196  }
197  }
198  }
199 
201  unsigned int order () const
202  {
203  return k;
204  }
205  };
206 
213  template<int k, int d>
215  {
216  public:
219  {
220  for (std::size_t i=0; i<QkSize<k,d>::value; i++)
221  li[i] = LocalKey(0,0,i);
222  }
223 
225  std::size_t size () const
226  {
227  return QkSize<k,d>::value;
228  }
229 
231  const LocalKey& localKey (std::size_t i) const
232  {
233  return li[i];
234  }
235 
236  private:
237  std::vector<LocalKey> li;
238  };
239 
241  template<int k, int d, class LB>
243  {
244  public:
245 
247  template<typename F, typename C>
248  void interpolate (const F& f, std::vector<C>& out) const
249  {
250  typename LB::Traits::DomainType x;
251  typename LB::Traits::RangeType y;
252 
253  out.resize(QkSize<k,d>::value);
254 
255  for (int i=0; i<QkSize<k,d>::value; i++)
256  {
257  // convert index i to multiindex
258  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
259 
260  // Generate coordinate of the i-th Lagrange point
261  for (int j=0; j<d; j++)
262  x[j] = (1.0*alpha[j])/k;
263 
264  f.evaluate(x,y); out[i] = y;
265  }
266  }
267  };
268 
270  template<int d, class LB>
271  class QkLocalInterpolation<0,d,LB>
272  {
273  public:
275  template<typename F, typename C>
276  void interpolate (const F& f, std::vector<C>& out) const
277  {
278  typename LB::Traits::DomainType x(0.5);
279  typename LB::Traits::RangeType y;
280  f.evaluate(x,y);
281  out.resize(1);
282  out[0] = y;
283  }
284  };
285 
286  }
287 
290  template<class D, class R, int k, int d>
292  {
296 
297  public:
298  // static number of basis functions
300 
303  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
304 
308  : gt(GeometryTypes::cube(d))
309  {}
310 
313  const typename Traits::LocalBasisType& localBasis () const
314  {
315  return basis;
316  }
317 
320  const typename Traits::LocalCoefficientsType& localCoefficients () const
321  {
322  return coefficients;
323  }
324 
327  const typename Traits::LocalInterpolationType& localInterpolation () const
328  {
329  return interpolation;
330  }
331 
334  GeometryType type () const
335  {
336  return gt;
337  }
338 
340  {
341  return new QkDGLagrangeLocalFiniteElement(*this);
342  }
343 
344  private:
345  LocalBasis basis;
346  LocalCoefficients coefficients;
347  LocalInterpolation interpolation;
348  GeometryType gt;
349  };
350 
352 
357  template<class Geometry, class RF, int k>
359  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
360  QkDGLagrangeLocalFiniteElement<
361  typename Geometry::ctype, RF, k, Geometry::mydimension
362  >,
363  Geometry
364  >
365  {
367  typename Geometry::ctype, RF, k, Geometry::mydimension
368  > LFE;
369  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
370 
371  static const LFE lfe;
372 
373  public:
375  QkDGFiniteElementFactory() : Base(lfe) {}
376  };
377 
378  template<class Geometry, class RF, int k>
381 
382 }
383 
384 
385 #endif // DUNE_PDELAB_FINITEELEMENT_QKDGLAGRANGE_HH
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglagrange.hh:201
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglagrange.hh:174
unsigned int size() const
number of shape functions
Definition: qkdglagrange.hh:148
R dp(int i, D x) const
Definition: qkdglagrange.hh:94
Lagrange shape functions of order k on the reference cube.
Definition: qkdglagrange.hh:141
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglagrange.hh:327
Factory for global-valued QkDG elements.
Definition: qkdglagrange.hh:358
const Traits::LocalBasisType & localBasis() const
Definition: qkdglagrange.hh:313
Definition: qkdglagrange.hh:242
Lagrange polynomials of degree k at equidistant points as a class.
Definition: qkdglagrange.hh:81
QkDGLagrangeLocalFiniteElement()
Definition: qkdglagrange.hh:307
QkDGLagrangeLocalFiniteElement * clone() const
Definition: qkdglagrange.hh:339
R x(int i)
Definition: qkdglagrange.hh:110
Definition: qkdglagrange.hh:291
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
R dp(int i, D x)
Definition: qkdglagrange.hh:64
Definition: qkdglagrange.hh:24
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglagrange.hh:145
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:117
R p(int i, D x)
Definition: qkdglagrange.hh:54
R p(int i, D x) const
Definition: qkdglagrange.hh:85
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglagrange.hh:248
QkDGFiniteElementFactory()
default constructor
Definition: qkdglagrange.hh:375
const LocalKey & localKey(std::size_t i) const
get i&#39;th index
Definition: qkdglagrange.hh:231
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglagrange.hh:154
QkDGLocalCoefficients()
Standard constructor.
Definition: qkdglagrange.hh:218
Layout map for Q1 elements.
Definition: qkdglagrange.hh:214
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglagrange.hh:320
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglagrange.hh:276
GeometryType type() const
Definition: qkdglagrange.hh:334
std::size_t size() const
number of coefficients
Definition: qkdglagrange.hh:225
Definition: qkdglagrange.hh:21
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglagrange.hh:303