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