Dune Core Modules (2.9.1)

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