DUNE PDELab (git)

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