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 
12 namespace 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.80.0 (May 4, 22:30, 2024)