Dune Core Modules (2.6.0)

p23dlocalbasis.hh
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_P2_3DLOCALBASIS_HH
4#define DUNE_P2_3DLOCALBASIS_HH
5
6#include <numeric>
7
8#include <dune/common/fmatrix.hh>
9
10#include <dune/localfunctions/common/localbasis.hh>
11
12namespace Dune
13{
24 template<class D, class R>
25 class P23DLocalBasis
26 {
27 public:
29 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
30 Dune::FieldMatrix<R,1,3> > Traits;
31
33 unsigned int size () const
34 {
35 return 10;
36 }
37
39 inline void evaluateFunction (const typename Traits::DomainType& in,
40 std::vector<typename Traits::RangeType>& out) const
41 {
42 out.resize(10);
43
44 int coeff;
45 R a[2], b[3], c[3];
46
47 // case 0:
48 coeff=2;
49 a[0]=1.0;
50 a[1]=0.5;
51 b[0]=-1.0;
52 b[1]=-1.0;
53 b[2]=-1.0;
54 c[0]=-1.0;
55 c[1]=-1.0;
56 c[2]=-1.0;
57
58 out[0] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
59
60 // case 1:
61 coeff=2;
62 a[0]=0.0;
63 a[1]=-0.5;
64 b[0]=1.0;
65 b[1]=0.0;
66 b[2]=0.0;
67 c[0]=1.0;
68 c[1]=0.0;
69 c[2]=0.0;
70
71 out[1] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
72
73 // case 2:
74 coeff=2;
75 a[0]=0.0;
76 a[1]=-0.5;
77 b[0]=0.0;
78 b[1]=1.0;
79 b[2]=0.0;
80 c[0]=0.0;
81 c[1]=1.0;
82 c[2]=0.0;
83
84 out[2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
85
86 // case 3:
87 coeff=2;
88 a[0]=0.0;
89 a[1]=-0.5;
90 b[0]=0.0;
91 b[1]=0.0;
92 b[2]=1.0;
93 c[0]=0.0;
94 c[1]=0.0;
95 c[2]=1.0;
96
97 out[3] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
98
99 // case 4:
100 coeff=4;
101 a[0]=0.0;
102 a[1]=1.0;
103 b[0]=1.0;
104 b[1]=0.0;
105 b[2]=0.0;
106 c[0]=-1.0;
107 c[1]=-1.0;
108 c[2]=-1.0;
109
110 out[4] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
111
112 // case 5:
113 coeff=4;
114 a[0]=0.0;
115 a[1]=0.0;
116 b[0]=1.0;
117 b[1]=0.0;
118 b[2]=0.0;
119 c[0]=0.0;
120 c[1]=1.0;
121 c[2]=0.0;
122
123 out[5] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
124
125 // case 6:
126 coeff=4;
127 a[0]=0.0;
128 a[1]=1.0;
129 b[0]=0.0;
130 b[1]=1.0;
131 b[2]=0.0;
132 c[0]=-1.0;
133 c[1]=-1.0;
134 c[2]=-1.0;
135
136 out[6] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
137
138 // case 7:
139 coeff=4;
140 a[0]=0.0;
141 a[1]=1.0;
142 b[0]=0.0;
143 b[1]=0.0;
144 b[2]=1.0;
145 c[0]=-1.0;
146 c[1]=-1.0;
147 c[2]=-1.0;
148
149 out[7] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
150
151 // case 8:
152 coeff=4;
153 a[0]=0.0;
154 a[1]=0.0;
155 b[0]=1.0;
156 b[1]=0.0;
157 b[2]=0.0;
158 c[0]=0.0;
159 c[1]=0.0;
160 c[2]=1.0;
161
162 out[8] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
163
164 // case 9:
165 coeff=4;
166 a[0]=0.0;
167 a[1]=0.0;
168 b[0]=0.0;
169 b[1]=1.0;
170 b[2]=0.0;
171 c[0]=0.0;
172 c[1]=0.0;
173 c[2]=1.0;
174
175 out[9] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1] + b[2]*in[2]) * (a[1] + c[0]*in[0] + c[1]*in[1] + c[2]*in[2]);
176
177 }
178
180 inline void
181 evaluateJacobian (const typename Traits::DomainType& in, // position
182 std::vector<typename Traits::JacobianType>& out) const // return value
183 {
184 out.resize(10);
185
186 R aa[3][3], bb[3][3];
187 // case 0:
188 //x derivative
189 aa[0][0]=-3.0;
190 bb[0][0]=4.0;
191 bb[1][0]=4.0;
192 bb[2][0]=4.0;
193 //y derivative
194 aa[0][1]=-3.0;
195 bb[0][1]=4.0;
196 bb[1][1]=4.0;
197 bb[2][1]=4.0;
198 // z derivative
199 aa[0][2]=-3.0;
200 bb[0][2]=4.0;
201 bb[1][2]=4.0;
202 bb[2][2]=4.0;
203
204 out[0][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
205 out[0][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
206 out[0][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
207
208 //case 1:
209 //x derivative
210 aa[0][0]=-1.0;
211 bb[0][0]=4.0;
212 bb[1][0]=0.0;
213 bb[2][0]=0.0;
214 //y derivative
215 aa[0][1]=0.0;
216 bb[0][1]=0.0;
217 bb[1][1]=0.0;
218 bb[2][1]=0.0;
219 // z derivative
220 aa[0][2]=0.0;
221 bb[0][2]=0.0;
222 bb[1][2]=0.0;
223 bb[2][2]=0.0;
224
225 out[1][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
226 out[1][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
227 out[1][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
228
229 // case 2:
230 //x derivative
231 aa[0][0]=0.0;
232 bb[0][0]=0.0;
233 bb[1][0]=0.0;
234 bb[2][0]=0.0;
235 //y derivative
236 aa[0][1]=-1.0;
237 bb[0][1]=0.0;
238 bb[1][1]=4.0;
239 bb[2][1]=0.0;
240 // z derivative
241 aa[0][2]=0.0;
242 bb[0][2]=0.0;
243 bb[1][2]=0.0;
244 bb[2][2]=0.0;
245
246 out[2][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
247 out[2][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
248 out[2][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
249
250 // case 3:
251 //x derivative
252 aa[0][0]=0.0;
253 bb[0][0]=0.0;
254 bb[1][0]=0.0;
255 bb[2][0]=0.0;
256 //y derivative
257 aa[0][1]=0.0;
258 bb[0][1]=0.0;
259 bb[1][1]=0.0;
260 bb[2][1]=0.0;
261 // z derivative
262 aa[0][2]=-1.0;
263 bb[0][2]=0.0;
264 bb[1][2]=0.0;
265 bb[2][2]=4.0;
266
267 out[3][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
268 out[3][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
269 out[3][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
270
271 // case 4:
272 //x derivative
273 aa[0][0]=4.0;
274 bb[0][0]=-8.0;
275 bb[1][0]=-4.0;
276 bb[2][0]=-4.0;
277 //y derivative
278 aa[0][1]=0.0;
279 bb[0][1]=-4.0;
280 bb[1][1]=0.0;
281 bb[2][1]=0.0;
282 // z derivative
283 aa[0][2]=0.0;
284 bb[0][2]=-4.0;
285 bb[1][2]=0.0;
286 bb[2][2]=0.0;
287
288 out[4][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
289 out[4][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
290 out[4][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
291
292 // case 5:
293 //x derivative
294 aa[0][0]=0.0;
295 bb[0][0]=0.0;
296 bb[1][0]=4.0;
297 bb[2][0]=0.0;
298 //y derivative
299 aa[0][1]=0.0;
300 bb[0][1]=4.0;
301 bb[1][1]=0.0;
302 bb[2][1]=0.0;
303 // z derivative
304 aa[0][2]=0.0;
305 bb[0][2]=0.0;
306 bb[1][2]=0.0;
307 bb[2][2]=0.0;
308
309 out[5][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
310 out[5][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
311 out[5][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
312
313 // case 6:
314 //x derivative
315 aa[0][0]=0.0;
316 bb[0][0]=0.0;
317 bb[1][0]=-4.0;
318 bb[2][0]=0.0;
319 //y derivative
320 aa[0][1]=4.0;
321 bb[0][1]=-4.0;
322 bb[1][1]=-8.0;
323 bb[2][1]=-4.0;
324 // z derivative
325 aa[0][2]=0.0;
326 bb[0][2]=0.0;
327 bb[1][2]=-4.0;
328 bb[2][2]=0.0;
329
330 out[6][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
331 out[6][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
332 out[6][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
333
334 // case 7:
335 //x derivative
336 aa[0][0]=0.0;
337 bb[0][0]=0.0;
338 bb[1][0]=0.0;
339 bb[2][0]=-4.0;
340 //y derivative
341 aa[0][1]=0.0;
342 bb[0][1]=0.0;
343 bb[1][1]=0.0;
344 bb[2][1]=-4.0;
345 // z derivative
346 aa[0][2]=4.0;
347 bb[0][2]=-4.0;
348 bb[1][2]=-4.0;
349 bb[2][2]=-8.0;
350
351 out[7][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
352 out[7][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
353 out[7][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
354
355 //case 8:
356 //x derivative
357 aa[0][0]=0.0;
358 bb[0][0]=0.0;
359 bb[1][0]=0.0;
360 bb[2][0]=4.0;
361 //y derivative
362 aa[0][1]=0.0;
363 bb[0][1]=0.0;
364 bb[1][1]=0.0;
365 bb[2][1]=0.0;
366 // z derivative
367 aa[0][2]=0.0;
368 bb[0][2]=4.0;
369 bb[1][2]=0.0;
370 bb[2][2]=0.0;
371
372 out[8][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
373 out[8][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
374 out[8][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
375
376 // case 9:
377 //x derivative
378 aa[0][0]=0.0;
379 bb[0][0]=0.0;
380 bb[1][0]=0.0;
381 bb[2][0]=0.0;
382 //y derivative
383 aa[0][1]=0.0;
384 bb[0][1]=0.0;
385 bb[1][1]=0.0;
386 bb[2][1]=4.0;
387 // z derivative
388 aa[0][2]=0.0;
389 bb[0][2]=0.0;
390 bb[1][2]=4.0;
391 bb[2][2]=0.0;
392
393 out[9][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
394 out[9][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
395 out[9][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
396
397 }
398
404 void partial(const std::array<unsigned int,3>& order,
405 const typename Traits::DomainType& in,
406 std::vector<typename Traits::RangeType>& out) const
407 {
408 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
409 if (totalOrder == 0) {
410 evaluateFunction(in, out);
411 } else if (totalOrder == 1) {
412 out.resize(size());
413
414 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
415 switch (direction) {
416 case 0:
417 out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
418 out[1] =-1.0 + 4.0*in[0];
419 out[2] = 0.0;
420 out[3] = 0.0;
421 out[4] = 4.0 - 4.0*(2.0*in[0] + in[1] + in[2]);
422 out[5] = 4.0*in[1];
423 out[6] =-4.0*in[1];
424 out[7] =-4.0*in[2];
425 out[8] = 4.0*in[2];
426 out[9] = 0.0;
427 break;
428 case 1:
429 out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
430 out[1] = 0.0;
431 out[2] =-1.0 + 4.0*in[1];
432 out[3] = 0.0;
433 out[4] =-4.0*in[0];
434 out[5] = 4.0*in[0];
435 out[6] = 4.0 - 4.0*(in[0] + 2.0*in[1] + in[2]);
436 out[7] =-4.0*in[2];
437 out[8] = 0.0;
438 out[9] = 4.0*in[2];
439 break;
440 case 2:
441 out[0] =-3.0 + 4.0*(in[0] + in[1] + in[2]);
442 out[1] = 0.0;
443 out[2] = 0.0;
444 out[3] =-1.0 + 4.0*in[2];
445 out[4] =-4.0*in[0];
446 out[5] = 0.0;
447 out[6] =-4.0*in[1];
448 out[7] = 4.0 - 4.0*(in[0] + in[1] + 2.0*in[2]);
449 out[8] = 4.0*in[0];
450 out[9] = 4.0*in[1];
451 break;
452 default:
453 DUNE_THROW(RangeError, "Component out of range.");
454 }
455 } else {
456 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
457 }
458 }
459
461 unsigned int order () const
462 {
463 return 2;
464 }
465 };
466}
467#endif
Dune::FieldMatrix
A dense n x m matrix.
Definition: fmatrix.hh:68
Dune::FieldVector
vector space out of a tensor product of fields.
Definition: fvector.hh:93
Dune::NotImplemented
Default exception for dummy implementations.
Definition: exceptions.hh:261
Dune::P23DLocalBasis
Quadratic Lagrange shape functions on the tetrahedron.
Definition: p23dlocalbasis.hh:26
Dune::P23DLocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: p23dlocalbasis.hh:461
Dune::P23DLocalBasis::size
unsigned int size() const
number of shape functions
Definition: p23dlocalbasis.hh:33
Dune::P23DLocalBasis::Traits
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
export type traits for function signature
Definition: p23dlocalbasis.hh:30
Dune::P23DLocalBasis::partial
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of any order of all shape functions.
Definition: p23dlocalbasis.hh:404
Dune::P23DLocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: p23dlocalbasis.hh:181
Dune::P23DLocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: p23dlocalbasis.hh:39
Dune::RangeError
Default exception class for range errors.
Definition: exceptions.hh:252
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:331
Dune
Dune namespace.
Definition: alignedallocator.hh:10
Dune::LocalBasisTraits
Type traits for LocalBasisVirtualInterface.
Definition: localbasis.hh:32
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)