dune-pdelab  2.7-git
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 <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 #include <dune/localfunctions/common/localinterpolation.hh>
15 
16 namespace Dune
17 {
18  namespace QkStuff
19  {
20 
22  template<class D, class R, int k>
24  {
25  R xi_gl[k+1];
26  R w_gl[k+1];
27 
28  public:
30  {
31  int matched_order=-1;
32  int matched_size=-1;
33  for (int order=1; order<=40; order++)
34  {
35  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,order,Dune::QuadratureType::GaussLobatto);
36  if (rule.size()==k+1)
37  {
38  matched_order = order;
39  matched_size = rule.size();
40  //std::cout << "GL: input order=" << order << " delivered=" << rule.order() << " size=" << rule.size()<< std::endl;
41  break;
42  }
43  }
44  if (matched_order<0) DUNE_THROW(Dune::Exception,"could not find Gauss Lobatto rule of appropriate size");
45  const Dune::QuadratureRule<D,1>& rule = Dune::QuadratureRules<D,1>::rule(Dune::GeometryType::cube,matched_order,Dune::QuadratureType::GaussLobatto);
46  size_t count=0;
47  for (typename Dune::QuadratureRule<D,1>::const_iterator it=rule.begin(); it!=rule.end(); ++it)
48  {
49  size_t group=count/2;
50  size_t member=count%2;
51  size_t newj;
52  if (member==1) newj=group; else newj=k-group;
53  xi_gl[newj] = it->position()[0];
54  w_gl[newj] = it->weight();
55  count++;
56  }
57  for (int j=0; j<matched_size/2; j++)
58  if (xi_gl[j]>0.5)
59  {
60  R temp=xi_gl[j];
61  xi_gl[j] = xi_gl[k-j];
62  xi_gl[k-j] = temp;
63  temp=w_gl[j];
64  w_gl[j] = w_gl[k-j];
65  w_gl[k-j] = temp;
66  }
67  // for (int i=0; i<k+1; i++)
68  // std::cout << "i=" << i
69  // << " xi=" << xi_gl[i]
70  // << " w=" << w_gl[i]
71  // << std::endl;
72  }
73 
74  // ith Lagrange polynomial of degree k in one dimension
75  R p (int i, D x) const
76  {
77  R result(1.0);
78  for (int j=0; j<=k; j++)
79  if (j!=i) result *= (x-xi_gl[j])/(xi_gl[i]-xi_gl[j]);
80  return result;
81  }
82 
83  // derivative of ith Lagrange polynomial of degree k in one dimension
84  R dp (int i, D x) const
85  {
86  R result(0.0);
87 
88  for (int j=0; j<=k; j++)
89  if (j!=i)
90  {
91  R prod( 1.0/(xi_gl[i]-xi_gl[j]) );
92  for (int l=0; l<=k; l++)
93  if (l!=i && l!=j) prod *= (x-xi_gl[l])/(xi_gl[i]-xi_gl[l]);
94  result += prod;
95  }
96  return result;
97  }
98 
99  // get ith Lagrange point
100  R x (int i) const
101  {
102  return xi_gl[i];
103  }
104 
105  // get weight of ith Lagrange point
106  R w (int i) const
107  {
108  return w_gl[i];
109  }
110  };
111 
124  template<class D, class R, int k, int d>
126  {
127  enum{ n = QkSize<k,d>::value };
129 
130  public:
131  typedef LocalBasisTraits<D,d,Dune::FieldVector<D,d>,R,1,Dune::FieldVector<R,1>,Dune::FieldMatrix<R,1,d> > Traits;
132 
134  unsigned int size () const
135  {
136  return QkSize<k,d>::value;
137  }
138 
140  inline void evaluateFunction (const typename Traits::DomainType& in,
141  std::vector<typename Traits::RangeType>& out) const
142  {
143  out.resize(size());
144  for (size_t i=0; i<size(); i++)
145  {
146  // convert index i to multiindex
147  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
148 
149  // initialize product
150  out[i] = 1.0;
151 
152  // dimension by dimension
153  for (int j=0; j<d; j++)
154  out[i] *= poly.p(alpha[j],in[j]);
155  }
156  }
157 
159  inline void
160  evaluateJacobian (const typename Traits::DomainType& in, // position
161  std::vector<typename Traits::JacobianType>& out) const // return value
162  {
163  out.resize(size());
164 
165  // Loop over all shape functions
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  // Loop over all coordinate directions
172  for (int j=0; j<d; j++)
173  {
174  // Initialize: the overall expression is a product
175  // if j-th bit of i is set to -1, else 1
176  out[i][0][j] = poly.dp(alpha[j],in[j]);
177 
178  // rest of the product
179  for (int l=0; l<d; l++)
180  if (l!=j)
181  out[i][0][j] *= poly.p(alpha[l],in[l]);
182  }
183  }
184  }
185 
187  void partial(const std::array<unsigned int, Traits::dimDomain>& order,
188  const typename Traits::DomainType& in,
189  std::vector<typename Traits::RangeType>& out) const {
190  auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
191  if (totalOrder == 0) {
192  evaluateFunction(in, out);
193  } else {
194  DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
195  }
196  }
197 
199  unsigned int order () const
200  {
201  return k;
202  }
203  };
204 
206  template<int k, int d, class LB>
208  {
210 
211  public:
212 
214  template<typename F, typename C>
215  void interpolate (const F& ff, std::vector<C>& out) const
216  {
217  typename LB::Traits::DomainType x;
218  auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
219 
220  out.resize(QkSize<k,d>::value);
221 
222  for (int i=0; i<QkSize<k,d>::value; i++)
223  {
224  // convert index i to multiindex
225  Dune::FieldVector<int,d> alpha(multiindex<k,d>(i));
226 
227  // Generate coordinate of the i-th Lagrange point
228  for (int j=0; j<d; j++)
229  x[j] = poly.x(alpha[j]);
230 
231  out[i] = f(x);
232  }
233  }
234  };
235 
237  template<int d, class LB>
239  {
240  public:
242  template<typename F, typename C>
243  void interpolate (const F& ff, std::vector<C>& out) const
244  {
245  typename LB::Traits::DomainType x(0.5);
246  auto&& f = Impl::makeFunctionWithCallOperator<typename LB::Traits::DomainType>(ff);
247  out.resize(1);
248  out[0] = f(x);
249  }
250  };
251 
252  }
253 
256  template<class D, class R, int k, int d>
258  {
262 
263  public:
264  // static number of basis functions
266 
269  typedef LocalFiniteElementTraits<LocalBasis,LocalCoefficients,LocalInterpolation> Traits;
270 
274  : gt(GeometryTypes::cube(d))
275  {}
276 
279  const typename Traits::LocalBasisType& localBasis () const
280  {
281  return basis;
282  }
283 
286  const typename Traits::LocalCoefficientsType& localCoefficients () const
287  {
288  return coefficients;
289  }
290 
293  const typename Traits::LocalInterpolationType& localInterpolation () const
294  {
295  return interpolation;
296  }
297 
300  std::size_t size() const
301  {
302  return basis.size();
303  }
304 
307  GeometryType type () const
308  {
309  return gt;
310  }
311 
313  {
314  return new QkDGGLLocalFiniteElement(*this);
315  }
316 
317  private:
318  LocalBasis basis;
319  LocalCoefficients coefficients;
320  LocalInterpolation interpolation;
321  GeometryType gt;
322  };
323 
325 
330  template<class Geometry, class RF, int k>
332  public ScalarLocalToGlobalFiniteElementAdaptorFactory<
333  QkDGGLLocalFiniteElement<
334  typename Geometry::ctype, RF, k, Geometry::mydimension
335  >,
336  Geometry
337  >
338  {
339  typedef QkDGGLLocalFiniteElement<
340  typename Geometry::ctype, RF, k, Geometry::mydimension
341  > LFE;
342  typedef ScalarLocalToGlobalFiniteElementAdaptorFactory<LFE, Geometry> Base;
343 
344  static const LFE lfe;
345 
346  public:
348  QkDGGLFiniteElementFactory() : Base(lfe) {}
349  };
350 
351  template<class Geometry, class RF, int k>
352  const typename QkDGGLFiniteElementFactory<Geometry, RF, k>::LFE
353  QkDGGLFiniteElementFactory<Geometry, RF, k>::lfe;
354 
355 }
356 
357 #endif // DUNE_PDELAB_FINITEELEMENT_QKDGLOBATTO_HH
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
Definition: qkdglagrange.hh:26
Layout map for Q1 elements.
Definition: qkdglagrange.hh:232
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:24
R w(int i) const
Definition: qkdglobatto.hh:106
GaussLobattoLagrangePolynomials()
Definition: qkdglobatto.hh:29
R dp(int i, D x) const
Definition: qkdglobatto.hh:84
R p(int i, D x) const
Definition: qkdglobatto.hh:75
R x(int i) const
Definition: qkdglobatto.hh:100
Lagrange shape functions of order k on the reference cube.
Definition: qkdglobatto.hh:126
unsigned int order() const
Polynomial order of the shape functions.
Definition: qkdglobatto.hh:199
unsigned int size() const
number of shape functions
Definition: qkdglobatto.hh:134
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: qkdglobatto.hh:140
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition: qkdglobatto.hh:131
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: qkdglobatto.hh:160
void partial(const std::array< unsigned int, Traits::dimDomain > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivative of all shape functions.
Definition: qkdglobatto.hh:187
Definition: qkdglobatto.hh:208
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:215
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:243
Definition: qkdglobatto.hh:258
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:269
std::size_t size() const
Definition: qkdglobatto.hh:300
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:286
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:293
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:279
GeometryType type() const
Definition: qkdglobatto.hh:307
@ n
Definition: qkdglobatto.hh:265
QkDGGLLocalFiniteElement * clone() const
Definition: qkdglobatto.hh:312
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:273
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:338
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:348
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139