dune-pdelab  2.5-dev
qkdglobatto.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_QKDGLOBATTO_HH
5 #define DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_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 
18  template<class D, class R, int k>
20  {
21  R xi_gl[k+1];
22  R w_gl[k+1];
23 
24  public:
26  {
27  int matched_order=-1;
28  int matched_size=-1;
29  for (int order=1; order<=40; order++)
30  {
31  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
32  if (rule.size()==k+1)
33  {
34  matched_order = order;
35  matched_size = rule.size();
36  //std::cout << "GL: input order=" << order << " delivered=" << rule.order() << " size=" << rule.size()<< std::endl;
37  break;
38  }
39  }
40  if (matched_order<0) DUNE_THROW(Dune::Exception,"could not find Gauss Lobatto rule of appropriate size");
41  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
42  size_t count=0;
43  for (typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
44  {
45  size_t group=count/2;
46  size_t member=count%2;
47  size_t newj;
48  if (member==1) newj=group; else newj=k-group;
49  xi_gl[newj] = it->position()[0];
50  w_gl[newj] = it->weight();
51  count++;
52  }
53  for (int j=0; j<matched_size/2; j++)
54  if (xi_gl[j]>0.5)
55  {
56  R temp=xi_gl[j];
57  xi_gl[j] = xi_gl[k-j];
58  xi_gl[k-j] = temp;
59  temp=w_gl[j];
60  w_gl[j] = w_gl[k-j];
61  w_gl[k-j] = temp;
62  }
63  // for (int i=0; i<k+1; i++)
64  // std::cout << "i=" << i
65  // << " xi=" << xi_gl[i]
66  // << " w=" << w_gl[i]
67  // << std::endl;
68  }
69 
70  // ith Lagrange polynomial of degree k in one dimension
71  R p (int i, D x) const
72  {
73  R result(1.0);
74  for (int j=0; j<=k; j++)
75  if (j!=i) result *= (x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
76  return result;
77  }
78 
79  // derivative of ith Lagrange polynomial of degree k in one dimension
80  R dp (int i, D x) const
81  {
82  R result(0.0);
83 
84  for (int j=0; j<=k; j++)
85  if (j!=i)
86  {
87  R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
88  for (int l=0; l<=k; l++)
89  if (l!=i && l!=j) prod *= (x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
90  result += prod;
91  }
92  return result;
93  }
94 
95  // get ith Lagrange point
96  R x (int i) const
97  {
98  return xi_gl[i];
99  }
100 
101  // get weight of ith Lagrange point
102  R w (int i) const
103  {
104  return w_gl[i];
105  }
106  };
107 
120  template<class D, class R, int k, int d>
122  {
123  enum{ n = QkSize<k,d>::value };
125 
126  public:
127  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
128 
130  unsigned int size () const
131  {
132  return QkSize<k,d>::value;
133  }
134 
136  inline void evaluateFunction (const typename Traits::DomainType& in,
137  std::vector<typename Traits::RangeType>& out) const
138  {
139  out.resize(size());
140  for (size_t i=0; i<size(); i++)
141  {
142  // convert index i to multiindex
143  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
144 
145  // initialize product
146  out[i] = 1.0;
147 
148  // dimension by dimension
149  for (int j=0; j<d; j++)
150  out[i] *= poly.p(alpha[j],in[j]);
151  }
152  }
153 
155  inline void
156  evaluateJacobian (const typename Traits::DomainType& in, // position
157  std::vector<typename Traits::JacobianType>& out) const // return value
158  {
159  out.resize(size());
160 
161  // Loop over all shape functions
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  // Loop over all coordinate directions
168  for (int j=0; j<d; j++)
169  {
170  // Initialize: the overall expression is a product
171  // if j-th bit of i is set to -1, else 1
172  out[i][0][j] = poly.dp(alpha[j],in[j]);
173 
174  // rest of the product
175  for (int l=0; l<d; l++)
176  if (l!=j)
177  out[i][0][j] *= poly.p(alpha[l],in[l]);
178  }
179  }
180  }
181 
183  unsigned int order () const
184  {
185  return k;
186  }
187  };
188 
190  template<int k, int d, class LB>
192  {
194 
195  public:
196 
198  template<typename F, typename C>
199  void interpolate (const F& f, std::vector<C>& out) const
200  {
201  typename LB::Traits::DomainType x;
202  typename LB::Traits::RangeType y;
203 
204  out.resize(QkSize<k,d>::value);
205 
206  for (int i=0; i<QkSize<k,d>::value; i++)
207  {
208  // convert index i to multiindex
209  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
210 
211  // Generate coordinate of the i-th Lagrange point
212  for (int j=0; j<d; j++)
213  x[j] = poly.x(alpha[j]);
214 
215  f.evaluate(x,y); out[i] = y;
216  }
217  }
218  };
219 
221  template<int d, class LB>
223  {
224  public:
226  template<typename F, typename C>
227  void interpolate (const F& f, std::vector<C>& out) const
228  {
229  typename LB::Traits::DomainType x(0.5);
230  typename LB::Traits::RangeType y;
231  f.evaluate(x,y);
232  out.resize(1);
233  out[0] = y;
234  }
235  };
236 
237  }
238 
241  template<class D, class R, int k, int d>
243  {
247 
248  public:
249  // static number of basis functions
251 
254  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
255 
259  {
260  gt.makeCube(d);
261  }
262 
265  const typename Traits::LocalBasisType& localBasis () const
266  {
267  return basis;
268  }
269 
272  const typename Traits::LocalCoefficientsType& localCoefficients () const
273  {
274  return coefficients;
275  }
276 
279  const typename Traits::LocalInterpolationType& localInterpolation () const
280  {
281  return interpolation;
282  }
283 
286  GeometryType type () const
287  {
288  return gt;
289  }
290 
292  {
293  return new QkDGGLLocalFiniteElement(*this);
294  }
295 
296  private:
297  LocalBasis basis;
298  LocalCoefficients coefficients;
299  LocalInterpolation interpolation;
300  GeometryType gt;
301  };
302 
304 
309  template<class Geometry, class RF, int k>
311  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
312  QkDGGLLocalFiniteElement<
313  typename Geometry::ctype, RF, k, Geometry::mydimension
314  >,
315  Geometry
316  >
317  {
318  typedef QkDGGLLocalFiniteElement<
319  typename Geometry::ctype, RF, k, Geometry::mydimension
320  > LFE;
321  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
322 
323  static const LFE lfe;
324 
325  public:
327  QkDGGLFiniteElementFactory() : Base(lfe) {}
328  };
329 
330  template<class Geometry, class RF, int k>
333 
334 }
335 
336 #endif // DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:19
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglobatto.hh:127
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:254
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:227
unsigned int size() const
number of shape functions
Definition: qkdglobatto.hh:130
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:327
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglobatto.hh:156
Definition: qkdglagrange.hh:21
Definition: qkdglobatto.hh:191
R p(int i, D x) const
Definition: qkdglobatto.hh:71
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglobatto.hh:136
R x(int i) const
Definition: qkdglobatto.hh:96
Definition: qkdglobatto.hh:242
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:272
QkDGGLLocalFiniteElement * clone() const
Definition: qkdglobatto.hh:291
R dp(int i, D x) const
Definition: qkdglobatto.hh:80
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:199
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:265
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
R w(int i) const
Definition: qkdglobatto.hh:102
Layout map for Q1 elements.
Definition: qkdglagrange.hh:214
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:279
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:258
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:310
GaussLobattoLagrangePolynomials()
Definition: qkdglobatto.hh:25
GeometryType type() const
Definition: qkdglobatto.hh:286
Lagrange shape functions of order k on the reference cube.
Definition: qkdglobatto.hh:121
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglobatto.hh:183