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