DUNE PDELab (git)

qkdglobatto.hh
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>
14#include <dune/pdelab/finiteelement/qkdglagrange.hh>
15
16namespace 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 {
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");
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:
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
265 enum{ n = QkStuff::QkSize<k,d>::value };
266
270
274 : gt(GeometryTypes::cube(d))
275 {}
276
279 const typename Traits::LocalBasisType& localBasis () const
280 {
281 return basis;
282 }
283
287 {
288 return coefficients;
289 }
290
294 {
295 return interpolation;
296 }
297
300 std::size_t size() const
301 {
302 return basis.size();
303 }
304
308 {
309 return gt;
310 }
311
312 QkDGGLLocalFiniteElement* clone () const
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>
333 QkDGGLLocalFiniteElement<
334 typename Geometry::ctype, RF, k, Geometry::mydimension
335 >,
336 Geometry
337 >
338 {
341 > LFE;
343
344 static const LFE lfe;
345
346 public:
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
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
@ cube
Cube element in any nonnegative dimension.
Definition: type.hh:122
GridImp::ctype ctype
define type used for coordinates in grid module
Definition: geometry.hh:100
static constexpr int mydimension
geometry dimension
Definition: geometry.hh:94
Default exception for dummy implementations.
Definition: exceptions.hh:263
Factory for global-valued QkDG elements.
Definition: qkdglobatto.hh:338
QkDGGLFiniteElementFactory()
default constructor
Definition: qkdglobatto.hh:348
Definition: qkdglobatto.hh:258
const Traits::LocalInterpolationType & localInterpolation() const
Definition: qkdglobatto.hh:293
LocalFiniteElementTraits< LocalBasis, LocalCoefficients, LocalInterpolation > Traits
Definition: qkdglobatto.hh:269
std::size_t size() const
Definition: qkdglobatto.hh:300
const Traits::LocalBasisType & localBasis() const
Definition: qkdglobatto.hh:279
GeometryType type() const
Definition: qkdglobatto.hh:307
const Traits::LocalCoefficientsType & localCoefficients() const
Definition: qkdglobatto.hh:286
QkDGGLLocalFiniteElement()
Definition: qkdglobatto.hh:273
Lagrange polynomials at Gauss-Lobatto points.
Definition: qkdglobatto.hh:24
Layout map for Q1 elements.
Definition: qkdglagrange.hh:232
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
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
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:243
Definition: qkdglobatto.hh:208
void interpolate(const F &ff, std::vector< C > &out) const
Local interpolation of a function.
Definition: qkdglobatto.hh:215
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
Factory for ScalarLocalToGlobalFiniteElementAdaptor objects.
Definition: localtoglobaladaptors.hh:244
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr GeometryType cube(unsigned int dim)
Returns a GeometryType representing a hypercube of dimension dim.
Definition: type.hh:462
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:279
@ GaussLobatto
Gauss-Lobatto rules.
Definition: quadraturerules.hh:176
Dune namespace.
Definition: alignedallocator.hh:13
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:35
D DomainType
domain type
Definition: localbasis.hh:43
traits helper struct
Definition: localfiniteelementtraits.hh:13
LB LocalBasisType
Definition: localfiniteelementtraits.hh:16
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:20
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:24
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 23, 23:29, 2024)