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  : gt(GeometryTypes::cube(d))
260  {}
261 
264  const typename Traits::LocalBasisType& localBasis () const
265  {
266  return basis;
267  }
268 
271  const typename Traits::LocalCoefficientsType& localCoefficients () const
272  {
273  return coefficients;
274  }
275 
278  const typename Traits::LocalInterpolationType& localInterpolation () const
279  {
280  return interpolation;
281  }
282 
285  GeometryType type () const
286  {
287  return gt;
288  }
289 
291  {
292  return new QkDGGLLocalFiniteElement(*this);
293  }
294 
295  private:
296  LocalBasis basis;
297  LocalCoefficients coefficients;
298  LocalInterpolation interpolation;
299  GeometryType gt;
300  };
301 
303 
308  template<class Geometry, class RF, int k>
310  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
311  QkDGGLLocalFiniteElement<
312  typename Geometry::ctype, RF, k, Geometry::mydimension
313  >,
314  Geometry
315  >
316  {
317  typedef QkDGGLLocalFiniteElement<
318  typename Geometry::ctype, RF, k, Geometry::mydimension
319  > LFE;
320  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
321 
322  static const LFE lfe;
323 
324  public:
326  QkDGGLFiniteElementFactory() : Base(lfe) {}
327  };
328 
329  template<class Geometry, class RF, int k>
330  const typename QkDGGLFiniteElementFactory<Geometry, RF, k>::LFE
331  QkDGGLFiniteElementFactory<Geometry, RF, k>::lfe;
332 
333 }
334 
335 #endif // DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
Dune::QkStuff::QkDGLocalCoefficients
Layout map for Q1 elements.
Definition: qkdglagrange.hh:214
Dune::QkDGGLLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:264
Dune::QkDGGLFiniteElementFactory::QkDGGLFiniteElementFactory
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:326
Dune::QkDGGLLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:278
Dune
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Dune::QkDGGLFiniteElementFactory
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:309
Dune::QkDGGLLocalFiniteElement::clone
QkDGGLLocalFiniteElement * clone() const
Definition: qkdglobatto.hh:290
Dune::QkStuff::GaussLobattoLagrangePolynomials::dp
R dp(int i, D x) const
Definition: qkdglobatto.hh:80
Dune::QkStuff::GaussLobattoLagrangePolynomials
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:19
Dune::QkStuff::QkGLLocalBasis::size
unsigned int size() const
number of shape functions
Definition: qkdglobatto.hh:130
Dune::QkDGGLLocalFiniteElement::n
@ n
Definition: qkdglobatto.hh:250
Dune::QkDGGLLocalFiniteElement::type
GeometryType type() const
Definition: qkdglobatto.hh:285
Dune::QkDGGLLocalFiniteElement
Definition: qkdglobatto.hh:242
Dune::QkStuff::QkSize
Definition: qkdglagrange.hh:21
Dune::QkStuff::QkGLLocalInterpolation
Definition: qkdglobatto.hh:191
Dune::QkStuff::GaussLobattoLagrangePolynomials::p
R p(int i, D x) const
Definition: qkdglobatto.hh:71
Dune::QkDGGLLocalFiniteElement::QkDGGLLocalFiniteElement
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:258
Dune::QkStuff::GaussLobattoLagrangePolynomials::x
R x(int i) const
Definition: qkdglobatto.hh:96
Dune::QkStuff::QkGLLocalInterpolation::interpolate
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:199
Dune::QkStuff::QkGLLocalBasis
Lagrange shape functions of order k on the reference cube.
Definition: qkdglobatto.hh:121
Dune::QkStuff::QkGLLocalInterpolation< 0, d, LB >::interpolate
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:227
Dune::QkStuff::GaussLobattoLagrangePolynomials::GaussLobattoLagrangePolynomials
GaussLobattoLagrangePolynomials()
Definition: qkdglobatto.hh:25
Dune::QkStuff::QkGLLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglobatto.hh:136
value
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
Dune::QkStuff::QkGLLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglobatto.hh:156
Dune::QkDGGLLocalFiniteElement::Traits
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:254
Dune::QkStuff::QkGLLocalBasis::Traits
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglobatto.hh:127
Dune::QkStuff::QkGLLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglobatto.hh:183
Dune::QkStuff::GaussLobattoLagrangePolynomials::w
R w(int i) const
Definition: qkdglobatto.hh:102
Dune::QkDGGLLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:271