DUNE PDELab (2.7)

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