dune-localfunctions  2.2.1
p23dlocalbasis.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil -*-
2 #ifndef DUNE_P2_3DLOCALBASIS_HH
3 #define DUNE_P2_3DLOCALBASIS_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 10;
33  }
34 
36  inline void evaluateFunction (const typename Traits::DomainType& in,
37  std::vector<typename Traits::RangeType>& out) const
38  {
39  out.resize(10);
40 
41  int coeff;
42  R a[2], b[3], c[3];
43 
44  // case 0:
45  coeff=2;
46  a[0]=1.0;
47  a[1]=0.5;
48  b[0]=-1.0;
49  b[1]=-1.0;
50  b[2]=-1.0;
51  c[0]=-1.0;
52  c[1]=-1.0;
53  c[2]=-1.0;
54 
55  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]);
56 
57  // case 1:
58  coeff=2;
59  a[0]=0.0;
60  a[1]=-0.5;
61  b[0]=1.0;
62  b[1]=0.0;
63  b[2]=0.0;
64  c[0]=1.0;
65  c[1]=0.0;
66  c[2]=0.0;
67 
68  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]);
69 
70  // case 2:
71  coeff=2;
72  a[0]=0.0;
73  a[1]=-0.5;
74  b[0]=0.0;
75  b[1]=1.0;
76  b[2]=0.0;
77  c[0]=0.0;
78  c[1]=1.0;
79  c[2]=0.0;
80 
81  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]);
82 
83  // case 3:
84  coeff=2;
85  a[0]=0.0;
86  a[1]=-0.5;
87  b[0]=0.0;
88  b[1]=0.0;
89  b[2]=1.0;
90  c[0]=0.0;
91  c[1]=0.0;
92  c[2]=1.0;
93 
94  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]);
95 
96  // case 4:
97  coeff=4;
98  a[0]=0.0;
99  a[1]=1.0;
100  b[0]=1.0;
101  b[1]=0.0;
102  b[2]=0.0;
103  c[0]=-1.0;
104  c[1]=-1.0;
105  c[2]=-1.0;
106 
107  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]);
108 
109  // case 5:
110  coeff=4;
111  a[0]=0.0;
112  a[1]=0.0;
113  b[0]=1.0;
114  b[1]=0.0;
115  b[2]=0.0;
116  c[0]=0.0;
117  c[1]=1.0;
118  c[2]=0.0;
119 
120  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]);
121 
122  // case 6:
123  coeff=4;
124  a[0]=0.0;
125  a[1]=1.0;
126  b[0]=0.0;
127  b[1]=1.0;
128  b[2]=0.0;
129  c[0]=-1.0;
130  c[1]=-1.0;
131  c[2]=-1.0;
132 
133  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]);
134 
135  // case 7:
136  coeff=4;
137  a[0]=0.0;
138  a[1]=1.0;
139  b[0]=0.0;
140  b[1]=0.0;
141  b[2]=1.0;
142  c[0]=-1.0;
143  c[1]=-1.0;
144  c[2]=-1.0;
145 
146  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]);
147 
148  // case 8:
149  coeff=4;
150  a[0]=0.0;
151  a[1]=0.0;
152  b[0]=1.0;
153  b[1]=0.0;
154  b[2]=0.0;
155  c[0]=0.0;
156  c[1]=0.0;
157  c[2]=1.0;
158 
159  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]);
160 
161  // case 9:
162  coeff=4;
163  a[0]=0.0;
164  a[1]=0.0;
165  b[0]=0.0;
166  b[1]=1.0;
167  b[2]=0.0;
168  c[0]=0.0;
169  c[1]=0.0;
170  c[2]=1.0;
171 
172  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]);
173 
174  }
175 
177  inline void
178  evaluateJacobian (const typename Traits::DomainType& in, // position
179  std::vector<typename Traits::JacobianType>& out) const // return value
180  {
181  out.resize(10);
182 
183  R aa[3][3], bb[3][3];
184  // case 0:
185  //x derivative
186  aa[0][0]=-3.0;
187  bb[0][0]=4.0;
188  bb[1][0]=4.0;
189  bb[2][0]=4.0;
190  //y derivative
191  aa[0][1]=-3.0;
192  bb[0][1]=4.0;
193  bb[1][1]=4.0;
194  bb[2][1]=4.0;
195  // z derivative
196  aa[0][2]=-3.0;
197  bb[0][2]=4.0;
198  bb[1][2]=4.0;
199  bb[2][2]=4.0;
200 
201  out[0][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
202  out[0][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
203  out[0][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
204 
205  //case 1:
206  //x derivative
207  aa[0][0]=-1.0;
208  bb[0][0]=4.0;
209  bb[1][0]=0.0;
210  bb[2][0]=0.0;
211  //y derivative
212  aa[0][1]=0.0;
213  bb[0][1]=0.0;
214  bb[1][1]=0.0;
215  bb[2][1]=0.0;
216  // z derivative
217  aa[0][2]=0.0;
218  bb[0][2]=0.0;
219  bb[1][2]=0.0;
220  bb[2][2]=0.0;
221 
222  out[1][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
223  out[1][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
224  out[1][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
225 
226  // case 2:
227  //x derivative
228  aa[0][0]=0.0;
229  bb[0][0]=0.0;
230  bb[1][0]=0.0;
231  bb[2][0]=0.0;
232  //y derivative
233  aa[0][1]=-1.0;
234  bb[0][1]=0.0;
235  bb[1][1]=4.0;
236  bb[2][1]=0.0;
237  // z derivative
238  aa[0][2]=0.0;
239  bb[0][2]=0.0;
240  bb[1][2]=0.0;
241  bb[2][2]=0.0;
242 
243  out[2][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
244  out[2][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
245  out[2][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
246 
247  // case 3:
248  //x derivative
249  aa[0][0]=0.0;
250  bb[0][0]=0.0;
251  bb[1][0]=0.0;
252  bb[2][0]=0.0;
253  //y derivative
254  aa[0][1]=0.0;
255  bb[0][1]=0.0;
256  bb[1][1]=0.0;
257  bb[2][1]=0.0;
258  // z derivative
259  aa[0][2]=-1.0;
260  bb[0][2]=0.0;
261  bb[1][2]=0.0;
262  bb[2][2]=4.0;
263 
264  out[3][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
265  out[3][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
266  out[3][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
267 
268  // case 4:
269  //x derivative
270  aa[0][0]=4.0;
271  bb[0][0]=-8.0;
272  bb[1][0]=-4.0;
273  bb[2][0]=-4.0;
274  //y derivative
275  aa[0][1]=0.0;
276  bb[0][1]=-4.0;
277  bb[1][1]=0.0;
278  bb[2][1]=0.0;
279  // z derivative
280  aa[0][2]=0.0;
281  bb[0][2]=-4.0;
282  bb[1][2]=0.0;
283  bb[2][2]=0.0;
284 
285  out[4][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
286  out[4][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
287  out[4][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
288 
289  // case 5:
290  //x derivative
291  aa[0][0]=0.0;
292  bb[0][0]=0.0;
293  bb[1][0]=4.0;
294  bb[2][0]=0.0;
295  //y derivative
296  aa[0][1]=0.0;
297  bb[0][1]=4.0;
298  bb[1][1]=0.0;
299  bb[2][1]=0.0;
300  // z derivative
301  aa[0][2]=0.0;
302  bb[0][2]=0.0;
303  bb[1][2]=0.0;
304  bb[2][2]=0.0;
305 
306  out[5][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
307  out[5][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
308  out[5][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
309 
310  // case 6:
311  //x derivative
312  aa[0][0]=0.0;
313  bb[0][0]=0.0;
314  bb[1][0]=-4.0;
315  bb[2][0]=0.0;
316  //y derivative
317  aa[0][1]=4.0;
318  bb[0][1]=-4.0;
319  bb[1][1]=-8.0;
320  bb[2][1]=-4.0;
321  // z derivative
322  aa[0][2]=0.0;
323  bb[0][2]=0.0;
324  bb[1][2]=-4.0;
325  bb[2][2]=0.0;
326 
327  out[6][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
328  out[6][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
329  out[6][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
330 
331  // case 7:
332  //x derivative
333  aa[0][0]=0.0;
334  bb[0][0]=0.0;
335  bb[1][0]=0.0;
336  bb[2][0]=-4.0;
337  //y derivative
338  aa[0][1]=0.0;
339  bb[0][1]=0.0;
340  bb[1][1]=0.0;
341  bb[2][1]=-4.0;
342  // z derivative
343  aa[0][2]=4.0;
344  bb[0][2]=-4.0;
345  bb[1][2]=-4.0;
346  bb[2][2]=-8.0;
347 
348  out[7][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
349  out[7][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
350  out[7][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
351 
352  //case 8:
353  //x derivative
354  aa[0][0]=0.0;
355  bb[0][0]=0.0;
356  bb[1][0]=0.0;
357  bb[2][0]=4.0;
358  //y derivative
359  aa[0][1]=0.0;
360  bb[0][1]=0.0;
361  bb[1][1]=0.0;
362  bb[2][1]=0.0;
363  // z derivative
364  aa[0][2]=0.0;
365  bb[0][2]=4.0;
366  bb[1][2]=0.0;
367  bb[2][2]=0.0;
368 
369  out[8][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
370  out[8][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
371  out[8][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
372 
373  // case 9:
374  //x derivative
375  aa[0][0]=0.0;
376  bb[0][0]=0.0;
377  bb[1][0]=0.0;
378  bb[2][0]=0.0;
379  //y derivative
380  aa[0][1]=0.0;
381  bb[0][1]=0.0;
382  bb[1][1]=0.0;
383  bb[2][1]=4.0;
384  // z derivative
385  aa[0][2]=0.0;
386  bb[0][2]=0.0;
387  bb[1][2]=4.0;
388  bb[2][2]=0.0;
389 
390  out[9][0][0] = aa[0][0] + bb[0][0]*in[0] + bb[1][0]*in[1] + bb[2][0]*in[2];
391  out[9][0][1] = aa[0][1] + bb[0][1]*in[0] + bb[1][1]*in[1] + bb[2][1]*in[2];
392  out[9][0][2] = aa[0][2] + bb[0][2]*in[0] + bb[1][2]*in[1] + bb[2][2]*in[2];
393 
394  }
395 
397  unsigned int order () const
398  {
399  return 2;
400  }
401  };
402 }
403 #endif