Dune Core Modules (2.7.0)

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