Dune Core Modules (2.6.0)

pyramidp2localbasis.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_PYRAMID_P2_LOCALBASIS_HH
4#define DUNE_PYRAMID_P2_LOCALBASIS_HH
5
6#include <numeric>
7
8#include <dune/common/fmatrix.hh>
9
10#include <dune/localfunctions/common/localbasis.hh>
11
12namespace Dune
13{
28 template<class D, class R>
29 class PyramidP2LocalBasis
30 {
31 public:
33 typedef LocalBasisTraits<D,3,Dune::FieldVector<D,3>,R,1,Dune::FieldVector<R,1>,
34 Dune::FieldMatrix<R,1,3> > Traits;
35
37 unsigned int size () const
38 {
39 return 14;
40 }
41
43 inline void evaluateFunction (const typename Traits::DomainType& in,
44 std::vector<typename Traits::RangeType>& out) const
45 {
46 out.resize(14);
47
48 // transform to reference element with base [-1,1]^2
49 const R x = 2.0*in[0] + in[2] - 1.0;
50 const R y = 2.0*in[1] + in[2] - 1.0;
51 const R z = in[2];
52
53 if (x > y)
54 {
55 // vertices
56 out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
57 out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
58 out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
59 out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
60 out[4] = z*(2*z - 1);
61
62 // lower edges
63 out[5] = -0.5*(y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
64 out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
65 out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
66 out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
67
68 // upper edges
69 out[9] = z*(x + z - 1)*(y - z - 1);
70 out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
71 out[11] = -z*(y - z + 1)*(x + z - 1);
72 out[12] = z*(y - z + 1)*(x + z + 1);
73
74 // base face
75 out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
76 }
77 else
78 {
79 // vertices
80 out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
81 out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
82 out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
83 out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
84 out[4] = z*(2*z - 1);
85
86 // lower edges
87 out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
88 out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
89 out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
90 out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
91
92 // upper edges
93 out[9] = z*(y + z - 1)*(x - z - 1);
94 out[10] = -z*(x - z + 1)*(y + z - 1);
95 out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
96 out[12] = z*(x - z + 1)*(y + z + 1);
97
98 // base face
99 out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
100 }
101 }
102
104 inline void
105 evaluateJacobian (const typename Traits::DomainType& in, // position
106 std::vector<typename Traits::JacobianType>& out) const // return value
107 {
108 out.resize(14);
109
110 // transform to reference element with base [-1,1]^2
111 const R x = 2.0*in[0] + in[2] - 1.0;
112 const R y = 2.0*in[1] + in[2] - 1.0;
113 const R z = in[2];
114
115 // transformation of the gradient leads to a multiplication
116 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
117 if (x > y)
118 {
119 // vertices
120 out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
121 out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
122 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
123 + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
124 + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
125
126 out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
127 + (x + z)*(y - z)*(-y + z + 1)) - z);
128 out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
129 + (x + z)*(y - z)*(-(x + z + 1))) + z);
130 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
131 - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
132 - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
133 + (x + z)*(y - z)*(x - y + 2*z - 2))
134 - (x - y);
135
136 out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
137 out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
138 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
139 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
140 + (x + z)*(y - z)*(y - x - 2*z + 2));
141
142 out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
143 out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
144 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
145 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
146 + (y - z)*(x + z)*(y - x - 2*z));
147
148 out[4][0][0] = 0;
149 out[4][0][1] = 0;
150 out[4][0][2] = 4*z - 1;
151
152 // lower edges
153 out[5][0][0] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
154 + (y - z + 1)*(x + z - 1)*((y - 1) + z));
155 out[5][0][1] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
156 + (y - z + 1)*(x + z - 1)*((x + 1) - z));
157 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
158 - 0.5*((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
159 + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
160
161 out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
162 out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
163 + (y - z + 1)*((x + z + 1)*x + 2*z));
164 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
165 - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
166 + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
167
168 out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
169 + (x + z - 1)*((y - z - 1)*y + 2*z));
170 out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
171 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
172 - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
173 + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
174
175 out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
176 out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
177 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
178 - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
179
180 // upper edges
181 out[9][0][0] = 2*z*(y - z - 1);
182 out[9][0][1] = 2*z*(x + z - 1);
183 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
184 + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
185
186 out[10][0][0] = -2*z*(y - z - 1);
187 out[10][0][1] = -2*z*(x + z + 1);
188 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
189 - ((x + z + 1)*(y - z - 1) + 4*z)
190 - z*(-x + y - 2*z + 2);
191
192 out[11][0][0] = -2*z*(y - z + 1);
193 out[11][0][1] = -2*z*(x + z - 1);
194 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
195 - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
196
197 out[12][0][0] = 2*z*(y - z + 1);
198 out[12][0][1] = 2*z*(x + z + 1);
199 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
200 + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
201
202 // base face
203 out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
204 + (y - z + 1)*(x + z - 1)*(y - 1 + z));
205 out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
206 + (y - z + 1)*(x + z - 1)*(x + 1 - z));
207 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
208 + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
209 + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
210 }
211 else
212 {
213 // vertices
214 out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
215 out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
216 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
217 + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
218 + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
219
220 out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
221 out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
222 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
223 - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
224 + (x - z)*(y + z)*(-x + y + 2*z - 2));
225
226 out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
227 + (x - z)*(y + z)*(y + z + 1) + 4*z);
228 out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
229 + (x - z)*(y + z)*(x - z - 1) - 4*z);
230 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
231 + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
232 + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
233
234 out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
235 out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
236 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
237 + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
238 + (y + z)*(x - z)*(x - y - 2*z));
239
240 out[4][0][0] = 0;
241 out[4][0][1] = 0;
242 out[4][0][2] = 4*z - 1;
243
244 // lower edges
245 out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
246 out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
247 + (y + z - 1)*((x - z - 1)*x + 2*z));
248 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
249 - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
250 + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
251
252 out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
253 out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
254 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
255 - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
256
257 out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
258 out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
259 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
260 - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
261
262 out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
263 + (x - z + 1)*((y + z + 1)*y + 2*z));
264 out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
265 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
266 - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
267 + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
268
269 // upper edges
270 out[9][0][0] = 2*z*(y + z - 1);
271 out[9][0][1] = 2*z*(x - z - 1);
272 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
273 + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
274
275 out[10][0][0] = -2*z*(y + z - 1);
276 out[10][0][1] = -2*z*(x - z + 1);
277 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
278 - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
279
280 out[11][0][0] = -2*z*(y + z + 1);
281 out[11][0][1] = -2*z*(x - z - 1);
282 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
283 - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
284
285 out[12][0][0] = 2*z*(y + z + 1);
286 out[12][0][1] = 2*z*(x - z + 1);
287 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
288 + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
289
290 // base face
291 out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
292 + (x - z + 1)*(y + z - 1)*(y + 1 - z));
293 out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
294 + (x - z + 1)*(y + z - 1)*(x - 1 + z));
295 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
296 + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
297 + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
298 }
299 }
300
302 void partial (const std::array<unsigned int, 3>& order,
303 const typename Traits::DomainType& in, // position
304 std::vector<typename Traits::RangeType>& out) const // return value
305 {
306 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
307 if (totalOrder == 0) {
308 evaluateFunction(in, out);
309 } else if (totalOrder == 1) {
310 out.resize(size());
311
312 // transform to reference element with base [-1,1]^2
313 const R x = 2.0*in[0] + in[2] - 1.0;
314 const R y = 2.0*in[1] + in[2] - 1.0;
315 const R z = in[2];
316
317 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
318
319 // transformation of the gradient leads to a multiplication
320 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
321 if (x > y) {
322 switch (direction) {
323 case 0:
324 out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
325 out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
326 out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
327 out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
328 out[4] = 0;
329 out[5] = -((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((y - 1) + z));
330 out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
331 out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
332 out[8] = -(y - z + 1)*(2*x + z)*y;
333 out[9] = 2*z*(y - z - 1);
334 out[10] = -2*z*(y - z - 1);
335 out[11] = -2*z*(y - z + 1);
336 out[12] = 2*z*(y - z + 1);
337 out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
338 break;
339 case 1:
340 out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
341 out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
342 out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
343 out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
344 out[4] = 0;
345 out[5] = -((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*((x + 1) - z));
346 out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
347 out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
348 out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
349 out[9] = 2*z*(x + z - 1);
350 out[10] = -2*z*(x + z + 1);
351 out[11] = -2*z*(x + z - 1);
352 out[12] = 2*z*(x + z + 1);
353 out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
354 break;
355 case 2:
356 out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
357 out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
358 out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
359 out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
360 out[4] = 4*z - 1;
361 out[5] = -((y - z + 1)*(2*y - 3*x + z + 2*x*y + 6*x*z - 2*y*z + 2*x*x + 4*z*z - 3))/2;
362 out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
363 out[7] = z - z*(2*x + 1) - ((2*z - y*(z - y + 1))*(x + z - 1))/2 - ((2*x - y*(x + 1))*(x + z - 1))/2 + ((x + 1)*(x + z - 1)*(z - 2*y + 1))/2 + y*(x + 1)*(z - y + 1);
364 out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
365 out[9] = -(x + 3*z - 1)*(z - y + 1);
366 out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
367 out[11] = -(x + 3*z - 1)*(y - z + 1);
368 out[12] = (x + 3*z + 1)*(y - z + 1);
369 out[13] = (y - z + 1)*(2*y - 3*x + z + 2*x*y + 6*x*z - 2*y*z + 2*x*x + 4*z*z - 3);
370 break;
371 default:
372 DUNE_THROW(RangeError, "Component out of range.");
373 }
374 }
375 else // x <= y
376 {
377 switch (direction) {
378 case 0:
379 out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
380 out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
381 out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
382 out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
383 out[4] = 0;
384 out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
385 out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
386 out[7] = -y*(2*x - z)*(y + z - 1);
387 out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
388 out[9] = 2*z*(y + z - 1);
389 out[10] = -2*z*(y + z - 1);
390 out[11] = -2*z*(y + z + 1);
391 out[12] = 2*z*(y + z + 1);
392 out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
393 break;
394 case 1:
395 out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
396 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
397 out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
398 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
399 out[4] = 0;
400 out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
401 out[6] = -x*(2*y + z)*(x - z + 1);
402 out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
403 out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
404 out[9] = -2*z*(z - x + 1);
405 out[10] = -2*z*(x - z + 1);
406 out[11] = 2*z*(z - x + 1);
407 out[12] = 2*z*(x - z + 1);
408 out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
409 break;
410 case 2:
411 out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
412 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
413 out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
414 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
415 out[4] = 4*z - 1;
416 out[5] = z - z*(2*y + 1) - ((2*z - x*(z - x + 1))*(y + z - 1))/2 - ((2*y - x*(y + 1))*(y + z - 1))/2 + ((y + 1)*(y + z - 1)*(z - 2*x + 1))/2 + x*(y + 1)*(z - x + 1);
417 out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
418 out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
419 out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
420 out[9] = -(z - x + 1)*(y + 3*z - 1);
421 out[10] = -(x - z + 1)*(y + 3*z - 1);
422 out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
423 out[12] = (x - z + 1)*(y + 3*z + 1);
424 out[13] = (x - z + 1)*(2*x - 3*y + z + 2*x*y - 2*x*z + 6*y*z + 2*y*y + 4*z*z - 3);
425 break;
426 default:
427 DUNE_THROW(RangeError, "Component out of range.");
428 }
429 }
430 } else {
431 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
432 }
433 }
434
436 unsigned int order () const
437 {
438 return 2;
439 }
440 };
441}
442#endif
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::PyramidP2LocalBasis
Quadratic Lagrange shape functions on the pyramid.
Definition: pyramidp2localbasis.hh:30
Dune::PyramidP2LocalBasis::partial
void partial(const std::array< unsigned int, 3 > &order, const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate partial derivatives of all shape functions.
Definition: pyramidp2localbasis.hh:302
Dune::PyramidP2LocalBasis::order
unsigned int order() const
Polynomial order of the shape functions.
Definition: pyramidp2localbasis.hh:436
Dune::PyramidP2LocalBasis::size
unsigned int size() const
number of shape functions
Definition: pyramidp2localbasis.hh:37
Dune::PyramidP2LocalBasis::Traits
LocalBasisTraits< D, 3, Dune::FieldVector< D, 3 >, R, 1, Dune::FieldVector< R, 1 >, Dune::FieldMatrix< R, 1, 3 > > Traits
export type traits for function signature
Definition: pyramidp2localbasis.hh:34
Dune::PyramidP2LocalBasis::evaluateFunction
void evaluateFunction(const typename Traits::DomainType &in, std::vector< typename Traits::RangeType > &out) const
Evaluate all shape functions.
Definition: pyramidp2localbasis.hh:43
Dune::PyramidP2LocalBasis::evaluateJacobian
void evaluateJacobian(const typename Traits::DomainType &in, std::vector< typename Traits::JacobianType > &out) const
Evaluate Jacobian of all shape functions.
Definition: pyramidp2localbasis.hh:105
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 (Nov 24, 23:30, 2024)