DUNE PDELab (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{
23
24 // The traits provide static or dynamic order and size information
25 template<int compileTimeOrder>
26 struct LagrangePyramidOrderTraits;
27
28 // The traits provide static order and size information
29 template<int compileTimeOrder>
30 requires (compileTimeOrder >= 0)
31 struct LagrangePyramidOrderTraits<compileTimeOrder>
32 {
33 static constexpr bool is_static_order = true;
34
35 constexpr LagrangePyramidOrderTraits(int /*runTimeOrder*/ = compileTimeOrder)
36 {
37 static_assert((0 <= compileTimeOrder) and (compileTimeOrder <= 2), "LagrangePyramid: Only order 0,1, and 2 are supported");
38 }
39
43 static constexpr unsigned int order()
44 {
45 return compileTimeOrder;
46 }
47
51 static constexpr std::size_t size()
52 {
53 std::size_t result = 0;
54 for (int i=0; i<=compileTimeOrder; i++)
55 result += power(i+1,2);
56 return result;
57 }
58 };
59
60
61 // this specialization is for the dynamic order case
62 template<>
63 struct LagrangePyramidOrderTraits<-1>
64 {
65 unsigned int runTimeOrder_;
66 std::size_t size_;
67
68 static constexpr bool is_static_order = false;
69
70 constexpr explicit LagrangePyramidOrderTraits(int runTimeOrder)
71 : runTimeOrder_(runTimeOrder >= 0 ? (unsigned int)(runTimeOrder) : 0u)
72 , size_(0)
73 {
74 if ((runTimeOrder < 0) or (runTimeOrder > 2))
75 DUNE_THROW(Dune::InvalidStateException, "LagrangePyramid: Only order 0,1, and 2 are supported");
76 size_ = 0;
77 for (unsigned int i=0; i<=runTimeOrder_; i++)
78 size_ += power(i+1,2);
79 }
80
84 constexpr unsigned int order() const
85 {
86 return runTimeOrder_;
87 }
88
92 constexpr std::size_t size() const
93 {
94 return size_;
95 }
96 };
97
98
99
109 template<class D, class R, int compileTimeOrder>
110 class LagrangePyramidLocalBasis
111 : private LagrangePyramidOrderTraits<compileTimeOrder>
112 {
113 using OrderTraits = LagrangePyramidOrderTraits<compileTimeOrder>;
114
115 public:
116 using Traits = LocalBasisTraits<D,3,FieldVector<D,3>,R,1,FieldVector<R,1>,FieldMatrix<R,1,3> >;
117
119 constexpr LagrangePyramidLocalBasis() requires (OrderTraits::is_static_order)
120 : OrderTraits()
121 {}
122
124 explicit constexpr LagrangePyramidLocalBasis(OrderTraits orderTraits)
125 : OrderTraits(orderTraits)
126 {}
127
128 using OrderTraits::order;
129 using OrderTraits::size;
130
132 constexpr void evaluateFunction(const typename Traits::DomainType& in,
133 std::vector<typename Traits::RangeType>& out) const
134 {
135 out.resize(size());
136
137 if (order()==0)
138 {
139 out[0] = 1;
140 return;
141 }
142
143 if (order()==1)
144 {
145 if(in[0] > in[1])
146 {
147 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[1]);
148 out[1] = in[0]*(1-in[1])-in[2]*in[1];
149 out[2] = (1-in[0])*in[1]-in[2]*in[1];
150 out[3] = in[0]*in[1]+in[2]*in[1];
151 }
152 else
153 {
154 out[0] = (1-in[0])*(1-in[1])-in[2]*(1-in[0]);
155 out[1] = in[0]*(1-in[1])-in[2]*in[0];
156 out[2] = (1-in[0])*in[1]-in[2]*in[0];
157 out[3] = in[0]*in[1]+in[2]*in[0];
158 }
159
160 out[4] = in[2];
161
162 return;
163 }
164
165 if (order()==2)
166 {
167 // transform to reference element with base [-1,1]^2
168 const R x = 2.0*in[0] + in[2] - 1.0;
169 const R y = 2.0*in[1] + in[2] - 1.0;
170 const R z = in[2];
171
172 if (x > y)
173 {
174 // vertices
175 out[0] = 0.25*(x + z)*(x + z - 1)*(y - z - 1)*(y - z);
176 out[1] = -0.25*(x + z)*(y - z)*((x + z + 1)*(-y + z + 1) - 4*z) - z*(x - y);
177 out[2] = 0.25*(x + z)*(y - z)*(y - z + 1)*(x + z - 1);
178 out[3] = 0.25*(y - z)*(x + z)*(y - z + 1)*(x + z + 1);
179 out[4] = z*(2*z - 1);
180
181 // lower edges
182 out[5] = -0.5*(y - z + 1)*(x + z - 1)*(y - 1)*x;
183 out[6] = -0.5*(y - z + 1)*(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1));
184 out[7] = -0.5*(x + z - 1)*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1));
185 out[8] = -0.5*(y - z + 1)*(x + z - 1)*(x + 1)*y;
186
187 // upper edges
188 out[9] = z*(x + z - 1)*(y - z - 1);
189 out[10] = -z*((x + z + 1)*(y - z - 1) + 4*z);
190 out[11] = -z*(y - z + 1)*(x + z - 1);
191 out[12] = z*(y - z + 1)*(x + z + 1);
192
193 // base face
194 out[13] = (y - z + 1)*(x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1));
195 }
196 else
197 {
198 // vertices
199 out[0] = 0.25*(y + z)*(y + z - 1)*(x - z - 1)*(x - z);
200 out[1] = -0.25*(x - z)*(y + z)*(x - z + 1)*(-y - z + 1);
201 out[2] = 0.25*(x - z)*(y + z)*((x - z - 1)*(y + z + 1) + 4*z) + z*(x - y);
202 out[3] = 0.25*(y + z)*(x - z)*(x - z + 1)*(y + z + 1);
203 out[4] = z*(2*z - 1);
204
205 // lower edges
206 out[5] = -0.5*(y + z - 1)*(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1));
207 out[6] = -0.5*(x - z + 1)*(y + z - 1)*(y + 1)*x;
208 out[7] = -0.5*(x - z + 1)*(y + z - 1)*(x - 1)*y;
209 out[8] = -0.5*(x - z + 1)*(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1));
210
211 // upper edges
212 out[9] = z*(y + z - 1)*(x - z - 1);
213 out[10] = -z*(x - z + 1)*(y + z - 1);
214 out[11] = -z*((y + z + 1)*(x - z - 1) + 4*z);
215 out[12] = z*(x - z + 1)*(y + z + 1);
216
217 // base face
218 out[13] = (x - z + 1)*(y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1));
219 }
220
221 return;
222 }
223
224 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateFunction for order " << order());
225 }
226
232 constexpr void evaluateJacobian(const typename Traits::DomainType& in,
233 std::vector<typename Traits::JacobianType>& out) const
234 {
235 out.resize(size());
236
237 if (order()==0)
238 {
239 std::fill(out[0][0].begin(), out[0][0].end(), 0);
240 return;
241 }
242
243 if (order()==1)
244 {
245 if(in[0] > in[1])
246 {
247 out[0][0] = {-1 + in[1], -1 + in[0] + in[2], -1 + in[1]};
248 out[1][0] = { 1 - in[1], -in[0] - in[2], -in[1]};
249 out[2][0] = { -in[1], 1 - in[0] - in[2], -in[1]};
250 out[3][0] = { in[1], in[0] + in[2], in[1]};
251 }
252 else
253 {
254 out[0][0] = {-1 + in[1] + in[2], -1 + in[0], -1 + in[0]};
255 out[1][0] = { 1 - in[1] - in[2], -in[0], -in[0]};
256 out[2][0] = { -in[1] - in[2], 1 - in[0], -in[0]};
257 out[3][0] = { in[1] + in[2], in[0], in[0]};
258 }
259
260 out[4][0] = {0, 0, 1};
261 return;
262 }
263
264 if (order()==2)
265 {
266 // transform to reference element with base [-1,1]^2
267 const R x = 2.0*in[0] + in[2] - 1.0;
268 const R y = 2.0*in[1] + in[2] - 1.0;
269 const R z = in[2];
270
271 // transformation of the gradient leads to a multiplication
272 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
273 if (x > y)
274 {
275 // vertices
276 out[0][0][0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
277 out[0][0][1] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
278 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
279 + 0.25*((2*x + 2*z - 1)*(y - z - 1)*(y - z)
280 + (x + z)*(x + z - 1)*(-2*y + 2*z + 1));
281
282 out[1][0][0] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
283 + (x + z)*(y - z)*(-y + z + 1)) - z);
284 out[1][0][1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
285 + (x + z)*(y - z)*(-(x + z + 1))) + z);
286 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
287 - 0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z)
288 - (x + z)*((x + z + 1)*(-y + z + 1) - 4*z)
289 + (x + z)*(y - z)*(x - y + 2*z - 2))
290 - (x - y);
291
292 out[2][0][0] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
293 out[2][0][1] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
294 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
295 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z - 1)
296 + (x + z)*(y - z)*(y - x - 2*z + 2));
297
298 out[3][0][0] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
299 out[3][0][1] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
300 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
301 + 0.25*((y - x - 2*z)*(y - z + 1)*(x + z + 1)
302 + (y - z)*(x + z)*(y - x - 2*z));
303
304 out[4][0][0] = 0;
305 out[4][0][1] = 0;
306 out[4][0][2] = 4*z - 1;
307
308 // lower edges
309 out[5][0][0] = -(y - z + 1)*(y - 1)*(2*x + z - 1);
310 out[5][0][1] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
311 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
312 + 0.5*(x + z - 1)*(y - 1)*x - 0.5*(y - z + 1)*(y - 1)*x;
313
314 out[6][0][0] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
315 out[6][0][1] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1)
316 + (y - z + 1)*((x + z + 1)*x + 2*z));
317 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
318 - 0.5*(-(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1))
319 + (y - z + 1)*(((y - 1)*x - 1) + 2*y + 1));
320
321 out[7][0][0] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
322 + (x + z - 1)*((y - z - 1)*y + 2*z));
323 out[7][0][1] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
324 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
325 - 0.5*(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1)
326 + (x + z - 1)*((-(x + 1)*y - 1) + 2*x + 1));
327
328 out[8][0][0] = -(y - z + 1)*(2*x + z)*y;
329 out[8][0][1] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
330 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
331 - 0.5*(-x + y - 2*z + 2)*(x + 1)*y;
332
333 // upper edges
334 out[9][0][0] = 2*z*(y - z - 1);
335 out[9][0][1] = 2*z*(x + z - 1);
336 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
337 + (x + z - 1)*(y - z - 1) + z*(-x + y - 2*z);
338
339 out[10][0][0] = -2*z*(y - z - 1);
340 out[10][0][1] = -2*z*(x + z + 1);
341 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
342 - ((x + z + 1)*(y - z - 1) + 4*z)
343 - z*(-x + y - 2*z + 2);
344
345 out[11][0][0] = -2*z*(y - z + 1);
346 out[11][0][1] = -2*z*(x + z - 1);
347 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
348 - (y - z + 1)*(x + z - 1) - z*(-x + y - 2*z + 2);
349
350 out[12][0][0] = 2*z*(y - z + 1);
351 out[12][0][1] = 2*z*(x + z + 1);
352 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
353 + (y - z + 1)*(x + z + 1) + z*(-x + y - 2*z);
354
355 // base face
356 out[13][0][0] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
357 + (y - z + 1)*(x + z - 1)*(y - 1 + z));
358 out[13][0][1] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1))
359 + (y - z + 1)*(x + z - 1)*(x + 1 - z));
360 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
361 + ((-x + y - 2*z + 2)*((y - 1)*(x + 1) + z*(x - y + z + 1))
362 + (y - z + 1)*(x + z - 1)*(x - y + 2*z + 1));
363 }
364 else
365 {
366 // vertices
367 out[0][0][0] = 0.5*(y + z)*(y + z - 1)*(2*x - 2*z - 1);
368 out[0][0][1] = 0.5*(2*y + 2*z - 1)*(x - z - 1)*(x - z);
369 out[0][0][2] = 0.5*(out[0][0][0] + out[0][0][1])
370 + 0.25*((2*y + 2*z - 1)*(x - z - 1)*(x - z)
371 + (y + z)*(y + z - 1)*(-2*x + 2*z + 1));
372
373 out[1][0][0] = -0.5*(y + z)*(2*x - 2*z + 1)*(-y - z + 1);
374 out[1][0][1] = -0.5*(x - z)*(x - z + 1)*(-2*y - 2*z + 1);
375 out[1][0][2] = 0.5*(out[1][0][0] + out[1][0][1])
376 - 0.25*((x - y - 2*z)*(x - z + 1)*(-y - z + 1)
377 + (x - z)*(y + z)*(-x + y + 2*z - 2));
378
379 out[2][0][0] = 0.5*((y + z)*((x - z - 1)*(y + z + 1) + 4*z)
380 + (x - z)*(y + z)*(y + z + 1) + 4*z);
381 out[2][0][1] = 0.5*((x - z)*((x - z - 1)*(y + z + 1) + 4*z)
382 + (x - z)*(y + z)*(x - z - 1) - 4*z);
383 out[2][0][2] = 0.5*(out[2][0][0] + out[2][0][1])
384 + 0.25*((x - y - 2*z)*((x - z - 1)*(y + z + 1) + 4*z)
385 + (x - z)*(y + z)*(x - y - 2*z + 2) + 4*(x - y));
386
387 out[3][0][0] = 0.5*(y + z)*(2*x - 2*z + 1)*(y + z + 1);
388 out[3][0][1] = 0.5*(x - z)*(x - z + 1)*(2*y + 2*z + 1);
389 out[3][0][2] = 0.5*(out[3][0][0] + out[3][0][1])
390 + 0.25*((x - y - 2*z)*(x - z + 1)*(y + z + 1)
391 + (y + z)*(x - z)*(x - y - 2*z));
392
393 out[4][0][0] = 0;
394 out[4][0][1] = 0;
395 out[4][0][2] = 4*z - 1;
396
397 // lower edges
398 out[5][0][0] = -(y + z - 1)*(2*x - z - 1)*(y + 1);
399 out[5][0][1] = -(((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1)
400 + (y + z - 1)*((x - z - 1)*x + 2*z));
401 out[5][0][2] = 0.5*(out[5][0][0] + out[5][0][1])
402 - 0.5*((((x - z - 1)*(y + 1)*x - z) + z*(2*y + 1))
403 + (y + z - 1)*((-(y + 1)*x - 1) + 2*y + 1));
404
405 out[6][0][0] = -(2*x - z + 1)*(y + z - 1)*(y + 1);
406 out[6][0][1] = -(x - z + 1)*(2*y + z)*x;
407 out[6][0][2] = 0.5*(out[6][0][0] + out[6][0][1])
408 - 0.5*(x - y - 2*z + 2)*(y + 1)*x;
409
410 out[7][0][0] = -(2*x - z)*(y + z - 1)*y;
411 out[7][0][1] = -(x - z + 1)*(2*y + z - 1)*(x - 1);
412 out[7][0][2] = 0.5*(out[7][0][0] + out[7][0][1])
413 - 0.5*(x - y - 2*z + 2)*(x - 1)*y;
414
415 out[8][0][0] = -(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1)
416 + (x - z + 1)*((y + z + 1)*y + 2*z));
417 out[8][0][1] = -(x - z + 1)*(2*y + z + 1)*(x - 1);
418 out[8][0][2] = 0.5*(out[8][0][0] + out[8][0][1])
419 - 0.5*(-(((y + z + 1)*(x - 1)*y - z) + z*(2*x + 1))
420 + (x - z + 1)*(((x - 1)*y - 1) + 2*x + 1));
421
422 // upper edges
423 out[9][0][0] = 2*z*(y + z - 1);
424 out[9][0][1] = 2*z*(x - z - 1);
425 out[9][0][2] = 0.5*(out[9][0][0] + out[9][0][1])
426 + (y + z - 1)*(x - z - 1) + z*(x - y - 2*z);
427
428 out[10][0][0] = -2*z*(y + z - 1);
429 out[10][0][1] = -2*z*(x - z + 1);
430 out[10][0][2] = 0.5*(out[10][0][0] + out[10][0][1])
431 - (x - z + 1)*(y + z - 1) - z*(x - y - 2*z + 2);
432
433 out[11][0][0] = -2*z*(y + z + 1);
434 out[11][0][1] = -2*z*(x - z - 1);
435 out[11][0][2] = 0.5*(out[11][0][0] + out[11][0][1])
436 - ((y + z + 1)*(x - z - 1) + 4*z) - z*(x - y - 2*z + 2);
437
438 out[12][0][0] = 2*z*(y + z + 1);
439 out[12][0][1] = 2*z*(x - z + 1);
440 out[12][0][2] = 0.5*(out[12][0][0] + out[12][0][1])
441 + (x - z + 1)*(y + z + 1) + z*(x - y - 2*z);
442
443 // base face
444 out[13][0][0] = 2*((y + z - 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
445 + (x - z + 1)*(y + z - 1)*(y + 1 - z));
446 out[13][0][1] = 2*((x - z + 1)*((y + 1)*(x - 1) - z*(x - y - z - 1))
447 + (x - z + 1)*(y + z - 1)*(x - 1 + z));
448 out[13][0][2] = 0.5*(out[13][0][0] + out[13][0][1])
449 + (x - y - 2*z + 2)*((y + 1)*(x - 1) - z*(x - y - z - 1))
450 + (x - z + 1)*(y + z - 1)*(-(x - y - 2*z - 1));
451 }
452
453 return;
454 }
455
456 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::evaluateJacobian for order " << order());
457 }
458
465 constexpr void partial(const std::array<unsigned int,3>& partialOrders,
466 const typename Traits::DomainType& in,
467 std::vector<typename Traits::RangeType>& out) const
468 {
469 auto totalOrder = std::accumulate(partialOrders.begin(), partialOrders.end(), 0);
470
471 out.resize(size());
472
473 if (totalOrder == 0)
474 {
475 evaluateFunction(in, out);
476 return;
477 }
478
479 if (order()==0)
480 {
481 out[0] = 0;
482 return;
483 }
484
485 if (order()==1)
486 {
487 if (totalOrder == 1)
488 {
489 auto const direction = std::distance(partialOrders.begin(), std::find(partialOrders.begin(), partialOrders.end(), 1));
490 if (in[0] > in[1])
491 {
492 switch (direction)
493 {
494 case 0:
495 out[0] = -1 + in[1];
496 out[1] = 1 - in[1];
497 out[2] = -in[1];
498 out[3] = in[1];
499 out[4] = 0;
500 break;
501 case 1:
502 out[0] = -1 + in[0] + in[2];
503 out[1] = -in[0] - in[2];
504 out[2] = 1 - in[0] - in[2];
505 out[3] = in[0]+in[2];
506 out[4] = 0;
507 break;
508 case 2:
509 out[0] = -1 + in[1];
510 out[1] = -in[1];
511 out[2] = -in[1];
512 out[3] = in[1];
513 out[4] = 1;
514 break;
515 default:
516 DUNE_THROW(RangeError, "Component out of range.");
517 }
518 }
519 else /* (in[0] <= in[1]) */
520 {
521 switch (direction)
522 {
523 case 0:
524 out[0] = -1 + in[1] + in[2];
525 out[1] = 1 - in[1] - in[2];
526 out[2] = -in[1] - in[2];
527 out[3] = in[1] + in[2];
528 out[4] = 0;
529 break;
530 case 1:
531 out[0] = -1 + in[0];
532 out[1] = -in[0];
533 out[2] = 1 - in[0];
534 out[3] = in[0];
535 out[4] = 0;
536 break;
537 case 2:
538 out[0] = -1 + in[0];
539 out[1] = -in[0];
540 out[2] = -in[0];
541 out[3] = in[0];
542 out[4] = 1;
543 break;
544 default:
545 DUNE_THROW(RangeError, "Component out of range.");
546 }
547 }
548 } else if (totalOrder == 2)
549 {
550 if ((partialOrders[0] == 1 && partialOrders[1] == 1) ||
551 (partialOrders[1] == 1 && partialOrders[2] == 1 && in[0] > in[1]) ||
552 (partialOrders[0] == 1 && partialOrders[2] == 1 && in[0] <=in[1]))
553 {
554 out = {1, -1, -1, 1, 0};
555 } else
556 {
557 out = {0, 0, 0, 0, 0};
558 }
559
560 } else
561 {
562 out = {0, 0, 0, 0, 0};
563 }
564
565 return;
566 }
567
568 if (order()==2)
569 {
570 if (totalOrder == 1)
571 {
572 // transform to reference element with base [-1,1]^2
573 const R x = 2.0*in[0] + in[2] - 1.0;
574 const R y = 2.0*in[1] + in[2] - 1.0;
575 const R z = in[2];
576
577 auto const direction = std::distance(partialOrders.begin(), std::find(partialOrders.begin(), partialOrders.end(), 1));
578
579 // transformation of the gradient leads to a multiplication
580 // with the Jacobian [2 0 0; 0 2 0; 1 1 1]
581 if (x > y)
582 {
583 switch (direction)
584 {
585 case 0:
586 out[0] = 0.5*(y - z - 1)*(y - z)*(2*x + 2*z - 1);
587 out[1] = 2*(-0.25*((y - z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-y + z + 1)) - z);
588 out[2] = 0.5*(y - z)*(y - z + 1)*(2*x + 2*z - 1);
589 out[3] = 0.5*(y - z)*(2*x + 2*z + 1)*(y - z + 1);
590 out[4] = 0;
591 out[5] = -(y - z + 1)*(2*x + z - 1)*(y - 1);
592 out[6] = -(y - z + 1)*(2*x + z + 1)*(y - 1);
593 out[7] = -(((y - z - 1)*(x + 1)*y - z) + z*(2*x + 1) + (x + z - 1)*((y - z - 1)*y + 2*z));
594 out[8] = -(y - z + 1)*(2*x + z)*y;
595 out[9] = 2*z*(y - z - 1);
596 out[10] = -2*z*(y - z - 1);
597 out[11] = -2*z*(y - z + 1);
598 out[12] = 2*z*(y - z + 1);
599 out[13] = 2*((y - z + 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(y - 1 + z));
600 break;
601 case 1:
602 out[0] = 0.5*(x + z)*(x + z - 1)*(2*y - 2*z - 1);
603 out[1] = 2*(-0.25*((x + z)*((x + z + 1)*(-y + z + 1) - 4*z) + (x + z)*(y - z)*(-(x + z + 1))) + z);
604 out[2] = 0.5*(x + z)*(2*y - 2*z + 1)*(x + z - 1);
605 out[3] = 0.5*(2*y - 2*z + 1)*(x + z)*(x + z + 1);
606 out[4] = 0;
607 out[5] = -(x + z - 1)*(y - 1)*x - (y - z + 1)*(x + z - 1)*x;
608 out[6] = -(((x + z + 1)*(y - 1)*x - z) + z*(2*y + 1) + (y - z + 1)*((x + z + 1)*x + 2*z));
609 out[7] = -(x + z - 1)*(2*y - z - 1)*(x + 1);
610 out[8] = -(2*y - z + 1)*(x + z - 1)*(x + 1);
611 out[9] = 2*z*(x + z - 1);
612 out[10] = -2*z*(x + z + 1);
613 out[11] = -2*z*(x + z - 1);
614 out[12] = 2*z*(x + z + 1);
615 out[13] = 2*((x + z - 1)*((y - 1)*(x + 1) + z*(x - y + z + 1)) + (y - z + 1)*(x + z - 1)*(x + 1 - z));
616 break;
617 case 2:
618 out[0] = -((y - z)*(2*x + 2*z - 1)*(z - y + 1))/2;
619 out[1] = ((y - z + 1)*(y - 2*x + z + 2*x*y - 2*x*z + 2*y*z - 2*z*z))/2;
620 out[2] = ((y - z)*(2*x + 2*z - 1)*(y - z + 1))/2;
621 out[3] = ((y - z)*(2*x + 2*z + 1)*(y - z + 1))/2;
622 out[4] = 4*z - 1;
623 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;
624 out[6] = -((y - z + 1)*(3*y - 2*x + z + 3*x*y + x*z + y*z + x*x - 1))/2;
625 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);
626 out[8] = -((y - z + 1)*(y + z + 3*x*y + x*z + y*z + x*x - 1))/2;
627 out[9] = -(x + 3*z - 1)*(z - y + 1);
628 out[10] = (x + z + 1)*(z - y + 1) - 2*y*z - 6*z + 2*z*z;
629 out[11] = -(x + 3*z - 1)*(y - z + 1);
630 out[12] = (x + 3*z + 1)*(y - z + 1);
631 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);
632 break;
633 default:
634 DUNE_THROW(RangeError, "Component out of range.");
635 }
636 }
637 else // x <= y
638 {
639 switch (direction)
640 {
641 case 0:
642 out[0] = -((y + z)*(2*z - 2*x + 1)*(y + z - 1))/2;
643 out[1] = ((y + z)*(2*x - 2*z + 1)*(y + z - 1))/2;
644 out[2] = -((y + z + 1)*(y - 3*z - 2*x*y - 2*x*z + 2*y*z + 2*z*z))/2;
645 out[3] = ((y + z)*(2*x - 2*z + 1)*(y + z + 1))/2;
646 out[4] = 0;
647 out[5] = (y + 1)*(y + z - 1)*(z - 2*x + 1);
648 out[6] = -(y + 1)*(2*x - z + 1)*(y + z - 1);
649 out[7] = -y*(2*x - z)*(y + z - 1);
650 out[8] = z - z*(2*x + 1) - (2*z + y*(y + z + 1))*(x - z + 1) - y*(x - 1)*(y + z + 1);
651 out[9] = 2*z*(y + z - 1);
652 out[10] = -2*z*(y + z - 1);
653 out[11] = -2*z*(y + z + 1);
654 out[12] = 2*z*(y + z + 1);
655 out[13] = 2*(y + z - 1)*(2*x - z + 2*x*y - 2*x*z + 2*z*z);
656 break;
657 case 1:
658 out[0] = -(x - z)*(y + z - 0.5)*(z - x + 1);
659 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
660 out[2] = -((z - x + 1)*(x + 3*z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
661 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
662 out[4] = 0;
663 out[5] = z - z*(2*y + 1) - (2*z - x*(z - x + 1))*(y + z - 1) + x*(y + 1)*(z - x + 1);
664 out[6] = -x*(2*y + z)*(x - z + 1);
665 out[7] = -(x - 1)*(x - z + 1)*(2*y + z - 1);
666 out[8] = -(x - 1)*(x - z + 1)*(2*y + z + 1);
667 out[9] = -2*z*(z - x + 1);
668 out[10] = -2*z*(x - z + 1);
669 out[11] = 2*z*(z - x + 1);
670 out[12] = 2*z*(x - z + 1);
671 out[13] = 2*(x - z + 1)*(2*x*y - z - 2*y + 2*y*z + 2*z*z);
672 break;
673 case 2:
674 out[0] = -((x - z)*(2*y + 2*z - 1)*(z - x + 1))/2;
675 out[1] = ((x - z)*(2*y + 2*z - 1)*(x - z + 1))/2;
676 out[2] = ((x - z + 1)*(x - 2*y + z + 2*x*y + 2*x*z - 2*y*z - 2*z*z))/2;
677 out[3] = ((x - z)*(2*y + 2*z + 1)*(x - z + 1))/2;
678 out[4] = 4*z - 1;
679 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);
680 out[6] = -((x - z + 1)*(x + z + 3*x*y + x*z + y*z + y*y - 1))/2;
681 out[7] = -((x - z + 1)*(3*x*y - 4*y - z - x + x*z + y*z + y*y + 1))/2;
682 out[8] = -((x - z + 1)*(3*x - 2*y + z + 3*x*y + x*z + y*z + y*y - 1))/2;
683 out[9] = -(z - x + 1)*(y + 3*z - 1);
684 out[10] = -(x - z + 1)*(y + 3*z - 1);
685 out[11] = (y + z + 1)*(z - x + 1) - 2*x*z - 6*z + 2*z*z;
686 out[12] = (x - z + 1)*(y + 3*z + 1);
687 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);
688 break;
689 default:
690 DUNE_THROW(RangeError, "Component out of range.");
691 }
692 }
693 } else {
694 DUNE_THROW(NotImplemented, "Desired derivative order is not implemented");
695 }
696
697 return;
698 }
699
700 DUNE_THROW(NotImplemented, "LagrangePyramidLocalBasis::partial for order " << OrderTraits::order());
701 }
702
703 };
704
709 template<int compileTimeOrder>
710 class LagrangePyramidLocalCoefficients
711 : private LagrangePyramidOrderTraits<compileTimeOrder>
712 {
713 using OrderTraits = LagrangePyramidOrderTraits<compileTimeOrder>;
714
715 public:
716
717 using OrderTraits::order;
718 using OrderTraits::size;
719
721 LagrangePyramidLocalCoefficients () requires (OrderTraits::is_static_order)
722 : LagrangePyramidLocalCoefficients(OrderTraits())
723 {}
724
726 explicit LagrangePyramidLocalCoefficients (OrderTraits orderTraits)
727 : OrderTraits(orderTraits)
728 , localKeys_(size())
729 {
730 if (order()==0)
731 {
732 localKeys_[0] = LocalKey(0,0,0);
733 return;
734 }
735
736 if (order()==1)
737 {
738 for (std::size_t i=0; i<size(); i++)
739 localKeys_[i] = LocalKey(i,3,0);
740 return;
741 }
742
743 if (order()==2)
744 {
745 // Vertex shape functions
746 localKeys_[0] = LocalKey(0,3,0);
747 localKeys_[1] = LocalKey(1,3,0);
748 localKeys_[2] = LocalKey(2,3,0);
749 localKeys_[3] = LocalKey(3,3,0);
750 localKeys_[4] = LocalKey(4,3,0);
751
752 // Edge shape functions
753 localKeys_[5] = LocalKey(0,2,0);
754 localKeys_[6] = LocalKey(1,2,0);
755 localKeys_[7] = LocalKey(2,2,0);
756 localKeys_[8] = LocalKey(3,2,0);
757 localKeys_[9] = LocalKey(4,2,0);
758 localKeys_[10] = LocalKey(5,2,0);
759 localKeys_[11] = LocalKey(6,2,0);
760 localKeys_[12] = LocalKey(7,2,0);
761
762 // base face shape function
763 localKeys_[13] = LocalKey(0,1,0);
764
765 return;
766 }
767
768 // No general case
769 DUNE_THROW(NotImplemented, "LagrangePyramidLocalCoefficients for order " << order());
770
771 }
772
774 const LocalKey& localKey (std::size_t i) const
775 {
776 return localKeys_[i];
777 }
778
779 private:
780 std::vector<LocalKey> localKeys_;
781 };
782
787 template<class D, class R, int compileTimeOrder>
788 class LagrangePyramidLocalInterpolation
789 : private LagrangePyramidOrderTraits<compileTimeOrder>
790 {
791 using OrderTraits = LagrangePyramidOrderTraits<compileTimeOrder>;
792 using Traits = typename LagrangePyramidLocalBasis<D,R,compileTimeOrder>::Traits;
793
794 public:
796 constexpr LagrangePyramidLocalInterpolation() requires (OrderTraits::is_static_order)
797 : OrderTraits()
798 {}
799
801 explicit constexpr LagrangePyramidLocalInterpolation(OrderTraits orderTraits)
802 : OrderTraits(orderTraits)
803 {}
804
812 template<typename F, typename C>
813 void interpolate (const F& f, std::vector<C>& out) const
814 {
815 const unsigned int k = OrderTraits::order();
816 using Domain = typename Traits::DomainType;
817 using DomainField = typename Traits::DomainFieldType;
818
819 out.resize(OrderTraits::size());
820
821 // Specialization for zero-order case
822 if (k==0)
823 {
824 auto center = ReferenceElements<DomainField,3>::general(GeometryTypes::pyramid).position(0,0);
825 out[0] = f(center);
826 return;
827 }
828
829 // Specialization for first-order case
830 if (k==1)
831 {
832 for (unsigned int i=0; i<OrderTraits::size(); i++)
833 {
834 auto vertex = ReferenceElements<DomainField,3>::general(GeometryTypes::pyramid).position(i,3);
835 out[i] = f(vertex);
836 }
837 return;
838 }
839
840 // Specialization for second-order case
841 if (k==2)
842 {
843 out[0] = f( Domain( {0.0, 0.0, 0.0} ) );
844 out[1] = f( Domain( {1.0, 0.0, 0.0} ) );
845 out[2] = f( Domain( {0.0, 1.0, 0.0} ) );
846 out[3] = f( Domain( {1.0, 1.0, 0.0} ) );
847 out[4] = f( Domain( {0.0, 0.0, 1.0} ) );
848 out[5] = f( Domain( {0.0, 0.5, 0.0} ) );
849 out[6] = f( Domain( {1.0, 0.5, 0.0} ) );
850 out[7] = f( Domain( {0.5, 0.0, 0.0} ) );
851 out[8] = f( Domain( {0.5, 1.0, 0.0} ) );
852 out[9] = f( Domain( {0.0, 0.0, 0.5} ) );
853 out[10] = f( Domain( {0.5, 0.0, 0.5} ) );
854 out[11] = f( Domain( {0.0, 0.5, 0.5} ) );
855 out[12] = f( Domain( {0.5, 0.5, 0.5} ) );
856 out[13] = f( Domain( {0.5, 0.5, 0.0} ) );
857
858 return;
859 }
860
861 DUNE_THROW(NotImplemented, "LagrangePyramidLocalInterpolation not implemented for order " << k);
862 }
863
864 };
865
866} } // namespace Dune::Impl
867
868namespace Dune
869{
900 template<class D, class R, int compileTimeOrder = -1>
901 class LagrangePyramidLocalFiniteElement
902 : private Impl::LagrangePyramidOrderTraits<compileTimeOrder>
903 {
904 using OrderTraits = Impl::LagrangePyramidOrderTraits<compileTimeOrder>;
905
906 public:
909 using Traits = LocalFiniteElementTraits<Impl::LagrangePyramidLocalBasis<D,R,compileTimeOrder>,
910 Impl::LagrangePyramidLocalCoefficients<compileTimeOrder>,
911 Impl::LagrangePyramidLocalInterpolation<D,R,compileTimeOrder>>;
912
914 constexpr LagrangePyramidLocalFiniteElement ()
915 : OrderTraits()
916 , basis_(*this)
917 , coefficients_(*this)
918 , interpolation_(*this)
919 {
920 static_assert(OrderTraits::is_static_order, "Default constructor only allowed for compile-time order >= 0");
921 static_assert((0 <= compileTimeOrder) and (compileTimeOrder <= 2), "LagrangePyramid: Only order 0,1, and 2 are supported");
922 }
923
925 explicit constexpr LagrangePyramidLocalFiniteElement (int runTimeOrder)
926 : OrderTraits(runTimeOrder)
927 , basis_(*this)
928 , coefficients_(*this)
929 , interpolation_(*this)
930 {
931 if ((runTimeOrder < 0) or (runTimeOrder > 2))
932 DUNE_THROW(Dune::InvalidStateException, "LagrangePyramid: Only run-time order 0,1, and 2 are supported");
933 if (compileTimeOrder >= 0 && compileTimeOrder != runTimeOrder)
934 DUNE_THROW(Dune::InvalidStateException, "LagrangePyramid: Compile-time order must be identical to run-time order!");
935 }
936
939 constexpr const typename Traits::LocalBasisType& localBasis () const
940 {
941 return basis_;
942 }
943
946 constexpr const typename Traits::LocalCoefficientsType& localCoefficients () const
947 {
948 return coefficients_;
949 }
950
953 constexpr const typename Traits::LocalInterpolationType& localInterpolation () const
954 {
955 return interpolation_;
956 }
957
959 using OrderTraits::size;
960
963 static constexpr GeometryType type ()
964 {
965 return GeometryTypes::pyramid;
966 }
967
968 private:
969 typename Traits::LocalBasisType basis_;
970 typename Traits::LocalCoefficientsType coefficients_;
971 typename Traits::LocalInterpolationType interpolation_;
972 };
973
974} // namespace Dune
975
976#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::InvalidStateException
Default exception if a function was called while the object is not in a valid state for that function...
Definition: exceptions.hh:375
Dune::LagrangePyramidLocalFiniteElement
Lagrange finite element for 3d pyramids with compile-time polynomial order.
Definition: lagrangepyramid.hh:903
Dune::LagrangePyramidLocalFiniteElement::localBasis
constexpr const Traits::LocalBasisType & localBasis() const
Returns the local basis, i.e., the set of shape functions.
Definition: lagrangepyramid.hh:939
Dune::LagrangePyramidLocalFiniteElement::localInterpolation
constexpr const Traits::LocalInterpolationType & localInterpolation() const
Returns object that evaluates degrees of freedom.
Definition: lagrangepyramid.hh:953
Dune::LagrangePyramidLocalFiniteElement::localCoefficients
constexpr const Traits::LocalCoefficientsType & localCoefficients() const
Returns the assignment of the degrees of freedom to the element subentities.
Definition: lagrangepyramid.hh:946
Dune::LagrangePyramidLocalFiniteElement::type
static constexpr GeometryType type()
The reference element that the local finite element is defined on.
Definition: lagrangepyramid.hh:963
Dune::LagrangePyramidLocalFiniteElement::LagrangePyramidLocalFiniteElement
constexpr LagrangePyramidLocalFiniteElement(int runTimeOrder)
Constructor for run-time order.
Definition: lagrangepyramid.hh:925
Dune::LagrangePyramidLocalFiniteElement::LagrangePyramidLocalFiniteElement
constexpr LagrangePyramidLocalFiniteElement()
Constructor for compile-time order.
Definition: lagrangepyramid.hh:914
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,...)
Definition: exceptions.hh:314
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::PDELab::interpolate
void interpolate(const F &f, const GFS &gfs, XG &xg)
interpolation from a given grid function
Definition: interpolate.hh:177
Dune::Hybrid::accumulate
constexpr T accumulate(Range &&range, T value, F &&f)
Accumulate values.
Definition: hybridutilities.hh:284
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 & Uni Heidelberg  |  generated with Hugo v0.111.3 (May 19, 22:31, 2026)