dune-pdelab  2.7-git
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 
16 namespace 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>
288  class QkLocalInterpolation<0,d,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:
399  QkDGFiniteElementFactory() : Base(lfe) {}
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
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:121
R p(int i, D x)
Definition: qkdglagrange.hh:58
R dp(int i, D x)
Definition: qkdglagrange.hh:68
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
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglagrange.hh:337
const Traits::LocalBasisType & localBasis() const
Definition: qkdglagrange.hh:330
@ n
Definition: qkdglagrange.hh:316
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglagrange.hh:344
std::size_t size() const
Definition: qkdglagrange.hh:351
QkDGLagrangeLocalFiniteElement * clone() const
Definition: qkdglagrange.hh:363
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglagrange.hh:320
GeometryType type() const
Definition: qkdglagrange.hh:358
QkDGLagrangeLocalFiniteElement()
Definition: qkdglagrange.hh:324
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