dune-pdelab  2.4-dev
qkdg.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_QkDG_LOCALFINITEELEMENT_HH
5 #define DUNE_QkDG_LOCALFINITEELEMENT_HH
6 
7 #include <dune/common/fvector.hh>
8 #include <dune/common/deprecated.hh>
9 
10 #include <dune/geometry/type.hh>
11 
12 #include <dune/localfunctions/common/localbasis.hh>
13 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
14 #include <dune/localfunctions/common/localkey.hh>
15 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
16 
17 #include "finiteelementmap.hh"
18 
19 namespace Dune
20 {
21 
22  namespace QkStuff
23  {
24  // This is the main class
25  // usage: QkSize<2,3>::value
26  // k is the polynomial degree,
27  // n is the space dimension
28  template<int k, int n>
29  struct QkSize
30  {
31  enum{
33  };
34  };
35 
36  template<>
37  struct QkSize<0,1>
38  {
39  enum{
41  };
42  };
43 
44  template<int k>
45  struct QkSize<k,1>
46  {
47  enum{
48  value=k+1
49  };
50  };
51 
52  template<int n>
53  struct QkSize<0,n>
54  {
55  enum{
57  };
58  };
59 
60  // ith Lagrange polynomial of degree k in one dimension
61  template<class D, class R, int k>
62  R p (int i, D x)
63  {
64  R result(1.0);
65  for (int j=0; j<=k; j++)
66  if (j!=i) result *= (k*x-j)/(i-j);
67  return result;
68  }
69 
70  // derivative of ith Lagrange polynomial of degree k in one dimension
71  template<class D, class R, int k>
72  R dp (int i, D x)
73  {
74  R result(0.0);
75 
76  for (int j=0; j<=k; j++)
77  if (j!=i)
78  {
79  R prod( (k*1.0)/(i-j) );
80  for (int l=0; l<=k; l++)
81  if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
82  result += prod;
83  }
84  return result;
85  }
86 
88  template<class D, class R, int k>
90  {
91  public:
92  // ith Lagrange polynomial of degree k in one dimension
93  R p (int i, D x) const
94  {
95  R result(1.0);
96  for (int j=0; j<=k; j++)
97  if (j!=i) result *= (k*x-j)/(i-j);
98  return result;
99  }
100 
101  // derivative of ith Lagrange polynomial of degree k in one dimension
102  R dp (int i, D x) const
103  {
104  R result(0.0);
105 
106  for (int j=0; j<=k; j++)
107  if (j!=i)
108  {
109  R prod( (k*1.0)/(i-j) );
110  for (int l=0; l<=k; l++)
111  if (l!=i && l!=j) prod *= (k*x-l)/(i-l);
112  result += prod;
113  }
114  return result;
115  }
116 
117  // get ith Lagrange point
118  R x (int i)
119  {
120  return i/((1.0)*(k+1));
121  }
122  };
123 
124  template<int k, int d>
125  Dune::FieldVector<int,d> multiindex (int i)
126  {
127  Dune::FieldVector<int,d> alpha;
128  for (int j=0; j<d; j++)
129  {
130  alpha[j] = i % (k+1);
131  i = i/(k+1);
132  }
133  return alpha;
134  }
135 
148  template<class D, class R, int k, int d>
150  {
151  enum{ n = QkSize<k,d>::value };
152  public:
153  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
154 
156  unsigned int size () const
157  {
158  return QkSize<k,d>::value;
159  }
160 
162  inline void evaluateFunction (const typename Traits::DomainType& in,
163  std::vector<typename Traits::RangeType>& out) const
164  {
165  out.resize(size());
166  for (size_t i=0; i<size(); i++)
167  {
168  // convert index i to multiindex
169  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
170 
171  // initialize product
172  out[i] = 1.0;
173 
174  // dimension by dimension
175  for (int j=0; j<d; j++)
176  out[i] *= p<D,R,k>(alpha[j],in[j]);
177  }
178  }
179 
181  inline void
182  evaluateJacobian (const typename Traits::DomainType& in, // position
183  std::vector<typename Traits::JacobianType>& out) const // return value
184  {
185  out.resize(size());
186 
187  // Loop over all shape functions
188  for (size_t i=0; i<size(); i++)
189  {
190  // convert index i to multiindex
191  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
192 
193  // Loop over all coordinate directions
194  for (int j=0; j<d; j++)
195  {
196  // Initialize: the overall expression is a product
197  // if j-th bit of i is set to -1, else 1
198  out[i][0][j] = dp<D,R,k>(alpha[j],in[j]);
199 
200  // rest of the product
201  for (int l=0; l<d; l++)
202  if (l!=j)
203  out[i][0][j] *= p<D,R,k>(alpha[l],in[l]);
204  }
205  }
206  }
207 
209  unsigned int order () const
210  {
211  return k;
212  }
213  };
214 
221  template<int k, int d>
223  {
224  public:
227  {
228  for (std::size_t i=0; i<QkSize<k,d>::value; i++)
229  li[i] = LocalKey(0,0,i);
230  }
231 
233  std::size_t size () const
234  {
235  return QkSize<k,d>::value;
236  }
237 
239  const LocalKey& localKey (std::size_t i) const
240  {
241  return li[i];
242  }
243 
244  private:
245  std::vector<LocalKey> li;
246  };
247 
249  template<int k, int d, class LB>
251  {
252  public:
253 
255  template<typename F, typename C>
256  void interpolate (const F& f, std::vector<C>& out) const
257  {
258  typename LB::Traits::DomainType x;
259  typename LB::Traits::RangeType y;
260 
261  out.resize(QkSize<k,d>::value);
262 
263  for (int i=0; i<QkSize<k,d>::value; i++)
264  {
265  // convert index i to multiindex
266  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
267 
268  // Generate coordinate of the i-th Lagrange point
269  for (int j=0; j<d; j++)
270  x[j] = (1.0*alpha[j])/k;
271 
272  f.evaluate(x,y); out[i] = y;
273  }
274  }
275  };
276 
278  template<int d, class LB>
279  class QkLocalInterpolation<0,d,LB>
280  {
281  public:
283  template<typename F, typename C>
284  void interpolate (const F& f, std::vector<C>& out) const
285  {
286  typename LB::Traits::DomainType x(0.5);
287  typename LB::Traits::RangeType y;
288  f.evaluate(x,y);
289  out.resize(1);
290  out[0] = y;
291  }
292  };
293 
294  }
295 
298  template<class D, class R, int k, int d>
300  {
304 
305  public:
306  // static number of basis functions
308 
311  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
312 
316  {
317  gt.makeCube(d);
318  }
319 
322  const typename Traits::LocalBasisType& localBasis () const
323  {
324  return basis;
325  }
326 
329  const typename Traits::LocalCoefficientsType& localCoefficients () const
330  {
331  return coefficients;
332  }
333 
336  const typename Traits::LocalInterpolationType& localInterpolation () const
337  {
338  return interpolation;
339  }
340 
343  GeometryType type () const
344  {
345  return gt;
346  }
347 
349  {
350  return new QkDGLocalFiniteElement(*this);
351  }
352 
353  private:
354  LocalBasis basis;
355  LocalCoefficients coefficients;
356  LocalInterpolation interpolation;
357  GeometryType gt;
358  };
359 
361 
366  template<class Geometry, class RF, int k>
368  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
369  QkDGLocalFiniteElement<
370  typename Geometry::ctype, RF, k, Geometry::mydimension
371  >,
372  Geometry
373  >
374  {
375  typedef QkDGLocalFiniteElement<
376  typename Geometry::ctype, RF, k, Geometry::mydimension
377  > LFE;
378  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
379 
380  static const LFE lfe;
381 
382  public:
384  QkDGFiniteElementFactory() : Base(lfe) {}
385  };
386 
387  template<class Geometry, class RF, int k>
388  const typename QkDGFiniteElementFactory<Geometry, RF, k>::LFE
389  QkDGFiniteElementFactory<Geometry, RF, k>::lfe;
390 
391 
395  template<class D, class R, int k, int d>
397 
398 
402  template<class D, class R, int d>
403  class QkCGLocalFiniteElement<D,R,1,d>
404  {
405  enum{ k = 1 };
407 
408  class QkCGLocalCoefficients
409  {
410  enum{ k = 1 };
411  public:
413  QkCGLocalCoefficients () : li(QkStuff::QkSize<k,d>::value)
414  {
415  for (std::size_t i=0; i<QkStuff::QkSize<k,d>::value; i++)
416  li[i] = LocalKey(i,d,0);
417  }
418 
420  std::size_t size () const
421  {
423  }
424 
426  const LocalKey& localKey (std::size_t i) const
427  {
428  return li[i];
429  }
430 
431  private:
432  std::vector<LocalKey> li;
433  };
434 
436 
437  public:
438  // static number of basis functions
440 
443  typedef LocalFiniteElementTraits<LocalBasis,QkCGLocalCoefficients,LocalInterpolation> Traits;
444 
448  {
449  gt.makeCube(d);
450  }
451 
454  const typename Traits::LocalBasisType& localBasis () const
455  {
456  return basis;
457  }
458 
461  const typename Traits::LocalCoefficientsType& localCoefficients () const
462  {
463  return coefficients;
464  }
465 
468  const typename Traits::LocalInterpolationType& localInterpolation () const
469  {
470  return interpolation;
471  }
472 
475  GeometryType type () const
476  {
477  return gt;
478  }
479 
481  {
482  return new QkCGLocalFiniteElement(*this);
483  }
484 
485  private:
486  LocalBasis basis;
487  QkCGLocalCoefficients coefficients;
488  LocalInterpolation interpolation;
489  GeometryType gt;
490  };
491 
495  template<class D, class R>
496  class QkCGLocalFiniteElement<D,R,2,2>
497  {
498  enum{ k = 2 };
499  enum{ d = 2 };
501 
502  class QkCGLocalCoefficients
503  {
504  enum{ k = 2 };
505  enum{ d = 2 };
506  public:
508  QkCGLocalCoefficients () : li(QkStuff::QkSize<k,d>::value)
509  {
510  li[0] = LocalKey(0,2,0);
511  li[1] = LocalKey(2,1,0);
512  li[2] = LocalKey(1,2,0);
513  li[3] = LocalKey(0,1,0);
514  li[4] = LocalKey(0,0,0);
515  li[5] = LocalKey(1,1,0);
516  li[6] = LocalKey(2,2,0);
517  li[7] = LocalKey(3,1,0);
518  li[8] = LocalKey(3,2,0);
519  }
520 
522  std::size_t size () const
523  {
525  }
526 
528  const LocalKey& localKey (std::size_t i) const
529  {
530  return li[i];
531  }
532 
533  private:
534  std::vector<LocalKey> li;
535  };
536 
538 
539  public:
540  // static number of basis functions
542 
545  typedef LocalFiniteElementTraits<LocalBasis,QkCGLocalCoefficients,LocalInterpolation> Traits;
546 
550  {
551  gt.makeCube(d);
552  }
553 
556  const typename Traits::LocalBasisType& localBasis () const
557  {
558  return basis;
559  }
560 
563  const typename Traits::LocalCoefficientsType& localCoefficients () const
564  {
565  return coefficients;
566  }
567 
570  const typename Traits::LocalInterpolationType& localInterpolation () const
571  {
572  return interpolation;
573  }
574 
577  GeometryType type () const
578  {
579  return gt;
580  }
581 
583  {
584  return new QkCGLocalFiniteElement(*this);
585  }
586 
587  private:
588  LocalBasis basis;
589  QkCGLocalCoefficients coefficients;
590  LocalInterpolation interpolation;
591  GeometryType gt;
592  };
593 
594 
598  template<class D, class R>
599  class QkCGLocalFiniteElement<D,R,2,3>
600  {
601  enum{ k = 2 };
602  enum{ d = 3 };
604 
605  class QkCGLocalCoefficients
606  {
607  enum{ k = 2 };
608  enum{ d = 3 };
609  public:
611  QkCGLocalCoefficients () : li(QkStuff::QkSize<k,d>::value)
612  {
613  li[24] = LocalKey(6,3,0); li[25] = LocalKey(11,2,0); li[26] = LocalKey(7,3,0);
614  li[21] = LocalKey(8,2,0); li[22] = LocalKey(5,1,0); li[23] = LocalKey(9,2,0);
615  li[18] = LocalKey(4,3,0); li[19] = LocalKey(10,2,0); li[20] = LocalKey(5,3,0);
616 
617  li[15] = LocalKey(2,2,0); li[16] = LocalKey(3,1,0); li[17] = LocalKey(3,2,0);
618  li[12] = LocalKey(0,1,0); li[13] = LocalKey(0,0,0); li[14] = LocalKey(1,1,0);
619  li[ 9] = LocalKey(0,2,0); li[10] = LocalKey(2,1,0); li[11] = LocalKey(1,2,0);
620 
621  li[ 6] = LocalKey(2,3,0); li[ 7] = LocalKey(7,2,0); li[ 8] = LocalKey(3,3,0);
622  li[ 3] = LocalKey(4,2,0); li[ 4] = LocalKey(4,1,0); li[ 5] = LocalKey(5,2,0);
623  li[ 0] = LocalKey(0,3,0); li[ 1] = LocalKey(6,2,0); li[ 2] = LocalKey(1,3,0);
624  }
625 
627  std::size_t size () const
628  {
630  }
631 
633  const LocalKey& localKey (std::size_t i) const
634  {
635  return li[i];
636  }
637 
638  private:
639  std::vector<LocalKey> li;
640  };
641 
643 
644  public:
645  // static number of basis functions
647 
650  typedef LocalFiniteElementTraits<LocalBasis,QkCGLocalCoefficients,LocalInterpolation> Traits;
651 
655  {
656  gt.makeCube(d);
657  }
658 
661  const typename Traits::LocalBasisType& localBasis () const
662  {
663  return basis;
664  }
665 
668  const typename Traits::LocalCoefficientsType& localCoefficients () const
669  {
670  return coefficients;
671  }
672 
675  const typename Traits::LocalInterpolationType& localInterpolation () const
676  {
677  return interpolation;
678  }
679 
682  GeometryType type () const
683  {
684  return gt;
685  }
686 
688  {
689  return new QkCGLocalFiniteElement(*this);
690  }
691 
692  private:
693  LocalBasis basis;
694  QkCGLocalCoefficients coefficients;
695  LocalInterpolation interpolation;
696  GeometryType gt;
697  };
698 
699 }
700 
701 namespace Dune {
702  namespace PDELab {
703 
706  template<class D, class R, int k, int d>
708  : public Dune::PDELab::SimpleLocalFiniteElementMap< Dune::QkDGLocalFiniteElement<D,R,k,d> >
709  {
710 
711  public:
712 
713  bool fixedSize() const
714  {
715  return true;
716  }
717 
718  std::size_t size(GeometryType gt) const
719  {
720  if (gt == GeometryType(GeometryType::cube,d))
722  else
723  return 0;
724  }
725 
726  std::size_t maxLocalSize() const
727  {
729  }
730 
731  };
732 
733  }
734 }
735 
736 #endif
std::size_t maxLocalSize() const
Definition: qkdg.hh:726
QkCGLocalFiniteElement()
Definition: qkdg.hh:549
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdg.hh:468
QkDGLocalFiniteElement * clone() const
Definition: qkdg.hh:348
QkDGLocalCoefficients()
Standard constructor.
Definition: qkdg.hh:226
GeometryType type() const
Definition: qkdg.hh:343
GeometryType type() const
Definition: qkdg.hh:682
Definition: qkdg.hh:32
Definition: qkdg.hh:299
QkCGLocalFiniteElement * clone() const
Definition: qkdg.hh:687
R x(int i)
Definition: qkdg.hh:118
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdg.hh:461
unsigned int size() const
number of shape functions
Definition: qkdg.hh:156
Definition: qkdg.hh:396
QkDGLocalFiniteElement()
Definition: qkdg.hh:315
QkCGLocalFiniteElement()
Definition: qkdg.hh:447
static const unsigned int value
Definition: gridfunctionspace/tags.hh:175
LocalFiniteElementTraits< LocalBasis, QkCGLocalCoefficients, LocalInterpolation > Traits
Definition: qkdg.hh:443
const Traits::LocalBasisType & localBasis() const
Definition: qkdg.hh:661
Layout map for Q1 elements.
Definition: qkdg.hh:222
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdg.hh:153
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdg.hh:329
QkDGFiniteElementFactory()
default constructor
Definition: qkdg.hh:384
bool fixedSize() const
Definition: qkdg.hh:713
const Traits::LocalBasisType & localBasis() const
Definition: qkdg.hh:556
std::size_t size() const
number of coefficients
Definition: qkdg.hh:233
LocalFiniteElementTraits< LocalBasis, QkCGLocalCoefficients, LocalInterpolation > Traits
Definition: qkdg.hh:545
GeometryType type() const
Definition: qkdg.hh:475
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdg.hh:563
Lagrange shape functions of order k on the reference cube.
Definition: qkdg.hh:149
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdg.hh:284
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdg.hh:182
R dp(int i, D x)
Definition: qkdg.hh:72
std::size_t size(GeometryType gt) const
Definition: qkdg.hh:718
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdg.hh:675
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdg.hh:668
const Traits::LocalBasisType & localBasis() const
Definition: qkdg.hh:322
Definition: adaptivity.hh:27
R dp(int i, D x) const
Definition: qkdg.hh:102
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:107
Dune::FieldVector< int, d > multiindex(int i)
Definition: qkdg.hh:125
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdg.hh:162
GeometryType type() const
Definition: qkdg.hh:577
R p(int i, D x)
Definition: qkdg.hh:62
QkCGLocalFiniteElement * clone() const
Definition: qkdg.hh:480
LocalFiniteElementTraits< LocalBasis, QkCGLocalCoefficients, LocalInterpolation > Traits
Definition: qkdg.hh:650
QkCGLocalFiniteElement * clone() const
Definition: qkdg.hh:582
Lagrange polynomials of degree k at equidistant points as a class.
Definition: qkdg.hh:89
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdg.hh:570
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdg.hh:311
const LocalKey & localKey(std::size_t i) const
get i'th index
Definition: qkdg.hh:239
R p(int i, D x) const
Definition: qkdg.hh:93
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdg.hh:336
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdg.hh:209
QkCGLocalFiniteElement()
Definition: qkdg.hh:654
Definition: qkdg.hh:29
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdg.hh:256
Factory for global-valued QkDG elements.
Definition: qkdg.hh:367
const Traits::LocalBasisType & localBasis() const
Definition: qkdg.hh:454