dune-pdelab 2.7-git
Loading...
Searching...
No Matches
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
15namespace 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 {
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:
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
For backward compatibility – Do not use this!
Definition adaptivity.hh:28
Definition qkdglagrange.hh:24
Layout map for Q1 elements.
Definition qkdglagrange.hh:230
Lagrange polynomials at Gauss-Lobatto points.
Definition qkdglobatto.hh:23
R w(int i) const
Definition qkdglobatto.hh:105
GaussLobattoLagrangePolynomials()
Definition qkdglobatto.hh:28
R dp(int i, D x) const
Definition qkdglobatto.hh:83
R p(int i, D x) const
Definition qkdglobatto.hh:74
R x(int i) const
Definition qkdglobatto.hh:99
Lagrange shape functions of order k on the reference cube.
Definition qkdglobatto.hh:125
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
unsigned int order() const
Polynomial order of the shape functions.
Definition qkdglobatto.hh:198
unsigned int size() const
number of shape functions
Definition qkdglobatto.hh:133
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition qkdglobatto.hh:139
LocalBasisTraits< D, d, Dune::FieldVector< D, d >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, d > > Traits
Definition qkdglobatto.hh:130
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition qkdglobatto.hh:159
Definition qkdglobatto.hh:207
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition qkdglobatto.hh:214
void interpolate(const F &f, std::vector< C > &out) const
Local interpolation of a function.
Definition qkdglobatto.hh:241
Definition qkdglobatto.hh:255
const Traits::LocalInterpolationType & localInterpolation() const
Definition qkdglobatto.hh:290
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition qkdglobatto.hh:266
std::size_t size() const
Definition qkdglobatto.hh:297
const Traits::LocalBasisType & localBasis() const
Definition qkdglobatto.hh:276
QkDGGLLocalFiniteElement * clone() const
Definition qkdglobatto.hh:309
GeometryType type() const
Definition qkdglobatto.hh:304
@ n
Definition qkdglobatto.hh:262
const Traits::LocalCoefficientsType & localCoefficients() const
Definition qkdglobatto.hh:283
QkDGGLLocalFiniteElement()
Definition qkdglobatto.hh:270
Factory for global-valued QkDG elements.
Definition qkdglobatto.hh:335
QkDGGLFiniteElementFactory()
default constructor
Definition qkdglobatto.hh:345