\n+
\n+
\n+
\n+
\n+
\n
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
65 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[1]);
\n-
66 out[1] = in[0]*(1-in[1])-in[2]*in[1];
\n-
67 out[2] = (1-in[0])*in[1]-in[2]*in[1];
\n-
68 out[3] = in[0]*in[1]+in[2]*in[1];
\n-
\n-
\n-
\n-
72 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[0]);
\n-
73 out[1] = in[0]*(1-in[1])-in[2]*in[0];
\n-
74 out[2] = (1-in[0])*in[1]-in[2]*in[0];
\n-
75 out[3] = in[0]*in[1]+in[2]*in[0];
\n-
\n+
\n+
\n+
\n+
\n+
62 std::vector<typename Traits::RangeType>& out)
const \n+
\n+
\n+
\n+
\n+
67 out[1] = 1-4*(in[0]-0.5)*(in[0]-0.5);
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
74 std::vector<typename Traits::JacobianType>& out)
const \n+
\n+
\n
\n-
\n-
\n-
\n-
\n+
\n+
79 out[1][0][0] = 4-8*in[0];
\n+
\n+
\n+
\n
\n-
\n-
\n-
\n-
86 const R x = 2.0*in[0] + in[2] - 1.0;
\n-
87 const R y = 2.0*in[1] + in[2] - 1.0;
\n-
\n-
\n-
\n-
\n-
\n-
93 out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
\n-
94 out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
\n-
95 out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
\n-
96 out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
\n-
\n-
\n-
\n-
100 out[5] = -0.5*(y - z + 1)*(x + z - 1)*(y - 1)*x;
\n-
101 out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
\n-
102 out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
\n-
103 out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
\n-
\n-
\n-
106 out[9] = z*(x + z - 1)*(y - z - 1);
\n-
107 out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
\n-
108 out[11] = -z*(y - z + 1)*(x + z - 1);
\n-
109 out[12] = z*(y - z + 1)*(x + z + 1);
\n-
\n-
\n-
112 out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
\n-
\n-
\n-
\n-
\n-
117 out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
\n-
118 out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
\n-
119 out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
\n-
120 out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
\n-
121 out[4] = z*(2*z - 1);
\n-
\n-
\n-
124 out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
\n-
125 out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
\n-
126 out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
\n-
127 out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
\n-
\n-
\n-
130 out[9] = z*(y + z - 1)*(x - z - 1);
\n-
131 out[10] = -z*(x - z + 1)*(y + z - 1);
\n-
132 out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
\n-
133 out[12] = z*(x - z + 1)*(y + z + 1);
\n-
\n-
\n-
136 out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
\n-
\n-
\n-
\n-
\n-
\n-
142 DUNE_THROW(NotImplemented,
"LagrangePyramidLocalBasis::evaluateFunction for order " << k);
\n-
\n-
\n-
\n-
151 std::vector<typename Traits::JacobianType>& out)
const \n-
\n-
\n-
\n-
\n-
\n-
\n-
158 std::fill(out[0][0].begin(), out[0][0].end(), 0);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
166 out[0][0] = {-1 + in[1], -1 + in[0] + in[2], -1 + in[1]};
\n-
167 out[1][0] = { 1 - in[1], -in[0] - in[2], -in[1]};
\n-
168 out[2][0] = { -in[1], 1 - in[0] - in[2], -in[1]};
\n-
169 out[3][0] = { in[1], in[0] + in[2], in[1]};
\n-
\n-
\n-
\n-
173 out[0][0] = {-1 + in[1] + in[2], -1 + in[0], -1 + in[0]};
\n-
174 out[1][0] = { 1 - in[1] - in[2], -in[0], -in[0]};
\n-
175 out[2][0] = { -in[1] - in[2], 1 - in[0], -in[0]};
\n-
176 out[3][0] = { in[1] + in[2], in[0], in[0]};
\n-
\n-
\n-
179 out[4][0] = {0, 0, 1};
\n-
\n-
\n-
\n-
\n-
\n-
\n-
186 const R x = 2.0*in[0] + in[2] - 1.0;
\n-
187 const R y = 2.0*in[1] + in[2] - 1.0;
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
195 out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
\n-
196 out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
\n-
197 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
\n-
198 + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
\n-
199 + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
\n-
\n-
201 out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
\n-
202 + (x + z)*(y - z)*(-y + z + 1)) - z);
\n-
203 out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
\n-
204 + (x + z)*(y - z)*(-(x + z + 1))) + z);
\n-
205 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
\n-
206 - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
\n-
207 - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
\n-
208 + (x + z)*(y - z)*(x - y + 2*z - 2))
\n-
\n-
\n-
211 out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
\n-
212 out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
\n-
213 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
\n-
214 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
\n-
215 + (x + z)*(y - z)*(y - x - 2*z + 2));
\n+
\n+
84 void partial (
const std::array<unsigned int, 1>& order,
\n+
\n+
86 std::vector<typename Traits::RangeType>& out)
const \n+
\n+
88 auto totalOrder = order[0];
\n+
89 if (totalOrder == 0) {
\n+
90 evaluateFunction(in, out);
\n+
91 }
else if (totalOrder == 1) {
\n+
\n+
\n+
\n+
\n+
96 }
else if (totalOrder == 2) {
\n+
\n+
\n+
\n+
\n+
\n+
\n+
103 out[0] = out[1] = out[2] = 0;
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n
\n-
258 out[10][0][0] = -2*z*(y - z - 1);
\n-
259 out[10][0][1] = -2*z*(x + z + 1);
\n-
260 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
\n-
261 - ((x + z + 1)*(y - z - 1) + 4*z)
\n-
262 - z*(-x + y - 2*z + 2);
\n+
\n
\n-
264 out[11][0][0] = -2*z*(y - z + 1);
\n-
265 out[11][0][1] = -2*z*(x + z - 1);
\n-
266 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
\n-
267 - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
\n-
\n-
269 out[12][0][0] = 2*z*(y - z + 1);
\n-
270 out[12][0][1] = 2*z*(x + z + 1);
\n-
271 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
\n-
272 + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
\n-
\n-
\n-
275 out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
\n-
276 + (y - z + 1)*(x + z - 1)*(y - 1 + z));
\n-
277 out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
\n-
278 + (y - z + 1)*(x + z - 1)*(x + 1 - z));
\n-
279 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
\n-
280 + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
\n-
281 + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
\n-
\n-
\n-
\n-
\n-
286 out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
\n-
287 out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
\n-
288 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
\n-
289 + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
\n-
290 + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
\n-
\n-
292 out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
\n-
293 out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
\n-
294 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
\n-
295 - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
\n-
296 + (x - z)*(y + z)*(-x + y + 2*z - 2));
\n-
\n-
298 out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
\n-
299 + (x - z)*(y + z)*(y + z + 1) + 4*z);
\n-
300 out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
\n-
301 + (x - z)*(y + z)*(x - z - 1) - 4*z);
\n-
302 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
\n-
303 + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
\n-
304 + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
\n-
\n-
306 out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
\n-
307 out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
\n-
308 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
\n-
309 + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
\n-
310 + (y + z)*(x - z)*(x - y - 2*z));
\n+
\n+
\n+
266 std::vector<typename Traits::RangeType>& out)
const \n+
\n+
\n+
\n+
270 out[0] = 1 - in[0] - in[1] - in[2];
\n+
271 out[1] = 4 * in[0] * (1 - in[0] - in[1] - in[2]);
\n+
\n+
273 out[3] = 4 * in[1] * (1 - in[0] - in[1] - in[2]);
\n+
274 out[4] = 4 * in[0] * in[1];
\n+
\n+
276 out[6] = 4 * in[2] * (1 - in[0] - in[1] - in[2]);
\n+
277 out[7] = 4 * in[0] * in[2];
\n+
278 out[8] = 4 * in[1] * in[2];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
284 std::vector<typename Traits::JacobianType>& out)
const \n+
\n+
\n+
\n+
288 out[0][0][0] = -1; out[0][0][1] = -1; out[0][0][2] = -1;
\n+
289 out[1][0][0] = 4-8*in[0]-4*in[1]-4*in[2]; out[1][0][1] = -4*in[0]; out[1][0][2] = -4*in[0];
\n+
290 out[2][0][0] = 1; out[2][0][1] = 0; out[2][0][2] = 0;
\n+
291 out[3][0][0] = -4*in[1]; out[3][0][1] = 4-4*in[0]-8*in[1]-4*in[2]; out[3][0][2] = -4*in[1];
\n+
292 out[4][0][0] = 4*in[1]; out[4][0][1] = 4*in[0]; out[4][0][2] = 0;
\n+
293 out[5][0][0] = 0; out[5][0][1] = 1; out[5][0][2] = 0;
\n+
294 out[6][0][0] = -4*in[2]; out[6][0][1] = -4*in[2]; out[6][0][2] = 4-4*in[0]-4*in[1]-8*in[2];
\n+
295 out[7][0][0] = 4*in[2]; out[7][0][1] = 0; out[7][0][2] = 4*in[0];
\n+
296 out[8][0][0] = 0; out[8][0][1] = 4*in[2]; out[8][0][2] = 4*in[1];
\n+
297 out[9][0][0] = 0; out[9][0][1] = 0; out[9][0][2] = 1;
\n+
\n+
\n+
\n+
\n+
301 void partial (
const std::array<unsigned int, 3>& order,
\n+
\n+
303 std::vector<typename Traits::RangeType>& out)
const \n+
\n+
305 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
\n+
306 if (totalOrder == 0) {
\n+
307 evaluateFunction(in, out);
\n+
308 }
else if (totalOrder == 1) {
\n+
\n+
310 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
\n
\n-
\n-
\n-
314 out[4][0][2] = 4*z - 1;
\n-
\n-
\n-
317 out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
\n-
318 out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
\n-
319 + (y + z - 1)*((x - z - 1)*x + 2*z));
\n-
320 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
\n-
321 - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
\n-
322 + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
\n-
\n-
324 out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
\n-
325 out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
\n-
326 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
\n-
327 - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
\n-
\n-
329 out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
\n-
330 out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
\n-
331 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
\n-
332 - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
\n-
\n-
334 out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
\n-
335 + (x - z + 1)*((y + z + 1)*y + 2*z));
\n-
336 out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
\n-
337 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
\n-
338 - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
\n-
339 + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
\n-
\n-
\n-
342 out[9][0][0] = 2*z*(y + z - 1);
\n-
343 out[9][0][1] = 2*z*(x - z - 1);
\n-
344 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
\n-
345 + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
\n-
\n-
347 out[10][0][0] = -2*z*(y + z - 1);
\n-
348 out[10][0][1] = -2*z*(x - z + 1);
\n-
349 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
\n-
350 - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
\n-
\n-
352 out[11][0][0] = -2*z*(y + z + 1);
\n-
353 out[11][0][1] = -2*z*(x - z - 1);
\n-
354 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
\n-
355 - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
\n+
\n+
\n+
\n+
315 out[1] = 4-8*in[0]-4*in[1]-4*in[2];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
329 out[3] = 4-4*in[0]-8*in[1]-4*in[2];
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
\n+
344 out[6] = 4-4*in[0]-4*in[1]-8*in[2];
\n+
\n+
\n+
\n+
\n+
\n+
350 DUNE_THROW(RangeError,
"Component out of range.");
\n+
\n+
\n+
353 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
\n+
\n+
\n+
\n
\n-
357 out[12][0][0] = 2*z*(y + z + 1);
\n-
358 out[12][0][1] = 2*z*(x - z + 1);
\n-
359 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
\n-
360 + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
\n-
\n-
\n-
363 out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
\n-
364 + (x - z + 1)*(y + z - 1)*(y + 1 - z));
\n-
365 out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
\n-
366 + (x - z + 1)*(y + z - 1)*(x - 1 + z));
\n-
367 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
\n-
368 + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
\n-
369 + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
\n-
\n-
\n-
\n-
\n-
\n-
375 DUNE_THROW(NotImplemented,
"LagrangePyramidLocalBasis::evaluateJacobian for order " << k);
\n-
\n-
\n-
384 void partial(
const std::array<unsigned int,3>& order,
\n-
\n-
386 std::vector<typename Traits::RangeType>& out)
const \n-
\n-
388 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
\n-
\n-
\n-
\n-
\n-
\n-
394 evaluateFunction(in, out);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
408 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
421 out[0] = -1 + in[0] + in[2];
\n-
422 out[1] = -in[0] - in[2];
\n-
423 out[2] = 1 - in[0] - in[2];
\n-
424 out[3] = in[0]+in[2];
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
435 DUNE_THROW(RangeError,
"Component out of range.");
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
443 out[0] = -1 + in[1] + in[2];
\n-
444 out[1] = 1 - in[1] - in[2];
\n-
445 out[2] = -in[1] - in[2];
\n-
446 out[3] = in[1] + in[2];
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
464 DUNE_THROW(RangeError,
"Component out of range.");
\n-
\n-
\n-
467 }
else if (totalOrder == 2)
\n-
\n-
469 if ((order[0] == 1 && order[1] == 1) ||
\n-
470 (order[1] == 1 && order[2] == 1 && in[0] > in[1]) ||
\n-
471 (order[0] == 1 && order[2] == 1 && in[0] <=in[1]))
\n-
\n-
473 out = {1, -1, -1, 1, 0};
\n-
\n-
\n-
476 out = {0, 0, 0, 0, 0};
\n-
\n-
\n-
\n-
\n-
481 out = {0, 0, 0, 0, 0};
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
492 const R x = 2.0*in[0] + in[2] - 1.0;
\n-
493 const R y = 2.0*in[1] + in[2] - 1.0;
\n-
\n-
\n-
496 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
505 out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
\n-
506 out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
\n-
507 out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
\n-
508 out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
\n-
\n-
510 out[5] = -(y - z + 1)*(2*x + z - 1)*(y - 1);
\n-
511 out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
\n-
512 out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
\n-
513 out[8] = -(y - z + 1)*(2*x + z)*y;
\n-
514 out[9] = 2*z*(y - z - 1);
\n-
515 out[10] = -2*z*(y - z - 1);
\n-
516 out[11] = -2*z*(y - z + 1);
\n-
517 out[12] = 2*z*(y - z + 1);
\n-
518 out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
\n-
\n-
\n-
521 out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
\n-
522 out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
\n-
523 out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
\n-
524 out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
\n-
\n-
526 out[5] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
\n-
527 out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
\n-
528 out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
\n-
529 out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
\n-
530 out[9] = 2*z*(x + z - 1);
\n-
531 out[10] = -2*z*(x + z + 1);
\n-
532 out[11] = -2*z*(x + z - 1);
\n-
533 out[12] = 2*z*(x + z + 1);
\n-
534 out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
\n-
\n-
\n-
537 out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
\n-
538 out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
\n-
539 out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
\n-
540 out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
\n-
\n-
542 out[5] = (-(y - z + 1)*(2*x + z - 1)*(y - 1) - (x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x + (x + z - 1)*(y - 1)*x - (y - z + 1)*(y - 1)*x)/2;
\n-
543 out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
\n-
544 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);
\n-
545 out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
\n-
546 out[9] = -(x + 3*z - 1)*(z - y + 1);
\n-
547 out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
\n-
548 out[11] = -(x + 3*z - 1)*(y - z + 1);
\n-
549 out[12] = (x + 3*z + 1)*(y - z + 1);
\n-
550 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);
\n-
\n-
\n-
553 DUNE_THROW(RangeError,
"Component out of range.");
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
561 out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
\n-
562 out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
\n-
563 out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
\n-
564 out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
\n-
\n-
566 out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
\n-
567 out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
\n-
568 out[7] = -y*(2*x - z)*(y + z - 1);
\n-
569 out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
\n-
570 out[9] = 2*z*(y + z - 1);
\n-
571 out[10] = -2*z*(y + z - 1);
\n-
572 out[11] = -2*z*(y + z + 1);
\n-
573 out[12] = 2*z*(y + z + 1);
\n-
574 out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
\n-
\n-
\n-
577 out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
\n-
578 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
\n-
579 out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
\n-
580 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
\n-
\n-
582 out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
\n-
583 out[6] = -x*(2*y + z)*(x - z + 1);
\n-
584 out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
\n-
585 out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
\n-
586 out[9] = -2*z*(z - x + 1);
\n-
587 out[10] = -2*z*(x - z + 1);
\n-
588 out[11] = 2*z*(z - x + 1);
\n-
589 out[12] = 2*z*(x - z + 1);
\n-
590 out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
\n-
\n-
\n-
593 out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
\n-
594 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
\n-
595 out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
\n-
596 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
\n-
\n-
598 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);
\n-
599 out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
\n-
600 out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
\n-
601 out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
\n-
602 out[9] = -(z - x + 1)*(y + 3*z - 1);
\n-
603 out[10] = -(x - z + 1)*(y + 3*z - 1);
\n-
604 out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
\n-
605 out[12] = (x - z + 1)*(y + 3*z + 1);
\n-
606 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);
\n-
\n-
\n-
609 DUNE_THROW(RangeError,
"Component out of range.");
\n-
\n-
\n-
\n-
613 DUNE_THROW(NotImplemented,
"Desired derivative order is not implemented");
\n-
\n-
\n-
\n-
\n-
\n-
619 DUNE_THROW(NotImplemented,
"LagrangePyramidLocalBasis::partial for order " << k);
\n-
\n-
\n-
623 static constexpr unsigned int order ()
\n-
\n-
\n-
\n-
\n-
\n-
633 template<
unsigned int k>
\n-
634 class LagrangePyramidLocalCoefficients
\n-
\n-
\n-
638 LagrangePyramidLocalCoefficients ()
\n-
\n-
\n-
\n-
\n-
643 localKeys_[0] = LocalKey(0,0,0);
\n-
\n-
\n-
\n-
\n-
\n-
649 for (std::size_t i=0; i<size(); i++)
\n-
650 localKeys_[i] = LocalKey(i,3,0);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
657 localKeys_[0] = LocalKey(0,3,0);
\n-
658 localKeys_[1] = LocalKey(1,3,0);
\n-
659 localKeys_[2] = LocalKey(2,3,0);
\n-
660 localKeys_[3] = LocalKey(3,3,0);
\n-
661 localKeys_[4] = LocalKey(4,3,0);
\n-
\n-
\n-
664 localKeys_[5] = LocalKey(0,2,0);
\n-
665 localKeys_[6] = LocalKey(1,2,0);
\n-
666 localKeys_[7] = LocalKey(2,2,0);
\n-
667 localKeys_[8] = LocalKey(3,2,0);
\n-
668 localKeys_[9] = LocalKey(4,2,0);
\n-
669 localKeys_[10] = LocalKey(5,2,0);
\n-
670 localKeys_[11] = LocalKey(6,2,0);
\n-
671 localKeys_[12] = LocalKey(7,2,0);
\n-
\n-
\n-
674 localKeys_[13] = LocalKey(0,1,0);
\n-
\n-
\n-
\n-
\n-
\n-
680 DUNE_THROW(NotImplemented,
"LagrangePyramidLocalCoefficients for order " << k);
\n-
\n-
\n-
\n-
685 static constexpr std::size_t size ()
\n-
\n-
687 std::size_t result = 0;
\n-
688 for (
unsigned int i=0; i<=k; i++)
\n-
689 result += power(i+1,2);
\n-
\n-
\n-
\n-
694 const LocalKey& localKey (std::size_t i)
const \n-
\n-
696 return localKeys_[i];
\n-
\n-
\n-
\n-
700 std::vector<LocalKey> localKeys_;
\n-
\n-
\n-
707 template<
class LocalBasis>
\n-
708 class LagrangePyramidLocalInterpolation
\n-
\n-
\n-
\n-
719 template<
typename F,
typename C>
\n-
720 void interpolate (
const F& f, std::vector<C>& out)
const \n-
\n-
722 constexpr auto k = LocalBasis::order();
\n-
723 using D =
typename LocalBasis::Traits::DomainType;
\n-
724 using DF =
typename LocalBasis::Traits::DomainFieldType;
\n-
\n-
726 out.resize(LocalBasis::size());
\n-
\n-
\n-
\n-
\n-
731 auto center = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(0,0);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
739 for (
unsigned int i=0; i<LocalBasis::size(); i++)
\n-
\n-
741 auto vertex = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(i,3);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
750 out[0] = f( D( {0.0, 0.0, 0.0} ) );
\n-
751 out[1] = f( D( {1.0, 0.0, 0.0} ) );
\n-
752 out[2] = f( D( {0.0, 1.0, 0.0} ) );
\n-
753 out[3] = f( D( {1.0, 1.0, 0.0} ) );
\n-
754 out[4] = f( D( {0.0, 0.0, 1.0} ) );
\n-
755 out[5] = f( D( {0.0, 0.5, 0.0} ) );
\n-
756 out[6] = f( D( {1.0, 0.5, 0.0} ) );
\n-
757 out[7] = f( D( {0.5, 0.0, 0.0} ) );
\n-
758 out[8] = f( D( {0.5, 1.0, 0.0} ) );
\n-
759 out[9] = f( D( {0.0, 0.0, 0.5} ) );
\n-
760 out[10] = f( D( {0.5, 0.0, 0.5} ) );
\n-
761 out[11] = f( D( {0.0, 0.5, 0.5} ) );
\n-
762 out[12] = f( D( {0.5, 0.5, 0.5} ) );
\n-
763 out[13] = f( D( {0.5, 0.5, 0.0} ) );
\n-
\n-
\n-
\n-
\n-
768 DUNE_THROW(NotImplemented,
"LagrangePyramidLocalInterpolation not implemented for order " << k);
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
807 template<
class D,
class R,
int k>
\n-
\n-
\n-
\n-
\n-
\n-
814 Impl::LagrangePyramidLocalCoefficients<k>,
\n-
815 Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > >;
\n-
\n-
\n-
\n-
\n-
\n-
\n-
828 return coefficients_;
\n-
\n-
\n-
\n-
\n-
\n-
\n-
835 return interpolation_;
\n-
\n-
\n-
\n-
\n-
839 static constexpr std::size_t
size ()
\n-
\n-
841 return Impl::LagrangePyramidLocalBasis<D,R,k>::size();
\n-
\n-
\n-
\n-
\n-
846 static constexpr GeometryType
type ()
\n-
\n-
848 return GeometryTypes::pyramid;
\n-
\n-
\n-
\n-
\n-
852 Impl::LagrangePyramidLocalBasis<D,R,k> basis_;
\n-
853 Impl::LagrangePyramidLocalCoefficients<k> coefficients_;
\n-
854 Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > interpolation_;
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n-
\n+
\n+
\n+
\n+
![\"doxygen\"/](\"doxygen.svg\")