Dune Core Modules (2.7.0)

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