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/localtoglobaladaptors.hh>
13 
14 namespace 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>
285  class QkLocalInterpolation<0,d,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:
395  QkDGFiniteElementFactory() : Base(lfe) {}
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
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdglagrange.hh:119
R p(int i, D x)
Definition: qkdglagrange.hh:56
R dp(int i, D x)
Definition: qkdglagrange.hh:66
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
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglagrange.hh:333
const Traits::LocalBasisType & localBasis() const
Definition: qkdglagrange.hh:326
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglagrange.hh:340
std::size_t size() const
Definition: qkdglagrange.hh:347
QkDGLagrangeLocalFiniteElement * clone() const
Definition: qkdglagrange.hh:359
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglagrange.hh:316
GeometryType type() const
Definition: qkdglagrange.hh:354
QkDGLagrangeLocalFiniteElement()
Definition: qkdglagrange.hh:320
@ n
Definition: qkdglagrange.hh:312
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