Dune Core Modules (2.6.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
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
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: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::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: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.111.3 (Jul 15, 22:36, 2024)