dune-pdelab  2.4.1
qkdggl.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 // Qk DG basis with Gauss Lobatto points
5 
6 #ifndef DUNE_QkDGGL_LOCALFINITEELEMENT_HH
7 #define DUNE_QkDGGL_LOCALFINITEELEMENT_HH
8 
9 #include <dune/common/fvector.hh>
10 #include <dune/common/deprecated.hh>
11 
12 #include <dune/geometry/type.hh>
13 #include <dune/geometry/quadraturerules.hh>
14 
15 #include <dune/localfunctions/common/localbasis.hh>
16 #include <dune/localfunctions/common/localfiniteelementtraits.hh>
17 #include <dune/localfunctions/common/localkey.hh>
18 #include <dune/localfunctions/common/localtoglobaladaptors.hh>
19 
21 
22 namespace Dune
23 {
24 
25  namespace QkStuff
26  {
27 
29  template<class D, class R, int k>
31  {
32  R xi_gl[k+1];
33  R w_gl[k+1];
34 
35  public:
37  {
38  int matched_order=-1;
39  int matched_size=-1;
40  for (int order=1; order<=40; order++)
41  {
42  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
43  if (rule.size()==k+1)
44  {
45  matched_order = order;
46  matched_size = rule.size();
47  //std::cout << "GL: input order=" << order << " delivered=" << rule.order() << " size=" << rule.size()<< std::endl;
48  break;
49  }
50  }
51  if (matched_order<0) DUNE_THROW(Dune::Exception,"could not find Gauss Lobatto rule of appropriate size");
52  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
53  size_t count=0;
54  for (typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
55  {
56  size_t group=count/2;
57  size_t member=count%2;
58  size_t newj;
59  if (member==1) newj=group; else newj=k-group;
60  xi_gl[newj] = it->position()[0];
61  w_gl[newj] = it->weight();
62  count++;
63  }
64  for (size_t j=0; j<matched_size/2; j++)
65  if (xi_gl[j]>0.5)
66  {
67  R temp=xi_gl[j];
68  xi_gl[j] = xi_gl[k-j];
69  xi_gl[k-j] = temp;
70  temp=w_gl[j];
71  w_gl[j] = w_gl[k-j];
72  w_gl[k-j] = temp;
73  }
74  // for (int i=0; i<k+1; i++)
75  // std::cout << "i=" << i
76  // << " xi=" << xi_gl[i]
77  // << " w=" << w_gl[i]
78  // << std::endl;
79  }
80 
81  // ith Lagrange polynomial of degree k in one dimension
82  R p (int i, D x) const
83  {
84  R result(1.0);
85  for (int j=0; j<=k; j++)
86  if (j!=i) result *= (x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
87  return result;
88  }
89 
90  // derivative of ith Lagrange polynomial of degree k in one dimension
91  R dp (int i, D x) const
92  {
93  R result(0.0);
94 
95  for (int j=0; j<=k; j++)
96  if (j!=i)
97  {
98  R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
99  for (int l=0; l<=k; l++)
100  if (l!=i && l!=j) prod *= (x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
101  result += prod;
102  }
103  return result;
104  }
105 
106  // get ith Lagrange point
107  R x (int i) const
108  {
109  return xi_gl[i];
110  }
111 
112  // get weight of ith Lagrange point
113  R w (int i) const
114  {
115  return w_gl[i];
116  }
117  };
118 
131  template<class D, class R, int k, int d>
133  {
134  enum{ n = QkSize<k,d>::value };
136 
137  public:
138  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
139 
141  unsigned int size () const
142  {
143  return QkSize<k,d>::value;
144  }
145 
147  inline void evaluateFunction (const typename Traits::DomainType& in,
148  std::vector<typename Traits::RangeType>& out) const
149  {
150  out.resize(size());
151  for (size_t i=0; i<size(); i++)
152  {
153  // convert index i to multiindex
154  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
155 
156  // initialize product
157  out[i] = 1.0;
158 
159  // dimension by dimension
160  for (int j=0; j<d; j++)
161  out[i] *= poly.p(alpha[j],in[j]);
162  }
163  }
164 
166  inline void
167  evaluateJacobian (const typename Traits::DomainType& in, // position
168  std::vector<typename Traits::JacobianType>& out) const // return value
169  {
170  out.resize(size());
171 
172  // Loop over all shape functions
173  for (size_t i=0; i<size(); i++)
174  {
175  // convert index i to multiindex
176  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
177 
178  // Loop over all coordinate directions
179  for (int j=0; j<d; j++)
180  {
181  // Initialize: the overall expression is a product
182  // if j-th bit of i is set to -1, else 1
183  out[i][0][j] = poly.dp(alpha[j],in[j]);
184 
185  // rest of the product
186  for (int l=0; l<d; l++)
187  if (l!=j)
188  out[i][0][j] *= poly.p(alpha[l],in[l]);
189  }
190  }
191  }
192 
194  unsigned int order () const
195  {
196  return k;
197  }
198  };
199 
201  template<int k, int d, class LB>
203  {
205 
206  public:
207 
209  template<typename F, typename C>
210  void interpolate (const F& f, std::vector<C>& out) const
211  {
212  typename LB::Traits::DomainType x;
213  typename LB::Traits::RangeType y;
214 
215  out.resize(QkSize<k,d>::value);
216 
217  for (int i=0; i<QkSize<k,d>::value; i++)
218  {
219  // convert index i to multiindex
220  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
221 
222  // Generate coordinate of the i-th Lagrange point
223  for (int j=0; j<d; j++)
224  x[j] = poly.x(alpha[j]);
225 
226  f.evaluate(x,y); out[i] = y;
227  }
228  }
229  };
230 
232  template<int d, class LB>
234  {
235  public:
237  template<typename F, typename C>
238  void interpolate (const F& f, std::vector<C>& out) const
239  {
240  typename LB::Traits::DomainType x(0.5);
241  typename LB::Traits::RangeType y;
242  f.evaluate(x,y);
243  out.resize(1);
244  out[0] = y;
245  }
246  };
247 
248  }
249 
252  template<class D, class R, int k, int d>
254  {
258 
259  public:
260  // static number of basis functions
262 
265  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
266 
270  {
271  gt.makeCube(d);
272  }
273 
276  const typename Traits::LocalBasisType& localBasis () const
277  {
278  return basis;
279  }
280 
283  const typename Traits::LocalCoefficientsType& localCoefficients () const
284  {
285  return coefficients;
286  }
287 
290  const typename Traits::LocalInterpolationType& localInterpolation () const
291  {
292  return interpolation;
293  }
294 
297  GeometryType type () const
298  {
299  return gt;
300  }
301 
303  {
304  return new QkDGGLLocalFiniteElement(*this);
305  }
306 
307  private:
308  LocalBasis basis;
309  LocalCoefficients coefficients;
310  LocalInterpolation interpolation;
311  GeometryType gt;
312  };
313 
315 
320  template<class Geometry, class RF, int k>
322  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
323  QkDGGLLocalFiniteElement<
324  typename Geometry::ctype, RF, k, Geometry::mydimension
325  >,
326  Geometry
327  >
328  {
329  typedef QkDGGLLocalFiniteElement<
330  typename Geometry::ctype, RF, k, Geometry::mydimension
331  > LFE;
332  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
333 
334  static const LFE lfe;
335 
336  public:
338  QkDGGLFiniteElementFactory() : Base(lfe) {}
339  };
340 
341 template<class Geometry, class RF, int k>
344 
345 }
346 
347 namespace Dune {
348  namespace PDELab {
349 
352  template<class D, class R, int k, int d>
354  : public Dune::PDELab::SimpleLocalFiniteElementMap< Dune::QkDGGLLocalFiniteElement<D,R,k,d> >
355  {
356 
357  public:
358 
359  bool fixedSize() const
360  {
361  return true;
362  }
363 
364  bool hasDOFs(int codim) const
365  {
366  return codim == 0;
367  }
368 
369  std::size_t size(GeometryType gt) const
370  {
371  if (gt == GeometryType(GeometryType::cube,d))
373  else
374  return 0;
375  }
376 
377  std::size_t maxLocalSize() const
378  {
380  }
381 
382  };
383 
384  }
385 }
386 
387 #endif
Definition: qkdg.hh:29
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdggl.hh:147
Definition: qkdggl.hh:253
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdggl.hh:265
Lagrange shape functions of order k on the reference cube.
Definition: qkdggl.hh:132
static const unsigned int value
Definition: gridfunctionspace/tags.hh:177
unsigned int size() const
number of shape functions
Definition: qkdggl.hh:141
R p(int i, D x) const
Definition: qkdggl.hh:82
GaussLobattoLagrangePolynomials()
Definition: qkdggl.hh:36
Definition: adaptivity.hh:27
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdggl.hh:290
simple implementation where all entities have the same finite element
Definition: finiteelementmap.hh:107
Factory for global-valued QkDG elements.
Definition: qkdggl.hh:321
Layout map for Q1 elements.
Definition: qkdg.hh:222
std::size_t size(GeometryType gt) const
Definition: qkdggl.hh:369
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdggl.hh:138
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdggl.hh:30
QkDGGLLocalFiniteElement * clone() const
Definition: qkdggl.hh:302
const Traits::LocalBasisType & localBasis() const
Definition: qkdggl.hh:276
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdggl.hh:210
R dp(int i, D x) const
Definition: qkdggl.hh:91
bool fixedSize() const
Definition: qkdggl.hh:359
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdggl.hh:338
GeometryType type() const
Definition: qkdggl.hh:297
bool hasDOFs(int codim) const
Definition: qkdggl.hh:364
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdggl.hh:167
R w(int i) const
Definition: qkdggl.hh:113
QkDGGLLocalFiniteElement()
Definition: qkdggl.hh:269
R x(int i) const
Definition: qkdggl.hh:107
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdggl.hh:238
Definition: qkdggl.hh:202
std::size_t maxLocalSize() const
Definition: qkdggl.hh:377
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdggl.hh:194
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdggl.hh:283