DUNE PDELab (2.7)

lagrangepyramid.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_LAGRANGE_LAGRANGEPYRAMID_HH
4#define DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPYRAMID_HH
5
6#include <array>
7#include <numeric>
8
9#include <dune/common/fmatrix.hh>
10#include <dune/common/fvector.hh>
11#include <dune/common/math.hh>
12
13#include <dune/geometry/referenceelements.hh>
14
15#include <dune/localfunctions/common/localbasis.hh>
16#include <dune/localfunctions/common/localfiniteelementtraits.hh>
17#include <dune/localfunctions/common/localinterpolation.hh>
18#include <dune/localfunctions/common/localkey.hh>
19
20namespace Dune { namespace Impl
21{
31 template<class D, class R, unsigned int k>
32 class LagrangePyramidLocalBasis
33 {
34 public:
35 using Traits = LocalBasisTraits<D,3,FieldVector<D,3>,R,1,FieldVector<R,1>,FieldMatrix<R,1,3> >;
36
39 static constexpr std::size_t size ()
40 {
41 std::size_t result = 0;
42 for (unsigned int i=0; i<=k; i++)
43 result += power(i+1,2);
44 return result;
45 }
46
48 void evaluateFunction(const typename Traits::DomainType& in,
49 std::vector<typename Traits::RangeType>& out) const
50 {
51 out.resize(size());
52
53 // Specialization for zero-order case
54 if (k==0)
55 {
56 out[0] = 1;
57 return;
58 }
59
60 if (k==1)
61 {
62 if(in[0] > in[1])
63 {
64 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[1]);
65 out[1] = in[0]*(1-in[1])-in[2]*in[1];
66 out[2] = (1-in[0])*in[1]-in[2]*in[1];
67 out[3] = in[0]*in[1]+in[2]*in[1];
68 }
69 else
70 {
71 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[0]);
72 out[1] = in[0]*(1-in[1])-in[2]*in[0];
73 out[2] = (1-in[0])*in[1]-in[2]*in[0];
74 out[3] = in[0]*in[1]+in[2]*in[0];
75 }
76
77 out[4] = in[2];
78
79 return;
80 }
81
82 if (k==2)
83 {
84 // transform to reference element with base [-1,1]^2
85 const R x = 2.0*in[0] + in[2] - 1.0;
86 const R y = 2.0*in[1] + in[2] - 1.0;
87 const R z = in[2];
88
89 if (x > y)
90 {
91 // vertices
92 out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
93 out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
94 out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
95 out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
96 out[4] = z*(2*z - 1);
97
98 // lower edges
99 out[5] = -0.5*(y - z + 1)*(x + z - 1)*(y - 1)*x;
100 out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
101 out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
102 out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
103
104 // upper edges
105 out[9] = z*(x + z - 1)*(y - z - 1);
106 out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
107 out[11] = -z*(y - z + 1)*(x + z - 1);
108 out[12] = z*(y - z + 1)*(x + z + 1);
109
110 // base face
111 out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
112 }
113 else
114 {
115 // vertices
116 out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
117 out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
118 out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
119 out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
120 out[4] = z*(2*z - 1);
121
122 // lower edges
123 out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
124 out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
125 out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
126 out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
127
128 // upper edges
129 out[9] = z*(y + z - 1)*(x - z - 1);
130 out[10] = -z*(x - z + 1)*(y + z - 1);
131 out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
132 out[12] = z*(x - z + 1)*(y + z + 1);
133
134 // base face
135 out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
136 }
137
138 return;
139 }
140
141 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateFunction for order " << k);
142 }
143
149 void evaluateJacobian(const typename Traits::DomainType& in,
150 std::vector<typename Traits::JacobianType>& out) const
151 {
152 out.resize(size());
153
154 // Specialization for k==0
155 if (k==0)
156 {
157 std::fill(out[0][0].begin(), out[0][0].end(), 0);
158 return;
159 }
160
161 if (k==1)
162 {
163 if(in[0] > in[1])
164 {
165 out[0][0] = {-1 + in[1], -1 + in[0] + in[2], -1 + in[1]};
166 out[1][0] = { 1 - in[1], -in[0] - in[2], -in[1]};
167 out[2][0] = { -in[1], 1 - in[0] - in[2], -in[1]};
168 out[3][0] = { in[1], in[0] + in[2], in[1]};
169 }
170 else
171 {
172 out[0][0] = {-1 + in[1] + in[2], -1 + in[0], -1 + in[0]};
173 out[1][0] = { 1 - in[1] - in[2], -in[0], -in[0]};
174 out[2][0] = { -in[1] - in[2], 1 - in[0], -in[0]};
175 out[3][0] = { in[1] + in[2], in[0], in[0]};
176 }
177
178 out[4][0] = {0, 0, 1};
179 return;
180 }
181
182 if (k==2)
183 {
184 // transform to reference element with base [-1,1]^2
185 const R x = 2.0*in[0] + in[2] - 1.0;
186 const R y = 2.0*in[1] + in[2] - 1.0;
187 const R z = in[2];
188
189 // transformation of the gradient leads to a multiplication
190 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
191 if (x > y)
192 {
193 // vertices
194 out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
195 out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
196 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
197 + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
198 + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
199
200 out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
201 + (x + z)*(y - z)*(-y + z + 1)) - z);
202 out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
203 + (x + z)*(y - z)*(-(x + z + 1))) + z);
204 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
205 - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
206 - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
207 + (x + z)*(y - z)*(x - y + 2*z - 2))
208 - (x - y);
209
210 out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
211 out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
212 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
213 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
214 + (x + z)*(y - z)*(y - x - 2*z + 2));
215
216 out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
217 out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
218 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
219 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
220 + (y - z)*(x + z)*(y - x - 2*z));
221
222 out[4][0][0] = 0;
223 out[4][0][1] = 0;
224 out[4][0][2] = 4*z - 1;
225
226 // lower edges
227 out[5][0][0] = -(y - z + 1)*(y - 1)*(2*x + z - 1);
228 out[5][0][1] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
229 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
230 + 0.5*(x + z - 1)*(y - 1)*x - 0.5*(y - z + 1)*(y - 1)*x;
231
232 out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
233 out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
234 + (y - z + 1)*((x + z + 1)*x + 2*z));
235 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
236 - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
237 + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
238
239 out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
240 + (x + z - 1)*((y - z - 1)*y + 2*z));
241 out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
242 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
243 - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
244 + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
245
246 out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
247 out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
248 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
249 - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
250
251 // upper edges
252 out[9][0][0] = 2*z*(y - z - 1);
253 out[9][0][1] = 2*z*(x + z - 1);
254 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
255 + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
256
257 out[10][0][0] = -2*z*(y - z - 1);
258 out[10][0][1] = -2*z*(x + z + 1);
259 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
260 - ((x + z + 1)*(y - z - 1) + 4*z)
261 - z*(-x + y - 2*z + 2);
262
263 out[11][0][0] = -2*z*(y - z + 1);
264 out[11][0][1] = -2*z*(x + z - 1);
265 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
266 - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
267
268 out[12][0][0] = 2*z*(y - z + 1);
269 out[12][0][1] = 2*z*(x + z + 1);
270 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
271 + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
272
273 // base face
274 out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
275 + (y - z + 1)*(x + z - 1)*(y - 1 + z));
276 out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
277 + (y - z + 1)*(x + z - 1)*(x + 1 - z));
278 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
279 + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
280 + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
281 }
282 else
283 {
284 // vertices
285 out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
286 out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
287 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
288 + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
289 + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
290
291 out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
292 out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
293 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
294 - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
295 + (x - z)*(y + z)*(-x + y + 2*z - 2));
296
297 out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
298 + (x - z)*(y + z)*(y + z + 1) + 4*z);
299 out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
300 + (x - z)*(y + z)*(x - z - 1) - 4*z);
301 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
302 + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
303 + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
304
305 out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
306 out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
307 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
308 + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
309 + (y + z)*(x - z)*(x - y - 2*z));
310
311 out[4][0][0] = 0;
312 out[4][0][1] = 0;
313 out[4][0][2] = 4*z - 1;
314
315 // lower edges
316 out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
317 out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
318 + (y + z - 1)*((x - z - 1)*x + 2*z));
319 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
320 - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
321 + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
322
323 out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
324 out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
325 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
326 - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
327
328 out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
329 out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
330 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
331 - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
332
333 out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
334 + (x - z + 1)*((y + z + 1)*y + 2*z));
335 out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
336 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
337 - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
338 + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
339
340 // upper edges
341 out[9][0][0] = 2*z*(y + z - 1);
342 out[9][0][1] = 2*z*(x - z - 1);
343 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
344 + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
345
346 out[10][0][0] = -2*z*(y + z - 1);
347 out[10][0][1] = -2*z*(x - z + 1);
348 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
349 - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
350
351 out[11][0][0] = -2*z*(y + z + 1);
352 out[11][0][1] = -2*z*(x - z - 1);
353 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
354 - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
355
356 out[12][0][0] = 2*z*(y + z + 1);
357 out[12][0][1] = 2*z*(x - z + 1);
358 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
359 + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
360
361 // base face
362 out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
363 + (x - z + 1)*(y + z - 1)*(y + 1 - z));
364 out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
365 + (x - z + 1)*(y + z - 1)*(x - 1 + z));
366 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
367 + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
368 + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
369 }
370
371 return;
372 }
373
374 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateJacobian for order " << k);
375 }
376
383 void partial(const std::array<unsigned int,3>& order,
384 const typename Traits::DomainType& in,
385 std::vector<typename Traits::RangeType>& out) const
386 {
387 auto totalOrder = std::accumulate(order.begin(), order.end(), 0);
388
389 out.resize(size());
390
391 if (totalOrder == 0)
392 {
393 evaluateFunction(in, out);
394 return;
395 }
396
397 if (k==0)
398 {
399 out[0] = 0;
400 return;
401 }
402
403 if (k==1)
404 {
405 if (totalOrder == 1)
406 {
407 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
408 if (in[0] > in[1])
409 {
410 switch (direction)
411 {
412 case 0:
413 out[0] = -1 + in[1];
414 out[1] = 1 - in[1];
415 out[2] = -in[1];
416 out[3] = in[1];
417 out[4] = 0;
418 break;
419 case 1:
420 out[0] = -1 + in[0] + in[2];
421 out[1] = -in[0] - in[2];
422 out[2] = 1 - in[0] - in[2];
423 out[3] = in[0]+in[2];
424 out[4] = 0;
425 break;
426 case 2:
427 out[0] = -1 + in[1];
428 out[1] = -in[1];
429 out[2] = -in[1];
430 out[3] = in[1];
431 out[4] = 1;
432 break;
433 default:
434 DUNE_THROW(RangeError, "Component out of range.");
435 }
436 }
437 else /* (in[0] <= in[1]) */
438 {
439 switch (direction)
440 {
441 case 0:
442 out[0] = -1 + in[1] + in[2];
443 out[1] = 1 - in[1] - in[2];
444 out[2] = -in[1] - in[2];
445 out[3] = in[1] + in[2];
446 out[4] = 0;
447 break;
448 case 1:
449 out[0] = -1 + in[0];
450 out[1] = -in[0];
451 out[2] = 1 - in[0];
452 out[3] = in[0];
453 out[4] = 0;
454 break;
455 case 2:
456 out[0] = -1 + in[0];
457 out[1] = -in[0];
458 out[2] = -in[0];
459 out[3] = in[0];
460 out[4] = 1;
461 break;
462 default:
463 DUNE_THROW(RangeError, "Component out of range.");
464 }
465 }
466 } else if (totalOrder == 2)
467 {
468 if ((order[0] == 1 && order[1] == 1) ||
469 (order[1] == 1 && order[2] == 1 && in[0] > in[1]) ||
470 (order[0] == 1 && order[2] == 1 && in[0] <=in[1]))
471 {
472 out = {1, -1, -1, 1, 0};
473 } else
474 {
475 out = {0, 0, 0, 0, 0};
476 }
477
478 } else
479 {
480 out = {0, 0, 0, 0, 0};
481 }
482
483 return;
484 }
485
486 if (k==2)
487 {
488 if (totalOrder == 1)
489 {
490 // transform to reference element with base [-1,1]^2
491 const R x = 2.0*in[0] + in[2] - 1.0;
492 const R y = 2.0*in[1] + in[2] - 1.0;
493 const R z = in[2];
494
495 auto const direction = std::distance(order.begin(), std::find(order.begin(), order.end(), 1));
496
497 // transformation of the gradient leads to a multiplication
498 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
499 if (x > y)
500 {
501 switch (direction)
502 {
503 case 0:
504 out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
505 out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
506 out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
507 out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
508 out[4] = 0;
509 out[5] = -(y - z + 1)*(2*x + z - 1)*(y - 1);
510 out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
511 out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
512 out[8] = -(y - z + 1)*(2*x + z)*y;
513 out[9] = 2*z*(y - z - 1);
514 out[10] = -2*z*(y - z - 1);
515 out[11] = -2*z*(y - z + 1);
516 out[12] = 2*z*(y - z + 1);
517 out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
518 break;
519 case 1:
520 out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
521 out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
522 out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
523 out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
524 out[4] = 0;
525 out[5] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
526 out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
527 out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
528 out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
529 out[9] = 2*z*(x + z - 1);
530 out[10] = -2*z*(x + z + 1);
531 out[11] = -2*z*(x + z - 1);
532 out[12] = 2*z*(x + z + 1);
533 out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
534 break;
535 case 2:
536 out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
537 out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
538 out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
539 out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
540 out[4] = 4*z - 1;
541 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;
542 out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
543 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);
544 out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
545 out[9] = -(x + 3*z - 1)*(z - y + 1);
546 out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
547 out[11] = -(x + 3*z - 1)*(y - z + 1);
548 out[12] = (x + 3*z + 1)*(y - z + 1);
549 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);
550 break;
551 default:
552 DUNE_THROW(RangeError, "Component out of range.");
553 }
554 }
555 else // x <= y
556 {
557 switch (direction)
558 {
559 case 0:
560 out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
561 out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
562 out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
563 out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
564 out[4] = 0;
565 out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
566 out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
567 out[7] = -y*(2*x - z)*(y + z - 1);
568 out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
569 out[9] = 2*z*(y + z - 1);
570 out[10] = -2*z*(y + z - 1);
571 out[11] = -2*z*(y + z + 1);
572 out[12] = 2*z*(y + z + 1);
573 out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
574 break;
575 case 1:
576 out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
577 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
578 out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
579 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
580 out[4] = 0;
581 out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
582 out[6] = -x*(2*y + z)*(x - z + 1);
583 out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
584 out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
585 out[9] = -2*z*(z - x + 1);
586 out[10] = -2*z*(x - z + 1);
587 out[11] = 2*z*(z - x + 1);
588 out[12] = 2*z*(x - z + 1);
589 out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
590 break;
591 case 2:
592 out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
593 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
594 out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
595 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
596 out[4] = 4*z - 1;
597 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);
598 out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
599 out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
600 out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
601 out[9] = -(z - x + 1)*(y + 3*z - 1);
602 out[10] = -(x - z + 1)*(y + 3*z - 1);
603 out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
604 out[12] = (x - z + 1)*(y + 3*z + 1);
605 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);
606 break;
607 default:
608 DUNE_THROW(RangeError, "Component out of range.");
609 }
610 }
611 } else {
612 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
613 }
614
615 return;
616 }
617
618 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::partial for order " << k);
619 }
620
622 static constexpr unsigned int order ()
623 {
624 return k;
625 }
626 };
627
632 template<unsigned int k>
633 class LagrangePyramidLocalCoefficients
634 {
635 public:
637 LagrangePyramidLocalCoefficients ()
638 : localKeys_(size())
639 {
640 if (k==0)
641 {
642 localKeys_[0] = LocalKey(0,0,0);
643 return;
644 }
645
646 if (k==1)
647 {
648 for (std::size_t i=0; i<size(); i++)
649 localKeys_[i] = LocalKey(i,3,0);
650 return;
651 }
652
653 if (k==2)
654 {
655 // Vertex shape functions
656 localKeys_[0] = LocalKey(0,3,0);
657 localKeys_[1] = LocalKey(1,3,0);
658 localKeys_[2] = LocalKey(2,3,0);
659 localKeys_[3] = LocalKey(3,3,0);
660 localKeys_[4] = LocalKey(4,3,0);
661
662 // Edge shape functions
663 localKeys_[5] = LocalKey(0,2,0);
664 localKeys_[6] = LocalKey(1,2,0);
665 localKeys_[7] = LocalKey(2,2,0);
666 localKeys_[8] = LocalKey(3,2,0);
667 localKeys_[9] = LocalKey(4,2,0);
668 localKeys_[10] = LocalKey(5,2,0);
669 localKeys_[11] = LocalKey(6,2,0);
670 localKeys_[12] = LocalKey(7,2,0);
671
672 // base face shape function
673 localKeys_[13] = LocalKey(0,1,0);
674
675 return;
676 }
677
678 // No general case
679 DUNE_THROW(NotImplemented, "LagrangePyramidLocalCoefficients for order " << k);
680
681 }
682
684 static constexpr std::size_t size ()
685 {
686 std::size_t result = 0;
687 for (unsigned int i=0; i<=k; i++)
688 result += power(i+1,2);
689 return result;
690 }
691
693 const LocalKey& localKey (std::size_t i) const
694 {
695 return localKeys_[i];
696 }
697
698 private:
699 std::vector<LocalKey> localKeys_;
700 };
701
706 template<class LocalBasis>
707 class LagrangePyramidLocalInterpolation
708 {
709 public:
710
718 template<typename F, typename C>
719 void interpolate (const F& ff, std::vector<C>& out) const
720 {
721 constexpr auto k = LocalBasis::order();
722 using D = typename LocalBasis::Traits::DomainType;
723 using DF = typename LocalBasis::Traits::DomainFieldType;
724
725 auto&& f = Impl::makeFunctionWithCallOperator<D>(ff);
726
727 out.resize(LocalBasis::size());
728
729 // Specialization for zero-order case
730 if (k==0)
731 {
732 auto center = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(0,0);
733 out[0] = f(center);
734 return;
735 }
736
737 // Specialization for first-order case
738 if (k==1)
739 {
740 for (unsigned int i=0; i<LocalBasis::size(); i++)
741 {
742 auto vertex = ReferenceElements<DF,3>::general(GeometryTypes::pyramid).position(i,3);
743 out[i] = f(vertex);
744 }
745 return;
746 }
747
748 // Specialization for second-order case
749 if (k==2)
750 {
751 out[0] = f( D( {0.0, 0.0, 0.0} ) );
752 out[1] = f( D( {1.0, 0.0, 0.0} ) );
753 out[2] = f( D( {0.0, 1.0, 0.0} ) );
754 out[3] = f( D( {1.0, 1.0, 0.0} ) );
755 out[4] = f( D( {0.0, 0.0, 1.0} ) );
756 out[5] = f( D( {0.0, 0.5, 0.0} ) );
757 out[6] = f( D( {1.0, 0.5, 0.0} ) );
758 out[7] = f( D( {0.5, 0.0, 0.0} ) );
759 out[8] = f( D( {0.5, 1.0, 0.0} ) );
760 out[9] = f( D( {0.0, 0.0, 0.5} ) );
761 out[10] = f( D( {0.5, 0.0, 0.5} ) );
762 out[11] = f( D( {0.0, 0.5, 0.5} ) );
763 out[12] = f( D( {0.5, 0.5, 0.5} ) );
764 out[13] = f( D( {0.5, 0.5, 0.0} ) );
765
766 return;
767 }
768
769 DUNE_THROW(NotImplemented, "LagrangePyramidLocalInterpolation not implemented for order " << k);
770 }
771
772 };
773
774} } // namespace Dune::Impl
775
776namespace Dune
777{
808 template<class D, class R, int k>
809 class LagrangePyramidLocalFiniteElement
810 {
811 public:
814 using Traits = LocalFiniteElementTraits<Impl::LagrangePyramidLocalBasis<D,R,k>,
815 Impl::LagrangePyramidLocalCoefficients<k>,
816 Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > >;
817
823 LagrangePyramidLocalFiniteElement() {}
824
827 const typename Traits::LocalBasisType& localBasis () const
828 {
829 return basis_;
830 }
831
834 const typename Traits::LocalCoefficientsType& localCoefficients () const
835 {
836 return coefficients_;
837 }
838
841 const typename Traits::LocalInterpolationType& localInterpolation () const
842 {
843 return interpolation_;
844 }
845
847 static constexpr std::size_t size ()
848 {
849 return Impl::LagrangePyramidLocalBasis<D,R,k>::size();
850 }
851
854 static constexpr GeometryType type ()
855 {
856 return GeometryTypes::pyramid;
857 }
858
859 private:
860 Impl::LagrangePyramidLocalBasis<D,R,k> basis_;
861 Impl::LagrangePyramidLocalCoefficients<k> coefficients_;
862 Impl::LagrangePyramidLocalInterpolation<Impl::LagrangePyramidLocalBasis<D,R,k> > interpolation_;
863 };
864
865} // namespace Dune
866
867#endif // DUNE_LOCALFUNCTIONS_LAGRANGE_LAGRANGEPYRAMID_HH
Dune::GeometryType
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:280
Dune::LagrangePyramidLocalFiniteElement
Lagrange finite element for 3d pyramids with compile-time polynomial order.
Definition: lagrangepyramid.hh:810
Dune::LagrangePyramidLocalFiniteElement::localCoefficients
const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangepyramid.hh:834
Dune::LagrangePyramidLocalFiniteElement::size
static constexpr std::size_t size()
The number of shape functions.
Definition: lagrangepyramid.hh:847
Dune::LagrangePyramidLocalFiniteElement::LagrangePyramidLocalFiniteElement
LagrangePyramidLocalFiniteElement()
Default constructor.
Definition: lagrangepyramid.hh:823
Dune::LagrangePyramidLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangepyramid.hh:854
Dune::LagrangePyramidLocalFiniteElement::localBasis
const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangepyramid.hh:827
Dune::LagrangePyramidLocalFiniteElement::localInterpolation
const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangepyramid.hh:841
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Dune::GeometryTypes::pyramid
constexpr GeometryType pyramid
GeometryType representing a 3D pyramid.
Definition: type.hh:828
Dune::GeometryTypes::vertex
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:798
Dune::PDELab::interpolate
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:174
Dune::Hybrid::accumulate
T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:290
math.hh
Some useful basic math stuff.
Dune
Dune namespace.
Definition: alignedallocator.hh:14
Dune::power
constexpr Mantissa power(Mantissa m, Exponent p)
Power method for integer exponents.
Definition: math.hh:73
Dune::Geo::ReferenceElements::general
static const ReferenceElement & general(const GeometryType &type)
get general reference elements
Definition: referenceelements.hh:196
Dune::LocalBasisTraits::DomainType
D DomainType
domain type
Definition: localbasis.hh:43
Dune::LocalFiniteElementTraits
traits helper struct
Definition: localfiniteelementtraits.hh:11
Dune::LocalFiniteElementTraits::LocalBasisType
LB LocalBasisType
Definition: localfiniteelementtraits.hh:14
Dune::LocalFiniteElementTraits::LocalCoefficientsType
LC LocalCoefficientsType
Definition: localfiniteelementtraits.hh:18
Dune::LocalFiniteElementTraits::LocalInterpolationType
LI LocalInterpolationType
Definition: localfiniteelementtraits.hh:22
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)