dune-localfunctions 3.0-git
pk2dlocalbasis.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#ifndef DUNE_PK2DLOCALBASIS_HH
4#define DUNE_PK2DLOCALBASIS_HH
5
6#include <numeric>
7
8#include <dune/common/fmatrix.hh>
9
11
12namespace Dune
13{
26 template<class D, class R, unsigned int k>
28 {
29 public:
30
32 enum {N = (k+1)*(k+2)/2};
33
37 enum {O = k};
38
39 typedef LocalBasisTraits<D,2,Dune::FieldVector<D,2>,R,1,Dune::FieldVector<R,1>,
40 Dune::FieldMatrix<R,1,2>, 2 > Traits;
41
44 {
45 for (unsigned int i=0; i<=k; i++)
46 pos_[i] = (1.0*i)/std::max(k,(unsigned int)1);
47 }
48
50 unsigned int size () const
51 {
52 return N;
53 }
54
56 inline void evaluateFunction (const typename Traits::DomainType& x,
57 std::vector<typename Traits::RangeType>& out) const
58 {
59 out.resize(N);
60 // specialization for k==0, not clear whether that is needed
61 if (k==0) {
62 out[0] = 1;
63 return;
64 }
65
66 int n=0;
67 for (unsigned int j=0; j<=k; j++)
68 for (unsigned int i=0; i<=k-j; i++)
69 {
70 out[n] = 1.0;
71 for (unsigned int alpha=0; alpha<i; alpha++)
72 out[n] *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
73 for (unsigned int beta=0; beta<j; beta++)
74 out[n] *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
75 for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
76 out[n] *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
77 n++;
78 }
79 }
80
82 inline void
83 evaluateJacobian (const typename Traits::DomainType& x, // position
84 std::vector<typename Traits::JacobianType>& out) const // return value
85 {
86 out.resize(N);
87
88 // specialization for k==0, not clear whether that is needed
89 if (k==0) {
90 out[0][0][0] = 0; out[0][0][1] = 0;
91 return;
92 }
93
94 int n=0;
95 for (unsigned int j=0; j<=k; j++)
96 for (unsigned int i=0; i<=k-j; i++)
97 {
98 // x_0 derivative
99 out[n][0][0] = 0.0;
100 R factor=1.0;
101 for (unsigned int beta=0; beta<j; beta++)
102 factor *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
103 for (unsigned int a=0; a<i; a++)
104 {
105 R product=factor;
106 for (unsigned int alpha=0; alpha<i; alpha++)
107 if (alpha==a)
108 product *= D(1)/(pos_[i]-pos_[alpha]);
109 else
110 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
111 for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
112 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
113 out[n][0][0] += product;
114 }
115 for (unsigned int c=i+j+1; c<=k; c++)
116 {
117 R product=factor;
118 for (unsigned int alpha=0; alpha<i; alpha++)
119 product *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
120 for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
121 if (gamma==c)
122 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
123 else
124 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
125 out[n][0][0] += product;
126 }
127
128 // x_1 derivative
129 out[n][0][1] = 0.0;
130 factor = 1.0;
131 for (unsigned int alpha=0; alpha<i; alpha++)
132 factor *= (x[0]-pos_[alpha])/(pos_[i]-pos_[alpha]);
133 for (unsigned int b=0; b<j; b++)
134 {
135 R product=factor;
136 for (unsigned int beta=0; beta<j; beta++)
137 if (beta==b)
138 product *= D(1)/(pos_[j]-pos_[beta]);
139 else
140 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
141 for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
142 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
143 out[n][0][1] += product;
144 }
145 for (unsigned int c=i+j+1; c<=k; c++)
146 {
147 R product=factor;
148 for (unsigned int beta=0; beta<j; beta++)
149 product *= (x[1]-pos_[beta])/(pos_[j]-pos_[beta]);
150 for (unsigned int gamma=i+j+1; gamma<=k; gamma++)
151 if (gamma==c)
152 product *= -D(1)/(pos_[gamma]-pos_[i]-pos_[j]);
153 else
154 product *= (pos_[gamma]-x[0]-x[1])/(pos_[gamma]-pos_[i]-pos_[j]);
155 out[n][0][1] += product;
156 }
157
158 n++;
159 }
160
161 }
162
168 void partial(const std::array<unsigned int,2>& order,
169 const typename Traits::DomainType& in,
170 std::vector<typename Traits::RangeType>& out) const
171 {
172 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
173
174 switch (totalOrder)
175 {
176 case 0:
177 evaluateFunction(in,out);
178 break;
179 case 1:
180 {
181 std::array<int,1> directions;
182 directions[0] = std::find(order.begin(), order.end(), 1)-order.begin();
183 evaluate<1>(directions, in, out);
184 break;
185 }
186 case 2:
187 {
188 std::array<int,2> directions;
189 unsigned int counter = 0;
190 auto nonconstOrder = order; // need a copy that I can modify
191 for (int i=0; i<2; i++)
192 {
193 while (nonconstOrder[i])
194 {
195 directions[counter++] = i;
196 nonconstOrder[i]--;
197 }
198 }
199
200 evaluate<2>(directions, in, out);
201 break;
202 }
203 default:
204 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
205 }
206 }
207
209 template<unsigned int dOrder> //order of derivative
210 inline void evaluate(const std::array<int,dOrder>& directions, //direction of derivative
211 const typename Traits::DomainType& in, //position
212 std::vector<typename Traits::RangeType>& out) const //return value
213 {
214 out.resize(N);
215
216 if (dOrder > Traits::diffOrder)
217 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
218
219 if (dOrder==0)
220 evaluateFunction(in, out);
221 else if (dOrder==1)
222 {
223 int n=0;
224 for (unsigned int j=0; j<=k; j++)
225 for (unsigned int i=0; i<=k-j; i++, n++)
226 {
227 out[n] = 0.0;
228 for (unsigned int no1=0; no1 < k; no1++)
229 {
230 R factor = lagrangianFactorDerivative(directions[0], no1, i, j, in);
231 for (unsigned int no2=0; no2 < k; no2++)
232 if (no1 != no2)
233 factor *= lagrangianFactor(no2, i, j, in);
234
235 out[n] += factor;
236 }
237 }
238 }
239 else if (dOrder==2)
240 {
241 // specialization for k<2, not clear whether that is needed
242 if (k<2) {
243 std::fill(out.begin(), out.end(), 0.0);
244 return;
245 }
246
247 //f = prod_{i} f_i -> dxa dxb f = sum_{i} {dxa f_i sum_{k \neq i} dxb f_k prod_{l \neq k,i} f_l
248 int n=0;
249 for (unsigned int j=0; j<=k; j++)
250 for (unsigned int i=0; i<=k-j; i++, n++)
251 {
252 R res = 0.0;
253
254 for (unsigned int no1=0; no1 < k; no1++)
255 {
256 R factor1 = lagrangianFactorDerivative(directions[0], no1, i, j, in);
257 for (unsigned int no2=0; no2 < k; no2++)
258 {
259 if (no1 == no2)
260 continue;
261 R factor2 = factor1*lagrangianFactorDerivative(directions[1], no2, i, j, in);
262 for (unsigned int no3=0; no3 < k; no3++)
263 {
264 if (no3 == no1 || no3 == no2)
265 continue;
266 factor2 *= lagrangianFactor(no3, i, j, in);
267 }
268 res += factor2;
269 }
270
271 }
272 out[n] = res;
273 }
274 }
275 }
276
278 unsigned int order () const
279 {
280 return k;
281 }
282
283 private:
285 typename Traits::RangeType lagrangianFactor(const int no, const int i, const int j, const typename Traits::DomainType& x) const
286 {
287 if ( no < i)
288 return (x[0]-pos_[no])/(pos_[i]-pos_[no]);
289 if (no < i+j)
290 return (x[1]-pos_[no-i])/(pos_[j]-pos_[no-i]);
291 return (pos_[no+1]-x[0]-x[1])/(pos_[no+1]-pos_[i]-pos_[j]);
292 }
293
297 typename Traits::RangeType lagrangianFactorDerivative(const int direction, const int no, const int i, const int j, const typename Traits::DomainType& x) const
298 {
299 if ( no < i)
300 return (direction == 0) ? 1.0/(pos_[i]-pos_[no]) : 0;
301
302 if (no < i+j)
303 return (direction == 0) ? 0: 1.0/(pos_[j]-pos_[no-i]);
304
305 return -1.0/(pos_[no+1]-pos_[i]-pos_[j]);
306 }
307
308 D pos_[k+1]; // positions on the interval
309 };
310
311}
312#endif
Definition brezzidouglasmarini1cube2d.hh:14
Type traits for LocalBasisVirtualInterface.
Definition localbasis.hh:38
R RangeType
range type
Definition localbasis.hh:61
@ diffOrder
number of partial derivatives supported
Definition localbasis.hh:74
D DomainType
domain type
Definition localbasis.hh:49
Lagrange shape functions of arbitrary order on the reference triangle.
Definition pk2dlocalbasis.hh:28
void evaluate(const std::array< int, dOrder > &directions, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate higher derivatives of all shape functions.
Definition pk2dlocalbasis.hh:210
@ O
Definition pk2dlocalbasis.hh:37
void evaluateFunction(const typename Traits::DomainType &x, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition pk2dlocalbasis.hh:56
unsigned int size() const
number of shape functions
Definition pk2dlocalbasis.hh:50
void partial(const std::array< unsigned int, 2 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition pk2dlocalbasis.hh:168
void evaluateJacobian(const typename Traits::DomainType &x, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition pk2dlocalbasis.hh:83
LocalBasisTraits< D, 2, Dune::FieldVector< D, 2 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 2 >, 2 > Traits
Definition pk2dlocalbasis.hh:40
@ N
Definition pk2dlocalbasis.hh:32
Pk2DLocalBasis()
Standard constructor.
Definition pk2dlocalbasis.hh:43
unsigned int order() const
Polynomial order of the shape functions.
Definition pk2dlocalbasis.hh:278