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