DUNE-FEM (unstable)

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