Dune Core Modules (2.6.0)

prismp2localbasis.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_PRISM_P2_LOCALBASIS_HH
4#define DUNE_PRISM_P2_LOCALBASIS_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 PrismP2LocalBasis
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 18;
36 }
37
39 inline void evaluateFunction (const typename Traits::DomainType& in,
40 std::vector<typename Traits::RangeType>& out) const
41 {
42 out.resize(18);
43
44
45 int coeff;
46 R a[2], b[2], c[2], a1d, b1d, c1d;
47
48
49 // lower triangle:
50 coeff= 2;
51 a[0] = 1;
52 a[1] = 0.5;
53 b[0] = -1;
54 b[1] = -1;
55 c[0] = -1;
56 c[1] = -1;
57 a1d = 1;
58 b1d = -3;
59 c1d = 2;
60 out[0] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
61
62
63 coeff= 2;
64 a[0] = 0;
65 a[1] = -0.5;
66 b[0] = 1;
67 b[1] = 0;
68 c[0] = 1;
69 c[1] = 0;
70 a1d = 1;
71 b1d = -3;
72 c1d = 2;
73 out[1] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
74
75
76 coeff= 2;
77 a[0] = 0;
78 a[1] = -0.5;
79 b[0] = 0;
80 b[1] = 1;
81 c[0] = 0;
82 c[1] = 1;
83 a1d = 1;
84 b1d = -3;
85 c1d = 2;
86 out[2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
87
88 //upper triangle
89 coeff= 2;
90 a[0] = 1;
91 a[1] = 0.5;
92 b[0] = -1;
93 b[1] = -1;
94 c[0] = -1;
95 c[1] = -1;
96 a1d = 0;
97 b1d = -1;
98 c1d = 2;
99 out[3] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
100
101
102 coeff= 2;
103 a[0] = 0;
104 a[1] = -0.5;
105 b[0] = 1;
106 b[1] = 0;
107 c[0] = 1;
108 c[1] = 0;
109 a1d = 0;
110 b1d = -1;
111 c1d = 2;
112 out[4] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
113
114
115 coeff= 2;
116 a[0] = 0;
117 a[1] = -0.5;
118 b[0] = 0;
119 b[1] = 1;
120 c[0] = 0;
121 c[1] = 1;
122 a1d = 0;
123 b1d = -1;
124 c1d = 2;
125 out[5] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
126
127 // vertical edges
128 coeff= 2;
129 a[0] = 1;
130 a[1] = 0.5;
131 b[0] = -1;
132 b[1] = -1;
133 c[0] = -1;
134 c[1] = -1;
135 a1d = 0;
136 b1d = 4;
137 c1d = -4;
138 out[6] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
139
140
141 coeff= 2;
142 a[0] = 0;
143 a[1] = -0.5;
144 b[0] = 1;
145 b[1] = 0;
146 c[0] = 1;
147 c[1] = 0;
148 a1d = 0;
149 b1d = 4;
150 c1d = -4;
151 out[7] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
152
153
154 coeff= 2;
155 a[0] = 0;
156 a[1] = -0.5;
157 b[0] = 0;
158 b[1] = 1;
159 c[0] = 0;
160 c[1] = 1;
161 a1d = 0;
162 b1d = 4;
163 c1d = -4;
164 out[8] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
165
166 // lower triangle edges
167 coeff= 4;
168 a[0] = 0;
169 a[1] = 1;
170 b[0] = 1;
171 b[1] = 0;
172 c[0] = -1;
173 c[1] = -1;
174 a1d = 1;
175 b1d = -3;
176 c1d = 2;
177 out[9] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
178
179
180 coeff= 4;
181 a[0] = 0;
182 a[1] = 1;
183 b[0] = 0;
184 b[1] = 1;
185 c[0] = -1;
186 c[1] = -1;
187 a1d = 1;
188 b1d = -3;
189 c1d = 2;
190 out[10] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
191
192
193 coeff= 4;
194 a[0] = 0;
195 a[1] = 0;
196 b[0] = 1;
197 b[1] = 0;
198 c[0] = 0;
199 c[1] = 1;
200 a1d = 1;
201 b1d = -3;
202 c1d = 2;
203 out[11] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
204
205 // upper triangle edges
206 coeff= 4;
207 a[0] = 0;
208 a[1] = 1;
209 b[0] = 1;
210 b[1] = 0;
211 c[0] = -1;
212 c[1] = -1;
213 a1d = 0;
214 b1d = -1;
215 c1d = 2;
216 out[12] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
217
218
219 coeff= 4;
220 a[0] = 0;
221 a[1] = 1;
222 b[0] = 0;
223 b[1] = 1;
224 c[0] = -1;
225 c[1] = -1;
226 a1d = 0;
227 b1d = -1;
228 c1d = 2;
229 out[13] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
230
231
232 coeff= 4;
233 a[0] = 0;
234 a[1] = 0;
235 b[0] = 1;
236 b[1] = 0;
237 c[0] = 0;
238 c[1] = 1;
239 a1d = 0;
240 b1d = -1;
241 c1d = 2;
242 out[14] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
243
244 // quadrilateral sides
245 coeff= 4;
246 a[0] = 0;
247 a[1] = 1;
248 b[0] = 1;
249 b[1] = 0;
250 c[0] = -1;
251 c[1] = -1;
252 a1d = 0;
253 b1d = 4;
254 c1d = -4;
255 out[15] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
256
257
258 coeff= 4;
259 a[0] = 0;
260 a[1] = 1;
261 b[0] = 0;
262 b[1] = 1;
263 c[0] = -1;
264 c[1] = -1;
265 a1d = 0;
266 b1d = 4;
267 c1d = -4;
268 out[16] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
269
270
271 coeff= 4;
272 a[0] = 0;
273 a[1] = 0;
274 b[0] = 1;
275 b[1] = 0;
276 c[0] = 0;
277 c[1] = 1;
278 a1d = 0;
279 b1d = 4;
280 c1d = -4;
281 out[17] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1])*(a1d + in[2]*b1d + in[2]*in[2]*c1d);
282 }
283
285 inline void
286 evaluateJacobian (const typename Traits::DomainType& in, // position
287 std::vector<typename Traits::JacobianType>& out) const // return value
288 {
289 out.resize(18);
290
291
292
293 int coeff;
294 R a[2], b[2], c[2], aa[2], bb[2][2], a1d, b1d, c1d;
295
296
297 // lower triangle:
298 coeff= 2;
299 a[0] = 1;
300 a[1] = 0.5;
301 b[0] = -1;
302 b[1] = -1;
303 c[0] = -1;
304 c[1] = -1;
305 a1d = 1;
306 b1d = -3;
307 c1d = 2;
308 aa[0] = -3;
309 aa[1] = -3;
310 bb[0][0] = 4;
311 bb[0][1] = 4;
312 bb[1][0] = 4;
313 bb[1][1] = 4;
314 out[0][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
315 out[0][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
316 out[0][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
317
318
319
320 coeff= 2;
321 a[0] = 0;
322 a[1] = -0.5;
323 b[0] = 1;
324 b[1] = 0;
325 c[0] = 1;
326 c[1] = 0;
327 a1d = 1;
328 b1d = -3;
329 c1d = 2;
330 aa[0] = -1;
331 aa[1] = 0;
332 bb[0][0] = 4;
333 bb[0][1] = 0;
334 bb[1][0] = 0;
335 bb[1][1] = 0;
336 out[1][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
337 out[1][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
338 out[1][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
339
340
341 coeff= 2;
342 a[0] = 0;
343 a[1] = -0.5;
344 b[0] = 0;
345 b[1] = 1;
346 c[0] = 0;
347 c[1] = 1;
348 a1d = 1;
349 b1d = -3;
350 c1d = 2;
351 aa[0] = 0;
352 aa[1] = -1;
353 bb[0][0] = 0;
354 bb[0][1] = 0;
355 bb[1][0] = 0;
356 bb[1][1] = 4;
357 out[2][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
358 out[2][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
359 out[2][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
360
361
362 //upper triangle
363 coeff= 2;
364 a[0] = 1;
365 a[1] = 0.5;
366 b[0] = -1;
367 b[1] = -1;
368 c[0] = -1;
369 c[1] = -1;
370 a1d = 0;
371 b1d = -1;
372 c1d = 2;
373 aa[0] = -3;
374 aa[1] = -3;
375 bb[0][0] = 4;
376 bb[0][1] = 4;
377 bb[1][0] = 4;
378 bb[1][1] = 4;
379 out[3][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
380 out[3][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
381 out[3][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
382
383
384
385 coeff= 2;
386 a[0] = 0;
387 a[1] = -0.5;
388 b[0] = 1;
389 b[1] = 0;
390 c[0] = 1;
391 c[1] = 0;
392 a1d = 0;
393 b1d = -1;
394 c1d = 2;
395 aa[0] = -1;
396 aa[1] = 0;
397 bb[0][0] = 4;
398 bb[0][1] = 0;
399 bb[1][0] = 0;
400 bb[1][1] = 0;
401 out[4][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
402 out[4][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
403 out[4][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
404
405
406
407 coeff= 2;
408 a[0] = 0;
409 a[1] = -0.5;
410 b[0] = 0;
411 b[1] = 1;
412 c[0] = 0;
413 c[1] = 1;
414 a1d = 0;
415 b1d = -1;
416 c1d = 2;
417 aa[0] = 0;
418 aa[1] = -1;
419 bb[0][0] = 0;
420 bb[0][1] = 0;
421 bb[1][0] = 0;
422 bb[1][1] = 4;
423 out[5][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
424 out[5][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
425 out[5][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
426
427
428
429 // vertical edges
430 coeff= 2;
431 a[0] = 1;
432 a[1] = 0.5;
433 b[0] = -1;
434 b[1] = -1;
435 c[0] = -1;
436 c[1] = -1;
437 a1d = 0;
438 b1d = 4;
439 c1d = -4;
440 aa[0] = -3;
441 aa[1] = -3;
442 bb[0][0] = 4;
443 bb[0][1] = 4;
444 bb[1][0] = 4;
445 bb[1][1] = 4;
446 out[6][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
447 out[6][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
448 out[6][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
449
450
451
452 coeff= 2;
453 a[0] = 0;
454 a[1] = -0.5;
455 b[0] = 1;
456 b[1] = 0;
457 c[0] = 1;
458 c[1] = 0;
459 a1d = 0;
460 b1d = 4;
461 c1d = -4;
462 aa[0] = -1;
463 aa[1] = 0;
464 bb[0][0] = 4;
465 bb[0][1] = 0;
466 bb[1][0] = 0;
467 bb[1][1] = 0;
468 out[7][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
469 out[7][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
470 out[7][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
471
472
473
474 coeff= 2;
475 a[0] = 0;
476 a[1] = -0.5;
477 b[0] = 0;
478 b[1] = 1;
479 c[0] = 0;
480 c[1] = 1;
481 a1d = 0;
482 b1d = 4;
483 c1d = -4;
484 aa[0] = 0;
485 aa[1] = -1;
486 bb[0][0] = 0;
487 bb[0][1] = 0;
488 bb[1][0] = 0;
489 bb[1][1] = 4;
490 out[8][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
491 out[8][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
492 out[8][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
493
494
495
496
497 // lower triangle edges
498 coeff= 4;
499 a[0] = 0;
500 a[1] = 1;
501 b[0] = 1;
502 b[1] = 0;
503 c[0] = -1;
504 c[1] = -1;
505 a1d = 1;
506 b1d = -3;
507 c1d = 2;
508 aa[0] = 4;
509 aa[1] = 0;
510 bb[0][0] = -8;
511 bb[0][1] = -4;
512 bb[1][0] = -4;
513 bb[1][1] = 0;
514 out[9][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
515 out[9][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
516 out[9][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
517
518
519
520
521 coeff= 4;
522 a[0] = 0;
523 a[1] = 1;
524 b[0] = 0;
525 b[1] = 1;
526 c[0] = -1;
527 c[1] = -1;
528 a1d = 1;
529 b1d = -3;
530 c1d = 2;
531 aa[0] = 0;
532 aa[1] = 4; //changed from zero to 4
533 bb[0][0] = 0;
534 bb[0][1] = -4;
535 bb[1][0] = -4;
536 bb[1][1] = -8;
537 out[10][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
538 out[10][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
539 out[10][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
540
541
542
543 coeff= 4;
544 a[0] = 0;
545 a[1] = 0;
546 b[0] = 1;
547 b[1] = 0;
548 c[0] = 0;
549 c[1] = 1;
550 a1d = 1;
551 b1d = -3;
552 c1d = 2;
553 aa[0] = 0;
554 aa[1] = 0;
555 bb[0][0] = 0;
556 bb[0][1] = 4;
557 bb[1][0] = 4;
558 bb[1][1] = 0;
559 out[11][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
560 out[11][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
561 out[11][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
562
563
564
565 // upper triangle edges
566 coeff= 4;
567 a[0] = 0;
568 a[1] = 1;
569 b[0] = 1;
570 b[1] = 0;
571 c[0] = -1;
572 c[1] = -1;
573 a1d = 0;
574 b1d = -1;
575 c1d = 2;
576 aa[0] = 4;
577 aa[1] = 0;
578 bb[0][0] = -8;
579 bb[0][1] = -4;
580 bb[1][0] = -4;
581 bb[1][1] = 0;
582 out[12][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
583 out[12][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
584 out[12][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
585
586
587
588 coeff= 4;
589 a[0] = 0;
590 a[1] = 1;
591 b[0] = 0;
592 b[1] = 1;
593 c[0] = -1;
594 c[1] = -1;
595 a1d = 0;
596 b1d = -1;
597 c1d = 2;
598 aa[0] = 0;
599 aa[1] = 4;
600 bb[0][0] = 0;
601 bb[0][1] = -4;
602 bb[1][0] = -4;
603 bb[1][1] = -8;
604 out[13][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
605 out[13][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
606 out[13][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
607
608
609
610 coeff= 4;
611 a[0] = 0;
612 a[1] = 0;
613 b[0] = 1;
614 b[1] = 0;
615 c[0] = 0;
616 c[1] = 1;
617 a1d = 0;
618 b1d = -1;
619 c1d = 2;
620 aa[0] = 0;
621 aa[1] = 0;
622 bb[0][0] = 0;
623 bb[0][1] = 4;
624 bb[1][0] = 4;
625 bb[1][1] = 0;
626 out[14][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
627 out[14][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
628 out[14][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
629
630
631
632 // quadrilateral sides
633 coeff= 4;
634 a[0] = 0;
635 a[1] = 1;
636 b[0] = 1;
637 b[1] = 0;
638 c[0] = -1;
639 c[1] = -1;
640 a1d = 0;
641 b1d = 4;
642 c1d = -4;
643 aa[0] = 4;
644 aa[1] = 0;
645 bb[0][0] = -8;
646 bb[0][1] = -4;
647 bb[1][0] = -4;
648 bb[1][1] = 0;
649 out[15][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
650 out[15][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
651 out[15][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
652
653
654
655
656
657 coeff= 4;
658 a[0] = 0;
659 a[1] = 1;
660 b[0] = 0;
661 b[1] = 1;
662 c[0] = -1;
663 c[1] = -1;
664 a1d = 0;
665 b1d = 4;
666 c1d = -4;
667 aa[0] = 0;
668 aa[1] = 4;
669 bb[0][0] = 0;
670 bb[0][1] = -4;
671 bb[1][0] = -4;
672 bb[1][1] = -8;
673 out[16][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
674 out[16][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
675 out[16][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
676
677
678
679
680 coeff= 4;
681 a[0] = 0;
682 a[1] = 0;
683 b[0] = 1;
684 b[1] = 0;
685 c[0] = 0;
686 c[1] = 1;
687 a1d = 0;
688 b1d = 4;
689 c1d = -4;
690 aa[0] = 0;
691 aa[1] = 0;
692 bb[0][0] = 0;
693 bb[0][1] = 4;
694 bb[1][0] = 4;
695 bb[1][1] = 0;
696 out[17][0][0] = (aa[0] + bb[0][0]*in[0] + bb[1][0]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
697 out[17][0][1] = (aa[1] + bb[0][1]*in[0] + bb[1][1]*in[1]) * (a1d + in[2]*b1d + in[2]*in[2]*c1d);
698 out[17][0][2] = coeff * (a[0] + b[0]*in[0] + b[1]*in[1]) * (a[1] + c[0]*in[0] + c[1]*in[1]) * (b1d + 2*c1d*in[2]);
699
700 }
701
703 void partial (const std::array<unsigned int, 3>& order,
704 const typename Traits::DomainType& in, // position
705 std::vector<typename Traits::RangeType>& out) const // return value
706 {
707 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
708 if (totalOrder == 0) {
709 evaluateFunction(in, out);
710 } else if (totalOrder == 1) {
711 out.resize(size());
712
713 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
714 switch (direction) {
715 case 0:
716 out[0] = (-3 + 4*(in[0] + in[1])) * (1 - 3*in[2] + 2*in[2]*in[2]);
717 out[1] = (-1 + 4*in[0]) * (1 - 3*in[2] + 2*in[2]*in[2]);
718 out[2] = 0;
719 out[3] = (-3 + 4*(in[0] + in[1])) * (-in[2] + 2*in[2]*in[2]);
720 out[4] = (-1 + 4*in[0]) * (-in[2] + 2*in[2]*in[2]);
721 out[5] = 0;
722 out[6] = (-3 + 4*(in[0] + in[1])) * 4*(in[2] - in[2]*in[2]);
723 out[7] = (-1 + 4*in[0]) * 4*(in[2] - in[2]*in[2]);
724 out[8] = 0;
725 out[9] = (4 - 8*in[0] - 4*in[1]) * (1 - 3*in[2] + 2*in[2]*in[2]);
726 out[10] = -4*in[1] * (1 - 3*in[2] + 2*in[2]*in[2]);
727 out[11] = 4*in[1] * (1 - 3*in[2] + 2*in[2]*in[2]);
728 out[12] = (4 - 8*in[0] - 4*in[1]) * (-in[2] + 2*in[2]*in[2]);
729 out[13] = -4*in[1] * (-in[2] + 2*in[2]*in[2]);
730 out[14] = 4*in[1] * (-in[2] + 2*in[2]*in[2]);
731 out[15] = (4 - 8*in[0] - 4*in[1]) * (4*in[2] - 4*in[2]*in[2]);
732 out[16] = -4*in[1] * (4*in[2] - 4*in[2]*in[2]);
733 out[17] = 4*in[1] * (4*in[2] - 4*in[2]*in[2]);
734 break;
735
736 case 1:
737 out[0] = (-3 + 4*(in[0] + in[1])) * (1 - 3*in[2] + 2*in[2]*in[2]);
738 out[1] = 0;
739 out[2] = (-1 + 4*in[1]) * (1 - 3*in[2] + 2*in[2]*in[2]);
740 out[3] = (-3 + 4*(in[0] + in[1])) * (-in[2] + 2*in[2]*in[2]);
741 out[4] = 0;
742 out[5] = (-1 + 4*in[1]) * (-in[2] + 2*in[2]*in[2]);
743 out[6] = (-3 + 4*(in[0] + in[1])) * 4*(in[2] - in[2]*in[2]);
744 out[7] = 0;
745 out[8] = (-1 + 4*in[1]) * 4*(in[2] - in[2]*in[2]);
746 out[9] = -4*in[0] * (1 - 3*in[2] + 2*in[2]*in[2]);
747 out[10] = (4 - 4*in[0] - 8*in[1]) * (1 - 3*in[2] + 2*in[2]*in[2]);
748 out[11] = 4*in[0] * (1 - 3*in[2] + 2*in[2]*in[2]);
749 out[12] = -4*in[0] * (-in[2] + 2*in[2]*in[2]);
750 out[13] = (4 - 4*in[0] - 8*in[1]) * (-in[2] + 2*in[2]*in[2]);
751 out[14] = 4*in[0] * (-in[2] + 2*in[2]*in[2]);
752 out[15] = -4*in[0] * 4*(in[2] - in[2]*in[2]);
753 out[16] = (4 - 4*in[0] - 8*in[1]) * (4*in[2] - 4*in[2]*in[2]);
754 out[17] = 4*in[0] * (4*in[2] - 4*in[2]*in[2]);
755 break;
756
757 case 2:
758 out[0] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]) * (-3 + 4*in[2]);
759 out[1] = 2 * in[0] * (-0.5 + in[0]) * (-3 + 4*in[2]);
760 out[2] = 2 * in[1] * (-0.5 + in[1]) * (-3 + 4*in[2]);
761 out[3] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]) * (-1 + 4*in[2]);
762 out[4] = 2 * in[0] * (-0.5 + in[0]) * (-1 + 4*in[2]);
763 out[5] = 2 * in[1] * (-0.5 + in[1]) * (-1 + 4*in[2]);
764 out[6] = 2 * (1 - in[0] - in[1]) * (0.5 - in[0] - in[1]) * (4 - 8*in[2]);
765 out[7] = 2 * in[0] * (-0.5 + in[0]) * (4 - 8*in[2]);
766 out[8] = 2*in[1] * (-0.5 + in[1]) * (4 - 8*in[2]);
767 out[9] = 4*in[0] * (1 - in[0] - in[1]) * (-3 + 4*in[2]);
768 out[10] = 4*in[1] * (1 - in[0] - in[1]) * (-3 + 4*in[2]);
769 out[11] = 4*in[0]*in[1] * (-3 + 4*in[2]);
770 out[12] = 4 * in[0] * (1 - in[0] - in[1]) * (-1 + 4*in[2]);
771 out[13] = 4 * in[1] * (1 - in[0] - in[1]) * (-1 + 4*in[2]);
772 out[14] = 4*in[0]*in[1] * (-1 + 4*in[2]);
773 out[15] = 4 * in[0] * (1 - in[0] - in[1]) * (4 - 8*in[2]);
774 out[16] = 4 * in[1] * (1 - in[0] - in[1]) * (4 - 8*in[2]);
775 out[17] = 4*in[0]*in[1] * (4 - 8*in[2]);
776 break;
777
778 default:
779 DUNE_THROW(RangeError, "Component out of range.");
780 }
781 } else {
782 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
783 }
784 }
785
787 unsigned int order () const
788 {
789 return 2;
790 }
791 };
792}
793#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::PrismP2LocalBasis
Quadratic Lagrange shape functions on the prism.
Definition: prismp2localbasis.hh:26
Dune::PrismP2LocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: prismp2localbasis.hh:286
Dune::PrismP2LocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: prismp2localbasis.hh:39
Dune::PrismP2LocalBasis::size
unsigned int size() const
number of shape functions
Definition: prismp2localbasis.hh:33
Dune::PrismP2LocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: prismp2localbasis.hh:787
Dune::PrismP2LocalBasis::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 all shape functions.
Definition: prismp2localbasis.hh:703
Dune::PrismP2LocalBasis::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: prismp2localbasis.hh:30
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)